arXiv:1504.05158v1 [] 20 Apr 2015

Multi-swarm PSO algorithm for the Quadratic
Assignment Problem: a massive parallel
implementation on the OpenCL platform
Piotr Szwed and Wojciech Chmiel
AGH University of Science and Technology
Al. Mickiewicza 30
30 059 Kraków, Poland
{pszwed@agh.edu.pl, wch@agh.edu.pl}
November 16, 2017
Abstract
This paper presents a multi-swarm PSO algorithm for the Quadratic
Assignment Problem (QAP) implemented on OpenCL platform. Our
work was motivated by results of time efficiency tests performed for singleswarm algorithm implementation that showed clearly that the benefits of
a parallel execution platform can be fully exploited, if the processed population is large. The described algorithm can be executed in two modes:
with independent swarms or with migration. We discuss the algorithm
construction, as well as we report results of tests performed on several
problem instances from the QAPLIB library. During the experiments the
algorithm was configured to process large populations. This allowed us
to collect statistical data related to values of goal function reached by
individual particles. We use them to demonstrate on two test cases that
although single particles seem to behave chaotically during the optimization process, when the whole population is analyzed, the probability that
a particle will select a near-optimal solution grows.
Keywords: QAP, PSO, OpenCL, GPU calculation, particle swarm
optimization, mulit-swarm, discrete optimization

1

Introduction

Quadratic Assignment Problem (QAP) [19, 5] is a well known combinatorial
problem that can be used as optimization model in many areas [3, 23, 14]. QAP
may be formulated as follows: given a set of n facilities and n locations, the goal
is to find an assignment of facilities to unique locations that minimizes the sum
of flows between facilities multiplied by distances between their locations. As the

1

problem is NP hard [30], it can be solved optimally for small problem instances.
For larger problems (n>30), several heuristic algorithm were proposed [34, 24, 6].
One of the discussed methods [26, 21] is the Particle Swarm Optimization
(PSO). It attempts to find an optimal problem solution by moving a population
of particles in the search space. Each particle is characterized by two features its
position and velocity. Depending on a method variation, particles may exchange
information on their positions and reached values of goal functions [8].
In our recent work [33] we have developed PSO algorithm for the Quadratic
Assignment Problem on OpenCL platform. The algorithm was capable of processing one swarm, in which particles shared information about the global best
solution to update their search directions. Following typical patterns for GPU
based calculations, the implementation was a combination of parallel tasks (kernels) executed on GPU orchestrated by sequential operations run on the host
(CPU). Such organization of computations involves inevitable overhead related
to data transfer between the host and the GPU device. The time efficiency
test reported in [33] showed clearly that the benefits of a parallel execution
platform can be fully exploited, if processed populations are large, e.g. if they
comprise several hundreds or thousands particles. For smaller populations sequential algorithm implementation was superior both as regards the total swarm
processing time and the time required to process one particle. This suggested
a natural improvement of the previously developed algorithm: by scaling it up
to high numbers of particles organized into several swarms.
In this paper we discuss a multi-swarm implementation PSO algorithm for
the QAP problem on OpenCL platform. The algorithm can be executed in two
modes: with independent swarms, each maintaining its best solution, or with
migration between swarms. We describe the algorithm construction, as well
as we report tests performed on several problem instances from the QAPLIB
library [28]. Their results show advantages of massive parallel computing: the
obtained solutions are very close to optimal or best known for particular problem
instances.
The developed algorithm is not designed to exploit the problem specificity
(see for example [10]), as well as it is not intended to compete with supercomputer or grid based implementations providing exact solutions for the QAP [2].
On the contrary, we are targeting low-end GPU devices, which are present in
most laptops and workstations in everyday use, and accept near-optimal solutions.
During the tests the algorithm was configured to process large numbers of
particles (in the order of 10000). This allowed us to collect data related to
goal function values reached by individual particles and present such statistical measures as percentile ranks and probability mass functions for the whole
populations or selected swarms.
The paper is organized as follows: next Section 2 discusses the QAP problem,
as well as the PSO method. It is followed by Section 3, which describes the
adaptation of PSO to the QAP and the parallel implementation on the OpenCL
platform. Experiments performed and their results are presented in Section 4.
Section 5 provides concluding remarks.
2

2

Related works

2.1

Quadratic Assignment Problem

Quadratic Assignment Problem was introduced by Koopmans and Beckman in
1957 as a mathematical model describing assignment of economic activities to
a set of locations [19].
Let V = {1, ..., n} be a set of locations (nodes) linked by n2 arcs. Each
arc linking a pair of nodes (k, l) is attributed with a non-negative weight dkl
interpreted as a distance. Distances are usually presented in form of n × n
distance matrix D = [dkl ]. The next problem component is a set of facilities
N = {1, ..., n} and a n × n non-negative flow matrix F = [fij ], whose elements
describe flows between pairs of facilities (i, j).
The problem goal is to find an assignment π : N → V that minimizes the
total cost calculated as sum of flows fij between pairs of facilities (i, j) multiplied
by distances dπ(i)π(j) between pairs of locations (π(i), π(j)), to which they are
assigned. The permutation π can be encoded as n2 binary variables xki , where
k = π(i), what gives the following problem statement:
min

n X
n X
n X
n
X

fij dkl xki xlj

(1)

i=1 j=1 k=1 l=1

subject to:
Pn

xij = 1,

Pn

for 1 ≤ j ≤ n

xij = 1,

for 1 ≤ i ≤ n

i=1
j=1

(2)

xij ∈ {0, 1}
The n × n matrix X = [xki ] satisfying (2) is called permutation matrix.
In most cases matrix D and F are symmetric. Moreover, their diagonal
elements are often equal 0. Otherwise, the component fii dkk xki xki can be extracted as a linear part of the goal function interpreted as an installation cost
of i-th facility at k-th location .
QAP models found application in various areas including transportation
[3], scheduling, electronics (wiring problem), distributed computing, statistical
data analysis (reconstruction of destroyed soundtracks), balancing of turbine
running [23], chemistry , genetics [29], creating the control panels and manufacturing [14].
In 1976 Sahni and Gonzalez proved that the QAP is strongly N P-hard [30],
by showing that a hypothetical existence of a polynomial time algorithm for
solving the QAP would imply an existence of a polynomial time algorithm for
an N P-complete decision problem - the Hamiltonian cycle.
In many research works QAP is considered one of the most challenging optimization problem. This in particular regards problem instances gathered in a
publicly available and continuously updated QAPLIB library [28, 4]. A practical size limit for problems that can be solved with exact algorithms is about

3

n = 30 [16]. In many cases optimal solutions were found with branch and bound
algorithm requiring high computational power offered by computational grids
[2] or supercomputing clusters equipped with a few dozen of processor cores and
hundreds gigabytes of memory [15]. On the other hand, in [10] a very successful
approach exploiting the problem structure was reported. It allowed to solve
several hard problems from QAPLIB using very little resources.
A number of heuristic algorithms allowing to find a near-optimal solutions
for QAP were proposed. They include Genetic Algorithm [1], various versions
of Tabu search [34], Ant Colonies [32, 12] and Bees algorithm [11]. Another
method, being discussed further, is Particle Swarm Optimization [26, 21] .

2.2

Particle Swarm Optimization

The classical PSO algorithm [8] is an optimization method defined for continuous
domain. During the optimization process a number of particles move through
a search space and update their state at discrete time steps t = 1, 2, 3, . . . Each
particle is characterized by position x(t) and velocity v(t). A particle remembers
its best position reached so far pL (t), as well as it can use information about
the best solution found by the swarm pG (t).
The state equation for a particle is given by the formula (3). Coefficients
c1 , c2 , c3 ∈ [0, 1] are called respectively inertia, cognition (or self recognition)
and social factors, whereas r1 , r2 are random numbers uniformly distributed in
[0, 1]
v(t + 1) = c1 · v(t) + c2 · r2 (t) · (pL (t) − x(t)) + c3 · r3 (t) · (pG (t) − x(t))

x(t + 1) = x(t) + v(t)

(3)

An adaptation of the PSO method to a discrete domain necessities in giving
interpretation to the velocity concept, as well as defining equivalents of scalar
multiplication, subtraction and addition for arguments being solutions and velocities. Examples of such interpretations can be found in [7] for the TSP and
[26] for the QAP.
A PSO algorithm for solving QAP using similar representations of particle state was proposed by Liu et al. [21]. Although the approach presented
there was inspiring, the paper gives very little information on efficiency of the
developed algorithm.

2.3

GPU based calculations

Recently many computationally demanding applications has been redesigned to
exploit the capabilities offered by massively parallel computing GPU platforms.
They include such tasks as: physically based simulations, signal processing, ray
tracing, geometric computing and data mining [27]. Several attempts have been
also made to develop various population based optimization algorithms on GPUs
including: the particle swarm optimization [36], the ant colony optimization [35],

4

the genetic [22] and memetic algorithm [20]. The described implementations
benefit from capabilities offered by GPUs by processing whole populations by
fast GPU cores running in parallel.

3

Algorithm design and implementation

In this section we describe the algorithm design, in particular the adaptation
of Particle Swarm Optimization metaheuristic to the QAP problem, as well as
a specific algorithm implementation on OpenCL platform. As it was stated in
Section 2.2, the PSO uses generic concepts of position x and velocity v that can
be mapped to a particular problem in various ways. Designing an algorithm for
a GPU platform requires decisions on how to divide it into parts that are either
executed sequentially at the host side or in parallel on the device.

3.1

PSO adaptation for the QAP problem

A state of a particle is a pair (X, V ). In the presented approach both are
n × n matrices, where n is the problem size. The permutation matrix X = [xij ]
encodes an assignment of facilities to locations. Its elements xij are equal to 1,
if j-th facility is assigned to i-th location, and take value 0 otherwise.
A particle moves in the solution space following the direction given by the
velocity V . Elements vij have the following interpretation: if vij has high
positive value, then a procedure determining the next solution should favor an
assignment xij = 1. On the other hand, if vij ≤ 0, then xij = 0 should be
preferred.
The state of a particle reached in the t-th iteration will be denoted by
(X(t), V (t)). In each iteration it is updated according to formulas (4) and
(5).
V (t + 1) = Sv c1 · V (t) + c2 · r2 (t) · (P L (t) − X(t))+

c3 · r3 (t) · (P G (t) − X(t))
X(t + 1) = Sx (X(t) + V (t))

(4)

(5)

Coefficients r2 and r3 are random numbers from [0, 1] generated in each iteration for every particle separately. They are introduced to model a random
choice between movements in the previous direction (according to c1 – inertia), the best local solution (self recognition) or the global best solution (social
behavior).
All operators appearing in (4) and (5) are standard operators from the linear
algebra. Instead of redefining them for a particular problem, see e.g. [7], we
propose to use aggregation functions Sv and Sx that allow to adapt the algorithm
to particular needs of a discrete problem.

5

The function Sv is used to assure that velocities have reasonable values.
Initially, we thought that unconstrained growth of velocity can be a problem,
therefore we have implemented a function, which restricts the elements of V to
an interval [−vmax , vmax ]. This function is referred as raw in Table 1. However,
the experiments conducted showed, that in case of small inertia factor, e.g. c1 =
0.5, after a few iterations all velocities tend to 0 and in consequence all particles
converge to the best solution encountered earlier by the swarm. To avoid such
effect another function that additionally performs column normalization was
proposed.
Pn For each j-th column the sum of absolute values of the elements
nj = i=1 |vij | is calculated and then the following assignment is made: vij ←
vij /nj .
According to formula (5) a new particle position X(t + 1) is obtained by
aggregating the previous state components: X(t) and V (t). As elements of a
matrix X(t) + V (t) may take values from [−vmax , vmax + 1], the Sx function
is responsible for converting it into a valid permutation matrix satisfying (2),
i.e. having exactly one 1 in each row and column. Actually, Sv is rather a
procedure, than a function, as it incorporates some elements of random choice.
Three variants of Sx procedures were implemented:
1. GlobalM ax(X) – iteratively searches for xrc , a maximum element in a
matrix X, sets it to 1 and clears other elements in the row r and c.
2. P ickColumn(X) – picks a column c from X, selects a maximum element
xrc , replaces it by 1 and clears other elements in r and c.
3. SecondT arget(X, Z, d) – similar to GlobalM ax(X), however during the
first d iterations ignores elements xij , such that zij = 1. (As the parameter
Z a solution X from the last iteration is used.)
In several experiments, in which GlobalM ax aggregation procedure was applied, particles seemed to get stuck, even if their velocities were far from zero
[33]. We reproduce this effect on a small 3 × 3 example:

1 0 0
X= 0 0 1 
0 1 0


8 1 3
X +V = 0 4 6 
2 4 2




7
V = 0
2

1
Sx (X + V ) =  0
0

1
4
3
0
0
1


3
5 
2

0
1 
0

If GlobalM ax(X(t) + V (t)) is used for the described case, in subsequent
iterations it will hold: X(t + 1) = X(t), until another particle is capable of
changing (P G (t) − X(t))) component of formula (4) for velocity calculation.
A solution for this problem can be to move a particle to a secondary direction
(target), by ignoring d < n elements that are in the solution X(t) already set
to 1. This, depending on d, gives an opportunity to reach other solutions,
hopefully yielding smaller goal function values. If chosen elements are maximal
6

in the remaining matrix, denoted here as X d V , they are still reasonable
movement directions. For the discussed example possible values of X d V
matrices, where d = 1 and d = 2, are shown below (6). Elements of a new
solution are marked with circles, whereas the upper index indicates the iteration,
in which the element was chosen.
82
 0
1V =
2

8
3

0
V
=
2
2


X

X



X

8

0
2V =
23

1
4
43
1
4
41
1
41
4


3
61 
2

32
6 
2
or


2

3
6 
2

(6)

It can be observed that for d = 1 the value is exactly the same, as it would
result from the GlobalM ax, however setting d = 2 allows to reach a different
solution. The pseudocode of SecondT arget procedure is listed in Algorithm 1.

3.2

Migration

The intended goal of the migration mechanism is to improve the algorithm
exploration capabilities by exchanging information between swarms. Actually,
it is not a true migration, as particles do not move. Instead we modify stored
P G [k] solutions (global best solution for a k-th swarm) replacing it by randomly
picked solution from a swarm that performed better (see Algortithm 2).
The newly set P G [k] value influences the velocity vector for all particles in
k-th swarm according to the formula (4). It may happen that the goal function
value corresponding to the assigned solution P G [k]) is worse than the previous
one. It is accepted, as the migration is primarily designed to increase diversity
within swarms.
It should be mentioned that a naive approach consisting in copying best
P G values between the swarms would be incorrect. (Consider replacing line 5
of Algorithm 2 with: P G [sm−k−1 ] ← P G [sk ].) In such case during algorithm
stagnation spanning over several iterations: in the first iteration the best value
P G [1] would be cloned, in the second two copies would be created, in the third
four and so on. Finally, after k iterations 2k swarms would follow the same
direction. In the first group of experiments reported in Section 4 we used up to
250 swarms. It means that after 8 iterations all swarms would be dominated by
a single solution.

7

Algorithm 1 Aggregation procedure SecondT arget
Require: X = Z + V - new solution requiring normalization
Require: Z - previous solution
Require: depth - number of iterations, in which during selection of maximum
element the algorithm ignore positions, where corresponding element of Z
is equal to 1
1: procedure SecondTarget(X,Z,depth)
2:
R ← {1, . . . , n}
3:
C ← {1, . . . , n}
4:
for i in (1,n) do
5:
Calculate M , the set of maximum elements
6:
if i ≤ depth then
7:
Ignore elements xij such that zij = 1
8:
M←
{(r, c) : zrc 6= 1 ∧ ∀i∈R,j∈C,zij 6=1 (xrc ≥ xij )}
9:
else
10:
M ← {(r, c) : ∀i∈R,j∈C (xrc ≥ xij )}
11:
end if
12:
Randomly select (r, c) from M
13:
R ← R \ {r}
. Update the sets R and C
14:
C ← C \ {c}
15:
for i in (1, n) do
16:
xri ← 0
. Clear r-th row
17:
xic ← 0
. Clear c-th column
18:
end for
19:
xrc ← 1
. Assign 1 to the maximum element
20:
end for
21:
return X
22: end procedure

8

Algorithm 2 Migration procedure
Require: d - migration depth satisfying d < m/2, where m is the number of
swarms
Require: P G - table of best solutions for m swarms
Require: X - set of all solutions
1: procedure migration(d,P G ,X)
2:
Sort swarms according to their PkG values into
a sequence (s1 , s2 . . . , sm−2 , sm−1 )
3:
for k in (1, d) do
4:
Randomly choose a solution xkj
belonging to the swarm sk
5:
Assign P G [sm−k−1 ] ← xkj
6:
Update the best goal function value for
the swarm m − k − 1
7:
end for
8: end procedure

3.3

OpenCL algorithm implementation

OpenCL [18] is a standard providing a common language, programming interfaces and hardware abstraction for heterogeneous platforms including GPU,
multicore CPU, DSP and FPGA [31]. It allows to accelerate computations by
decomposing them into a set of parallel tasks (work items) operating on separate
data.
A program on OpenCL platform is decomposed into two parts: sequential
executed by the CPU host and parallel executed by multicore devices. Functions
executed on devices are called kernels. They are written in a language being a
variant of C with some restrictions related to keywords and datatypes. When
first time loaded, the kernels are automatically translated into the instruction
set of the target device. The whole process takes about 500ms.
OpenCL supports 1D, 2D or 3D organization of data (arrays, matrices and
volumes). Each data element is identified by 1 to 3 indices, e.g. d[i][j] for twodimensional arrays. A work item is a scheduled kernel instance, which obtains
a combination of data indexes within the data range. To give an example, a 2D
array of data of n × m size should be processed by n · m kernel instances, which
are assigned with a pair of indexes (i, j), 0 ≤ i < n and 0 ≤ j < m. Those
indexes are used to identify data items assigned to kernels.
Additionally, kernels can be organized into workgroups, e.g. corresponding
to parts of a matrix, and synchronize their operations within a group using
so called local barrier mechanism. However, workgroups suffer from several
platform restrictions related to number of work items and amount of accessible
memory.
OpenCL uses three types of memory: global (that is exchanged between the
host and the device), local for a work group and private for a work item.
In our implementation we used aparapi platform [17] that allows to write
9

Generate particles

Apply Sx

Calculate goal function

Update best solutions

Generate velocities

YES

STOP
NO
Apply Sx

Calculate goal function

Update best solutions

Apply migration

Update velocities and apply Sv

Figure 1: Functional blocks of OpenCL based algorithm implementation.
OpenCL programs directly in Java language. The platform comprises two parts:
an API and a runtime capable of converting Java bytecodes into OpenCL workloads. Hence, the host part of the program is executed on a Java virtual machine, and originally written in Java kernels are executed on an OpenCL enabled
device.
The basic functional blocks of the algorithm are presented in Fig. 1. Implemented kernels are marked with gray color. The code responsible for generation
of random particles is executed by the host. We have also decided to leave
the code for updating best solutions at the host side. Actually, it comprises a
number of native System.arraycopy() calls. This regards also the migration
procedure, which sorts swarms indexes in a table according to their P G value
and copies randomly picked entries.
Particle data comprise a number of matrices (see Fig. 2): X and Xnew –
solutions, P L – local best particle solution and V – velocity. They are all
stored in large continuous tables shared by all particles. Appropriate table part
belonging to a particle can be identified based on the particle id transferred to
a kernel. Moreover, while updating velocity, the particles reference a table P G
indexed by the swarm id.
An important decision related to OpenCL program design is related to data
10

#p

#p

X

#p

V

PL
#p

#p

#sw

goal function
values
Xnew

PG

rands

Figure 2: Global variables used in the algorithm implementation
ranges selection. The memory layout in Fig. 2 suggests 3D range, whose dimensions are: row, column and particle number. This can be applied for relatively
simple velocity or goal function calculation. However, the proposed algorithms
for Sx , see Algorithm 1, are far too complicated to be implemented as a simple
parallel work item. Finally we decided to use one dimension (particle id) for Sx
and goal function calculation, and two dimensions (particle id, swarm id) for
velocity kernels.
It should be mentioned that requirements of parallel processing limits applicability of object oriented design at the host side. In particular we avoid creating
particles or swarms with their own memory and then copying small chunks between the host and the device. Instead we rather use a flyweight design pattern
[13]: if a particle abstraction is needed, a single object can be configured to
see parts of large global arrays X, V as its own memory and perform required
operations, e.g. initialization with random values.

4

Experiments and results

In this section we report results of conducted experiments, which aimed at
establishing the optimization performance of the implemented algorithm, as
well as to collect data related to its statistical properties.

4.1

Optimization results

The algorithm was tested on several problem instances form the QAPLIB [28],
whose size ranged between 12 and 150. Their results are gathered in Table 1
and Table 2. The selection of algorithm configuration parameters (c1 , c2 and
c3 factors, as well as the kernels used) was based on previous results published
in [33]. In all cases the second target Sx aggregation kernel was applied (see
Algorithm 1), which in previous experiments occurred the most successful.
During all tests reported in Table 1, apart the last, the total numbers of
particles were large: 10000-12500. For the last case only 2500 particles were

11

used due to 1GB memory limit of the GPU device (AMD Radeon HD 6750M
card). In this case the consumed GPU memory ranged about 950 MB.
The results show that algorithm is capable of finding solutions with goal
function values are close to reference numbers listed in QAPLIB. The gap is
between 0% and 6.4% for the biggest case tai150b. We have repeated tests
for tai60b problem to compare the implemented multi-swarm algorithm with
the previous single-swarm version published in [33]. Gap values for the best
results obtained with the single swarm algorithm were around 7%-8%. For the
multi-swarm implementation discussed here the gaps were between 0.64% and
2.03%.
The goal of the second group of experiments was to test the algorithm configured to employ large numbers of particles (50000-10000) for well known esc32*
problem instances from the QAPLIB. Altough they were considered hard, all of
them have been were recently solved optimally with exact algorithms [25, 10].
The results are summarized in Table 2. We used the following parameters:
c1 = 0.8, c2 = 0.5 c3 = 0.5, velocity kernel: normalized, Sx kernel: second target. During nearly all experiments optimal values of goal functions were reached
in one algorithm run. Only the problem esc32a occurred difficult, therefore for
this case the number of particles, as well as the upper iteration limits were
increased to reach the optimal solution. What was somehow surprising, in all
cases solutions differing from those listed in QAPLIB were obtained. Unfortunately, our algorithm was not prepared to collect sets of optimal solutions, so
we are not able to provide detailed results on their numbers.
It can be seen that optimal solutions for problem instances esc32c–h were
found in relatively small numbers of iterations. In particular, for esc32e and
es32g, which are characterized by small values of goal functions, optimal solutions were found during the initialization or in the first iteration.
The disadvantage of the presented algorithm is that it uses internally matrix representation for solutions and velocities. In consequence the memory
consumption is proportional to n2 , where n is the problem size. The same regards the time complexity, which for goal function and Sx procedures can be
estimated as o(n3 ). This makes optimization of large problems time consuming
(e.g. even 400 sec for one iteration for tai150b). However, for for medium size
problem instances, the iteration times are much smaller, in spite of large populations used. For two runs of the algorithm bur26a reported in Table 1, where
during each iteration 12500 particles were processed, the average iteration time
was equal 1.73 sec. For 50000-10000 particles and problems of size n = 32 the
average iteration time reported in Table 2 was less than 3.7 seconds.

4.2

Statistical results

An obvious benefit of massive parallel computations is the capability of processing large populations (see Table 2). Such approach to optimization may
resemble a little bit a brutal force attack: the solution space is randomly sampled millions of times to hit the best solution. No doubt that such approach
can be more successful if combined with a correctly designed exploration mech12

Instance

chr12a
bur26a
bur26a
lipa50a
tai60a
tai60a
tai60b
tai60b
tai60b
tai60b
tai64c
esc64a
tai80a
tai80b
sko100a
tai100b
esc128
tai150b

No

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

Self recognition c2

Inertia c1

Total particles

Swarm size

Number of swarms

Size
0.5
0.5
0.5
0.5
0.5
0.5
0.3
0.5
0.3
0.3
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

Social factor c3

0.5
0.5
0.5
0.5
0.5
0.5
0.3
0.5
0.3
0.3
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

Norm
Raw
Raw
Norm
Norm
Raw
Norm
Norm
Norm
Norm
Norm
Raw
Raw
Raw
Raw
Raw
Raw
Raw

Velocity kernel

0.5
0.8
0.8
0.5
0.8
0.8
0.8
0.5
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8

0
33%
0
0
0
0
0
0
33%
33%
0
0
0
0
0
0
0
0

Migration factor

10000
12500
12500
10000
10000
15000
10000
10000
10000
10000
10000
10000
10000
10000
10000
10000
5000
2500

9 552
5 426 670
5 429 693
62 794
7 539 614
7 426 672
620 557 952
617 825 984
612 078 720
614 088 768
1 856 396
116
14 038 392
835 426 944
154 874
1 196 819 712
64
530 816 224

Reached goal

50
50
50
50
50
50
50
50
50
50
50
50
50
50
50
50
50
50

9 552
5 426 670
5 426 670
62 093
7 205 962
7 205 962
608 215 054
608 215 054
608 215 054
608 215 054
1 855 928
116
13 499 184
818 415 043
152 002
1 185 996 137
64
498 896 643

Reference value

200
250
250
200
200
300
200
200
200
200
200
200
200
200
200
200
100
50

0.00%
0.00%
0.06%
1.13%
4.63%
3.06%
2.03%
1.58%
0.64%
0.97%
0.03%
0.00%
3.99%
2.08%
1.89%
0.91%
0.00%
6.40%

Gap

12
26
26
50
60
60
60
60
60
60
64
64
80
80
100
100
128
150

Table 1: Results of tests for various instance from QAPLIB

21
156
189
1640
817
917
909
1982
2220
1619
228
71
1718
1509
1877
1980
1875
1894

Iteration

13

Table 2: Results of tests for esc32* instances from QAPLIB (problem size
n = 32). Reached optimal values are marked with asterisks.
Instance
esc32a
esc32a
esc32a
esc32b
esc32c
esc32d
esc32e
esc32g
esc32h

Swarms
50
10
50
50
50
50
50
50
50

Particles
1000
5000
2000
1000
1000
1000
1000
1000
1000

Total part.
100000
100000
100000
50000
50000
50000
50000
50000
50000

Goal
138
134
130∗
168∗
642∗
400∗
2∗
6∗
438∗

Iter
412
909
2407
684
22
75
0
1
77

Time/iter [ms]
3590.08
3636.76
3653.88
3637.84
3695.19
3675.32
3670.38
3625.17
3625.17

anism that directs the random search process towards solutions providing good
or near-optimal solutions. In this section we analyze collected statistical data
related to the algorithm execution to show that the optimization performance of
the algorithm can be attributed not only to large sizes of processed population,
but also to the implemented exploration mechanism.
PSO algorithm can be considered a stochastic process controlled by random
variables r2 (t) and r3 (t) appearing in its state equation (3). Such analysis for
continuous problems were conducted in [9]. On the other hand, the observable
algorithm outcomes, i.e. the values of goal functions f (xi (t)) for solutions xi ,
i = 1, . . . , n reached in consecutive time moments t ∈ {1, 2, 3, . . . } can be also
treated as random variables, whose distributions change over time t. Our intuition is that a correctly designed algorithm should result in a nonstationary
stochastic process {f (xi (t)) : t ∈ T }, characterized by growing probability that
next values of goal functions in the analyzed population are closer to the optimal
solution.
To demonstrate such behavior of the implemented algorithm we have collected detailed information on goal function values during two optimization task
for the problem instance bur26a reported in Table 1 (cases 2 and 3). For both of
them the algorithm was configured to use 250 swarms comprising 50 particles.
In the case 2 the migration mechanism was applied and the optimal solution
was found in the iteration 156, in the case 3 (without migration) a solution very
close to optimal (gap 0.06%) was reached in the iteration 189.
Fig. 3 shows values of goal function for two selected particles during run
3. The plots show typical QAP specificity. PSO and many other algorithms
perform a local neighborhood search. For the QAP the neighborhood is characterized by great variations of goal function values. Although mean values of
goal function decrease in first twenty or thirty iterations, the particles behave
randomly and nothing indicates that during subsequent iterations smaller values
of goal functions would be reached more often.
In Fig. 4 percentile ranks (75%, 50% 25% and 5%) for two swarms, which

14

6.2

·106
#4716
#643

goal

6
5.8
5.6
5.4
−20

0

20

40

60

80

100 120 140

160 180 200

220

iteration

Figure 3: Variations of goal function values for two particles exploring the solutions space during the optimization process (bur26a problem instance)
reached best values in cases 2 and 3 are presented. Although the case 3 is
characterized by less frequent changes of scores, than the case 2, probably this
effect can not be attributed to the migration applied. It should be mentioned
that for a swarm comprising 50 particles, the 0.05 percentile corresponds to just
two of them.
·106

goal

6

Best

pct-0.75
pct-0.50
pct-0.25
pct-0.05

5.8

5.6

5.4
−20

0

20

40

60

80

100 120 140

160 180 200

220

iteration

·106

goal

Min (opt)

pct-0.75
pct-0.50
pct-0.25
pct-0.05

6
5.8
5.6
5.4
−20

0

20

40

60

80

100 120 140

160 180 200

220

iteration

Figure 4: Two runs of bur26a optimization. Percentile ranks for 50 particles
belonging to the most successful swarms: without migration (above) and with
migration (below).
Collected percentile rank values for the whole population comprising 12500

15

particles are presented in Fig. 5. For both cases the plots are clearly separated. It can be also observed that solutions very close to optimal are practically reached between the iterations 20 (37.3 sec) and 40 (72.4 sec). For the
whole population the 0.05 percentile represents 625 particles. Starting with the
iteration 42 their score varies between 5.449048 · 106 and 5.432361 · 106 , i.e. by
about 0.3%.
·106

goal

6

Best

pct-0.75
pct-0.50
pct-0.25
pct-0.05

5.8

5.6

5.4
−20

0

20

40

60

80

100 120 140

·106

220

Min (opt)

6

goal

160 180 200

iteration

pct-0.75
pct-0.50
pct-0.25
pct-0.05

5.8

5.6

5.4
−20

0

20

40

60

80

100 120 140

160 180 200

220

iteration

Figure 5: Two runs of bur26a optimization. Percentile ranks for all 12500
particles: without migration (above) and with migration (below).
Fig. 6 shows, how the probability distribution (probability mass function
PMF) changed during the optimization process. In both cases the the optimization process starts with a normal distribution with the mean value about
594500. In the subsequent iterations the maximum of PMF grows and moves
towards smaller values of the goal function. There is no fundamental difference
between the two cases, however for the case 3 (with migration) maximal values
of PMF are higher. It can be also observed that in the iteration 30 (completed
in 56 seconds) the probability of hitting a good solution is quite high, more then
10%.
Interpretation of PMF for the two most successful swarms that reached best
values in the discussed cases is not that obvious. For the case without migration
(Fig. 7 above) there is a clear separation between the initial distribution and
the distribution reached in the iteration, which yielded the best result. In the
second case (with migration) a number of particles were concentrated around
local minima.
The presented data shows advantages of optimization performed on mas16

0.15
it#1
it#30
it#60
it#120
it#189

p(f )

0.1

5 · 10−2
0
5.4

5.6

5.8

6

6.2

6.4
·106

goal value f

p(f )

0.2
it#1
it#30
it#50
it#100
it#156

0.1

0
5.4

5.6

5.8

6
goal value f

6.2

6.4
·106

Figure 6: Probability mass functions for 12050 particles organized into 250 x
50 swarms during two runs: without migration (above) and with migration
(below).
sive parallel processing platforms. Due to high number of solutions analyzed
simultaneously, the algorithm that does not exploit the problem structure can
yield acceptable results in relatively small number of iterations (and time). For
a low-end GPU devices, which was used during the test, good enough results
were obtained after 56 seconds. It should be mentioned that for both presented
cases the maximum number of iterations was set to 200. With 12500 particles,
the ratio of potentially explored solutions to the whole solution space was equal
200 · 12500/26! = 6.2 · 10−21 .

5

Conclusions

In this paper we describe a multi-swarm PSO algorithm for solving the QAP
problem designed for the OpenCL platform. The algorithm is capable of processing in parallel large number of particles organized into several swarms that
either run independently or communicate with use of the migration mechanism.
Several solutions related to particle state representation and particle movement
were inspired by the work of Liu at al. [21], however, they were refined here to
provide better performance.
We tested the algorithm on several problem instances from the QAPLIB library obtaining good results (small gaps between reached solutions and reference
values). However, it seems that for problem instances of large sizes the selected
17

p(f )

0.6
it#1
it#189

0.4
0.2
0
5.4

5.5

5.6

5.7

5.8

5.9

6

6.1

6.2

6.3

6.4
·106

6

6.1

6.2

6.3

6.4
·106

goal value f

it#1
it#156

p(f )

0.2

0.1

0

5.4

5.5

5.6

5.7

5.8

5.9

goal value f

Figure 7: Probability mass functions for 50 particles belonging to the most successful swarms during two runs: without migration (above) and with migration
(below). One point represents an upper bound for 5 particles.
representation of solutions in form of permutation matrices hinders potential
benefits of parallel processing.
During the experiments the algorithm was configured to process large populations. This allowed us to collect statistical data related to goal function
values reached by individual particles. We used them to demonstrate on two
cases that although single particles seem to behave chaotically during the optimization process, when the whole population is analyzed, the probability that
a particle will select a near-optimal solution grows. This growth is significant
for a number of initial iterations, then its speed diminishes and finally reaches
zero.
Statistical analysis of experimental data collected during optimization process may help to tune the algorithm parameters, as well as to establish realistic
limits related to expected improvement of goal functions. This in particular
regards practical applications of optimization techniques, in which recurring
optimization problems appear, i.e. the problems with similar size, complexity
and structure. Such problems can be near-optimally solved in bounded time on
massive parallel computation platforms even, if low-end devices are used.

18

References
[1] Ahuja, R.K., Orlin, J.B., Tiwari, A.: A greedy genetic algorithm for the
quadratic assignment problem. Computers & Operations Research 27(10),
917–934 (2000)
[2] Anstreicher, K., Brixius, N., Goux, J.P., Linderoth, J.: Solving large
quadratic assignment problems on computational grids. Mathematical Programming 91(3), 563–588 (2002)
[3] Bermudez, R., Cole, M.H.: A genetic algorithm approach to door assignments in breakbulk terminals. Tech. Rep. MBTC-1102, Mack-Blackwell
Transportation Center, University of Arkansas, Fayetteville, Arkansas
(2001)
[4] Burkard, R.E., Karisch, S.E., Rendl, F.: QAPLIB - a Quadratic Assignment Problem library. Journal of Global Optimization 10(4), 391–403
(1997)
[5] Çela, E.: The quadratic assignment problem: theory and algorithms. Combinatorial Optimization, Springer, Boston (1998)
[6] Chmiel, W., Kadłuczka, P., Packanik, G.: Performance of swarm algorithms for permutation problems. Automatyka 15(2), 117–126 (2009)
[7] Clerc, M.: Discrete particle swarm optimization, illustrated by the traveling
salesman problem. In: New optimization techniques in engineering, pp.
219–239. Springer (2004)
[8] Eberhart, R., Kennedy, J.: A new optimizer using particle swarm theory.
In: Micro Machine and Human Science, 1995. MHS ’95., Proceedings of
the Sixth International Symposium on. pp. 39–43 (Oct 1995)
[9] Fernández Martínez, J., García Gonzalo, E.: The PSO family: deduction, stochastic analysis and comparison. Swarm Intelligence 3(4), 245–273
(2009)
[10] Fischetti, M., Monaci, M., Salvagnin, D.: Three ideas for the quadratic
assignment problem. Operations Research 60(4), 954–964 (2012)
[11] Fon, C.W., Wong, K.Y.: Investigating the performance of bees algorithm
in solving quadratic assignment problems. International Journal of Operational Research 9(3), 241–257 (2010)
[12] Gambardella, L.M., Taillard, E., Dorigo, M.: Ant colonies for the quadratic
assignment problem. Journal of the operational research society pp. 167–
176 (1999)
[13] Gamma, E., Helm, R., Johnson, R., Vlissides, J.: Design Patterns: Elements of Reusable Object-Oriented Software. Pearson Education (1994)
19

[14] Grötschel, M.: Discrete mathematics in manufacturing. In: Malley, R.E.O.
(ed.) ICIAM 1991: Proceedings of the Second International Conference on
Industrial and Applied Mathematics. pp. 119–145. SIAM (1991)
[15] Hahn, P., Roth, A., Saltzman, M., Guignard, M.: Memory-aware parallelized rlt3 for solving quadratic assignment problems. Optimization online (2013), http://www.optimization-online.org/DB_HTML/2013/12/
4144.html
[16] Hahn, P.M., Zhu, Y.R., Guignard, M., Smith, J.M.: Exact solution of
emerging quadratic assignment problems. International Transactions in Operational Research 17(5), 525–552 (2010)
[17] Howes, L., Munshi, A.: Aparapi - AMD. http://developer.amd.com/
tools-and-sdks/opencl-zone/aparapi/, online: last accessed: Jan 2015
[18] Howes, L., Munshi, A.: The OpenCL specification. https://www.khronos.
org/registry/cl/specs/opencl-2.0.pdf, online: last accessed: Jan
2015
[19] Koopmans, T.C., Beckmann, M.J.: Assignment problems and the location
of economic activities. Econometrica 25, 53–76 (1957)
[20] Krüger, F., Maitre, O., Jiménez, S., Baumes, L., Collet, P.: Generic local
search (memetic) algorithm on a single GPGPU chip. In: Tsutsui, S., Collet, P. (eds.) Massively Parallel Evolutionary Computation on GPGPUs,
pp. 63–81. Natural Computing Series, Springer Berlin Heidelberg (2013)
[21] Liu, H., Abraham, A., Zhang, J.: A particle swarm approach to quadratic
assignment problems. In: Saad, A., Dahal, K., Sarfraz, M., Roy, R. (eds.)
Soft Computing in Industrial Applications, Advances in Soft Computing,
vol. 39, pp. 213–222. Springer Berlin Heidelberg (2007)
[22] Maitre, O.: Genetic programming on GPGPU cards using EASEA. In:
Tsutsui, S., Collet, P. (eds.) Massively Parallel Evolutionary Computation
on GPGPUs, pp. 227–248. Natural Computing Series, Springer Berlin Heidelberg (2013)
[23] Mason, A., Rönnqvist, M.: Solution methods for the balancing of jet turbines. Computers & OR 24(2), 153–167 (1997)
[24] Misevicius, A.: An implementation of the iterated tabu search algorithm
for the quadratic assignment problem. OR Spectrum 34(3), 665–690 (2012)
[25] Nyberg, A., Westerlund, T.: A new exact discrete linear reformulation
of the quadratic assignment problem. European Journal of Operational
Research 220(2), 314–319 (2012)
[26] Onwubolu, G.C., Sharma, A.: Particle swarm optimization for the assignment of facilities to locations. In: New Optimization Techniques in Engineering, pp. 567–584. Springer (2004)
20

[27] Owens, J.D., Luebke, D., Govindaraju, N., Harris, M., Krüger, J., Lefohn,
A.E., Purcell, T.J.: A survey of general-purpose computation on graphics
hardware. In: Computer graphics forum. vol. 26, pp. 80–113. Wiley Online
Library (2007)
[28] Peter Hahn and Miguel Anjos: QAPLIB home page. http://anjos.mgi.
polymtl.ca/qaplib/, online: last accessed: Jan 2015
[29] Phillips, A.T., Rosen, J.B.: A quadratic assignment formulation of the
molecular conformation problem. JOURNAL OF GLOBAL OPTIMIZATION 4, 229–241 (1994)
[30] Sahni, S., Gonzalez, T.: P-complete approximation problems. J. ACM
23(3), 555–565 (1976)
[31] Stone, J.E., Gohara, D., Shi, G.: Opencl: A parallel programming standard
for heterogeneous computing systems. Computing in science & engineering
12(3), 66 (2010)
[32] Stützle, T., Dorigo, M.: Aco algorithms for the quadratic assignment problem. New ideas in optimization pp. 33–50 (1999)
[33] Szwed, P., Chmiel, W., Kadłuczka, P.:
OpenCL implementation of PSO algorithm for the Quadratic Assignment Problem.
In: Rutkowski, L., Korytkowski, M., Scherer, R., Tadeusiewicz,
R., Zadeh, L.A., Zurada, J.M. (eds.) Artificial Intelligence and
Soft Computing, Lecture Notes in Computer Science, vol. Accepted for ICAISC’2015 Conference. Springer International Publishing (2015), http://home.agh.edu.pl/~pszwed/en/lib/exe/fetch.php?
media=papers:draft-icaics-2015-pso-qap-opencl.pdf
[34] Taillard, E.D.: Comparison of iterative searches for the quadratic assignment problem. Location Science 3(2), 87 – 105 (1995)
[35] Tsutsui, S., Fujimoto, N.: ACO with tabu search on GPUs for fast solution
of the QAP. In: Tsutsui, S., Collet, P. (eds.) Massively Parallel Evolutionary Computation on GPGPUs, pp. 179–202. Natural Computing Series,
Springer Berlin Heidelberg (2013)
[36] Zhou, Y., Tan, Y.: GPU-based parallel particle swarm optimization. In:
Evolutionary Computation, 2009. CEC’09. IEEE Congress on. pp. 1493–
1500. IEEE (2009)

21

