ADMM for the SDP relaxation of the QAP ∗
Danilo Elias Oliveira

†

Henry Wolkowicz

‡

Yangyang Xu

§

arXiv:1512.05448v1 [math.OC] 17 Dec 2015

December 18, 2015

Abstract
The semidefinite programming SDP relaxation has proven to be extremely strong for many hard
discrete optimization problems. This is in particular true for the quadratic assignment problem QAP,
arguably one of the hardest NP-hard discrete optimization problems. There are several difficulties that
arise in efficiently solving the SDP relaxation, e.g., increased dimension; inefficiency of the current primaldual interior point solvers in terms of both time and accuracy; and difficulty and high expense in adding
cutting plane constraints.
We propose using the alternating direction method of multipliers ADMM to solve the SDP relaxation.
This first order approach allows for inexpensive iterations, a method of cheaply obtaining low rank solutions, as well a trivial way of adding cutting plane inequalities. When compared to current approaches
and current best available bounds we obtain remarkable robustness, efficiency and improved bounds.

Keywords: Quadratic assignment problem, semidefinite programming relaxation, alternating direction
method of moments, large scale.
Classification code: 90C22, 90B80, 90C46, 90-08

Contents
1 Introduction

2

2 A New Derivation for the SDP Relaxation

2

3 A New ADMM Algorithm for the
3.1 Lower bound . . . . . . . . . . .
3.2 Feasible solution of QAP . . . .
3.3 Low-rank solution . . . . . . . .
3.4 Different choices for V, Vb . . . . .

SDP Relaxation
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .

4 Numerical experiments

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

5
6
7
8
8
9

5 Concluding Remarks

10

Index

10

∗ This

work is partially supported by NSERC and AFOSR.
of Combinatorics and Optimization, University of Waterloo.
‡ Dept. of Combinatorics and Optimization, University of Waterloo. Research supported by The Natural Sciences and
Engineering Research Council of Canada and by AFOSR. Email: hwolkowicz@uwaterloo.ca
§ Institute for Mathematics and its Applications (IMA), University of Minnesota
† Dept.

1

List of Tables
1

1

QAP Instances I and II. Requested tolerance 1e − 5. . . . . . . . . . . . . . . . . . . . . . . .

10

Introduction

The quadratic assignment problem (QAP), in the trace formulation is
p∗X := min hAXB − 2C, Xi,
X∈Πn

(1.1)

where A, B ∈ Sn are real symmetric n × n matrices, C is a real n × n matrix, h· , ·i denotes the trace inner
product, hY, Xi = trace Y X > , and Πn denotes the set of n × n permutation matrices. A typical objective of
the QAP is to assign n facilities to n locations while minimizing total cost. The assignment cost is the sum
of costs using the flows in Aij between a pair of facilities i, j multiplied by the distance in Bst between their
assigned locations s, t and adding on the location costs of a facility i in a position s given in Cis .
It is well known that the QAP is an NP-hard problem and that problems with size as moderate as
n = 30 still remain difficult to solve. Solution techniques rely on calculating efficient lower bounds. An
important tool for finding lower bounds is the work in [13] that provides a semidefinite programmming
(SDP), relaxation of (1.1). The methods of choice for SDP are based on a primal-dual interior-point, p-d ip, approach. These methods cannot solve large problems, have difficulty in obtaining high accuracy solutions
and cannot properly exploit sparsity. Moreover, it is very expensive to add on nonnegativity and cutting
plane constraints. The current state for finding bounds and solving QAP is given in e.g., [1, 2, 4, 7, 9].
In this paper we study an alternating direction method of multipliers (ADMM), for solving the SDP relaxation
of the QAP. We compare this with the best known results given in [9] and with the best known bounds
found at SDPLIB [5]. and with a p-d i-p methods based on the so-called HKM direction. We see that the
ADMM method is significantly faster and obtains high accuracy solutions. In addition there are advantages
in obtaining low rank SDP solutions that provide better feasible approximations for the QAP for upper
bounds. Finally, it is trivial to add nonnegativity and rounding constraints while iterating so as to obtain
significantly stronger bounds and also maintain sparsity during the iterations.
We note that previous success for ADMM for SDP in presented in [12]. A detailed survey article for
ADMM can be found in [3].

2

A New Derivation for the SDP Relaxation

We start the derivation from the following equivalent quadratically constrained quadratic problem
min hAXB − 2C, Xi
X

s.t. Xij Xik = 0, Xji Xki = 0, ∀i, ∀j 6= k,
2
Xij
− Xij = 0, ∀i, j,
n
n
X
X
2
2
Xij
− 1 = 0, ∀j,
Xij
− 1 = 0, ∀i.
i=1

(2.1)

j=1

Remark 2.1. Note that the quadratic orthogonality constraints X > X = I, XX > = I, and the linear row
and column sum constraints Xe = e, X > e = e can all be linearly represented using linear combinations of
those in (2.1).
In addition, the first set of constraints, the elementwise orthogonality of the row and columns of X, are
referred to as the gangster constraints. They are particularly strong constraints and enable many of the other
constraints to be redundant. In fact, after the facial reduction done below, many of these constraints also
become redundant. (See the definition of the index set J below.)
2

The Lagrangian for (2.1) is
L0 (X, U, V, W, u, v) =hAXB − 2C, Xi +

n X
X

(i)

Ujk Xij Xik +

i=1 j6=k

+

n
X

uj

j=1

n
X

!
2
Xij
−1

n X
X

(i)

Vjk Xji Xki +

X

i=1 j6=k

2
Wij (Xij
− Xij )

i,j



n
n
X
X
2
+
vi 
Xij
− 1 .

i=1

i=1

j=1

The dual problem is a maximization of the dual functional d0 ,
max d0 (U, V, W, u, v) := min L0 (X, U, V, W, u, v).

(2.2)

X

To simplify the dual problem, we homogenize the X terms in L0 by multiplying a unit scalar x0 to degree-1
terms and adding the single constraint x20 = 1 to the Lagrangian. We let
L1 (X, x0 , U, V, W, w0 , u, v) =hAXB − 2x0 C, Xi +

n X
X

(i)

Ujk Xij Xik +

i=1 j6=k

+

n
X
j=1

uj

n
X

!
2
Xij
−1

n X
X

(i)

Vjk Xji Xki +

i=1 j6=k

X

2
Wij (Xij
− x0 Xij )

i,j



n
n
X
X
2
+
vi 
Xij
− 1 + w0 (x20 − 1).

i=1

i=1

j=1

This homogenization technique is the same as that in [13]. The new dual problem is
max d1 (U, V, W, w0 , u, v) := min L1 (X, x0 , U, V, W, w0 , u, v).
X,x0

(2.3)

Note that d1 ≤ d0 . Hence, our relaxation still yields a lower bound to (2.1). In fact, the relaxations give
the same lower bound. This follows from strong duality of the trust region subproblem as shown in [13].
Let x = vec(X), y = [x0 ; x], and w = vec(W ), where x, w is the vectorization, columnwise, of X and W ,
respectively. Then
L1 (X, x0 , U, V, W, w0 , u, v) = y > [LQ + B1 (U ) + B2 (V ) + Arrow(w, w0 ) + K1 (u) + K2 (v)] y − e> (u + v) − w0 ,
where
K1 (u) = blkdiag(0, u ⊗ I), K2 (v) = blkdiag(0, I ⊗ v),


w0
− 21 w>
Arrow(w, w0 ) =
− 12 w Diag(w)
and
B1 (U ) = blkdiag(0, Ũ ),

B2 (V ) = blkdiag(0, Ṽ ).

Here, Ũ and Ṽ are n × n block matrices. Ũ has zero diagonal blocks and the (j, k)-th off-diagonal block to
(1)
(n)
be the diagonal matrix Diag(Ujk , . . . , Ujk ) for all j 6= k, and Ṽ has zero off-diagonal blocks and the i-th


(i)
(i)
0
V12 · · · V1n
 (i)
(i) 
0
· · · V2n 
 V21
diagonal block to be 
..
.. 
..
. Hence, the dual problem (2.3) is
 ..
.
 .
.
. 
(i)
(i)
Vn1 Vn2 · · ·
0
max − e> (u + v) − w0

(2.4)

s.t. LQ + B1 (U ) + B2 (V ) + Arrow(w, w0 ) + K1 (u) + K2 (v)  0.

3

Taking the dual of (2.4), we have the SDP relaxation of (2.1):
min hLQ , Y i
s.t. GJ (Y ) = E00 , diag(Ȳ ) = y0 ,
n
X
trace(Ỹii ) = 1, ∀i,
Ỹii = I,

(2.5)

i=1

Y  0,
where Ỹij is an n × n matrix for each (i, j), and we have assumed the block structure

Y =

y00
y0

y0>
Ȳ


;

Ȳ made of n × n block matrices Ỹ = (Ỹij ).

(2.6)

The index set J and the gangster operator GJ are defined properly below in Definition 2.1. (By abuse of
notation this is done after the facial reduction which results in a smaller J.)
Remark 2.2. If one more feasible quadratic constraint q(X) can be added to (2.1) and q(X) cannot be
linearly represented by those in (2.1), the relaxation following the same derivation as above can be tighter.
We conjecture that no more such q(X) exists, and thus (2.5) is the tightest among all Lagrange dual relaxation
from a quadratically constrained program like (2.1). However, this does not mean that more linear inequality
constraints cannot be added, i.e., linear cuts.
Theorem 2.1 ( [13]). The matrix Y is feasible for (2.5) if, and only if, it is feasible for (3.1).
2

As above, let x = vec X ∈ Rn be the vectorization of X by column. Y is the original matrix variable of


>
1
1
the SDP relaxation before the facial reduction. It can be motivated from the lifting Y =
.
vec X
vec X
The SDP relaxation of QAP presented in [13] uses facial reduction to guarantee strict feasibility. The
SDP obtained is
p∗R := minR hLQ , V̂ RV̂ > i
(2.7)
s.t. GJ (V̂ RV̂ > ) = E00
R  0,
where the so-called gangster operator, GJ , fixes all elements indexed by J and zeroes out all others,




0
−vec(C)>
1
0
LQ =
,
V̂ = 1
−vec(C)
B⊗A
ne V ⊗ V

(2.8)

with e being the vector of all ones, of appropriate
dimension and V ∈ Rn×(n−1) being a basis matrix of the
"
#
In−1
2
orthogonal complement of e, e.g., V =
. We let Y = V̂ RV̂ > ∈ Sn +1 .
−e
Lemma 2.1 ( [13]). The matrix R̂ defined by
"
#
1
0
(n−1)2 +1
R̂ :=
∈ S++
1
0 n2 (n−1) (nIn−1 − En−1 ) ⊗ (nIn−1 − En−1 )
is (strictly) feasible for (2.7).
2

2

Definition 2.1. The gangster operator GJ : Sn +1 → Sn +1 and is defined by

Yij if (i, j) ∈ J or (j, i) ∈ J
GJ (Y )ij =
0
otherwise
4

By abuse of notation, we let the same symbol denote the projection onto R|J| . We get the two equivalent
primal constraints:
GJ (V̂ RV̂ > ) = E00 ∈ Sn

2

+1

GJ (V̂ RV̂ > ) = GJ (E00 ) ∈ R|J| .

;
2

Therefore, the dual variable for the first form is Y ∈ Sn +1 . However, the dual variable for the second form
2
is y ∈ R|J| with the adjoint now yielding Y = GJ∗ (y) ∈ Sn +1 obtained by symmetrization and filling in the
missing elements with zeros.
The gangster index set, J is defined to be (00) union the set of of indices i < j in the matrix Ȳ in (2.6)
corresponding to:
1. the off-diagonal elements in the n diagonal blocks;
2. the diagonal elements in the off-diagonal blocks except for the last column of off-diagonal blocks and
also not the (n−2), (n−1) off-diagonal block. (These latter off-diagonal block constraints are redundant
after the facial reduction.)
We note that the gangster operator is self-adjoint, GJ∗ = GJ . Therefore, the dual of (2.7) can be written
as the following.
d∗Y := max hE00 , Y i
(= Y00 )
Y
(2.9)
>
>
s.t. V̂ GJ (Y )V̂  V̂ LQ V̂
Again by abuse of notation, using the same symbol twice, we get the two equivalent dual constraints:
V̂ > GJ (Y )V̂  V̂ > LQ V̂ ;

V̂ > GJ∗ (y)V̂  V̂ > LQ V̂ .
2

As above, the dual variable for the first form is Y ∈ Sn +1 and for the second form is y ∈ R|J| . We have
used G ∗ for the second form to emphasize that only the first form is self-adjoint.
Lemma 2.2 ( [13]). The matrices Ŷ ,Ẑ, with M > 0 sufficiently large, defined by
"
#
n
0
(n−1)2 +1
(n−1)2 +1
Ŷ := M
.
, Ẑ := V̂ > LQ V̂ − V̂ > GJ (Ŷ )V̂ ∈ S++
∈ S++
0 In ⊗ (In − En )
and are (strictly) feasible for (2.9).

3

A New ADMM Algorithm for the SDP Relaxation

We can write (2.7) equivalently as
min hLQ , Y i, s.t. GJ (Y ) = E00 , Y = V̂ RV̂ > , R  0.
R,Y

(3.1)

The augmented Lagrange of (3.1) is
LA (R, Y, Z) = hLQ , Y i + hZ, Y − V̂ RV̂ > i +

β
kY − V̂ RV̂ > k2F .
2

(3.2)

Recall that (R, Y, Z) are the primal reduced, primal, and dual variables respectively. We denote (R, Y, Z)
n
as the current iterate. We let Srn
+ denote the matrices in S+ with rank at most r. Our new algorithm is an
application of the alternating direction method of multipliers ADMM, that uses the augmented Lagrangian
in (3.2) and performs the following updates for (R+ , Y+ , Z+ ):
R+ = arg min LA (R, Y, Z),

(3.3a)

Y+ = arg min LA (R+ , Y, Z),

(3.3b)

R∈Srn
+
Y ∈Pi

Z+ = Z + γ · β(Y+ − V̂ R+ V̂ > ),
5

(3.3c)

where the simplest case for the polyhedral constraints Pi is the linear manifold from the gangster constraints:
P1 = {Y ∈ Sn

2

+1

: GJ (Y ) = E00 }

We use this notation as we add additional simple polyhedral constraints. The second case is the polytope:
P2 = P1 ∩ {0 ≤ Y ≤ 1}.
Let V̂ be normalized such that V̂ > V̂ = I. Then if r = n, the R-subproblem can be explicitly solved by
R+

=

arg minR0 hZ, Y − V̂ RV̂ > i + β2 kY − V̂ RV̂ > k2F
2

arg minR0 Y − V̂ RV̂ > + β1 Z
F

= arg minR0 R − V̂ > Y + β1 Z V̂

 
= PS+ V̂ > Y + β1 Z V̂ ,
=

2

(3.4)

F

where S+ denotes the SDP cone, and PS+ is the projection to S+ . For any symmetric matrix W , we have
>
PS+ (W ) = U+ Σ+ U+
,

where (U+ , Σ+ ) contains the positive eigenpairs of W and (U− , Σ− ) the negative eigenpairs.
If i = 1 in (3.3b), the Y -subproblem also has closed-form solution:
Y+ = arg min hLQ , Y i + hZ, Y − V̂ R+ V̂ > i +
GJ (Y )=E00

LQ + Z
β
GJ (Y )=E00


LQ + Z
>
=E00 + GJ c V̂ R+ V̂ −
β

= arg min

β
kY − V̂ R+ V̂ > k2F
2

2

Y − V̂ R+ V̂ > +

F

(3.5)

The advantage of using ADMM is that its complexity only slightly increases while we add more constraints to (2.7) to tighten the SDP relaxation. If 0 ≤ V̂ RV̂ > ≤ 1 is added in (2.7), then we have constraint
0 ≤ Y ≤ 1 in (3.1) and reach to the problem
p∗RY := min hLQ , Y i, s.t. GJ (Y ) = E00 , 0 ≤ Y ≤ 1, Y = V̂ RV̂ > , R  0.
R,Y

(3.6)

The ADMM for solving (3.6) has the same R-update and Z-update as those in (3.3), and the Y -update is
changed to



LQ + Z 
>
c
Y+ = E00 + min 1, max 0, GJ V̂ R+ V̂ −
.
(3.7)
β
With nonnegativity constraint, the less-than-one constraint is redundant but makes the algorithm converge
faster.

3.1

Lower bound

If we solve (2.7) or (3.1) exactly or to a very high accuracy, we get a lower bound of the original QAP.
However, the problem size of (2.7) or (3.1) can be extremely large, and thus having an exact or highly
accurate solution may take extremely long time. In the following, we provide an inexpensive way to get a
lower bound from the output of our algorithm that solves (3.1) to a moderate accuracy. Let (Rout , Y out , Z out )
be the output of the ADMM for (3.6).

6

Lemma 3.1. Let
R := {R  0},

Z := {Z : V̂ > Z V̂  0}.

Y := {Y : GJ (Y ) = E00 , 0 ≤ Y ≤ 1},

Define the ADMM dual function
g(Z) := min {hLQ + Z, Y i}.
Y ∈Y

Then the dual problem of ADMM (3.6) is defined as follows and satisfies weak duality.
d∗Z

:=

max g(Z)

≤

p∗R .

Z∈Z

Proof. The dual problem of (3.6) can be derived as
d∗Z := max
Z

min

R∈R,Y ∈Y

hLQ , Y i + hZ, Y − V̂ RV̂ > i

= max min hLQ , Y i + hZ, Y i + min hZ, −V̂ RV̂ > i
Y ∈Y

Z

R∈R

= max min hLQ , Y i + hZ, Y i + min hV̂ > Z V̂ , −Ri
Y ∈Y

Z

R∈R

= max min hLQ + Z, Y i,
Z∈Z Y ∈Y

= max g(Z)
Z∈Z

Weak duality follows in the usual way by exchanging the max and min.
For any Z ∈ Z, we have g(Z) is a lower bound of (3.6) and thus of the original QAP. We use the dual

function value of the projection g PZ (Z out ) as the lower bound, and next we show how to get PZ (Z̃) for
any symmetric matrix Z̃.
Let V̂⊥ be the orthonormal basis of the null space of V̂ . Then V̄ = (V̂ , V̂⊥ ) is an orthogonal matrix. Let


W11 W12
V̄ > Z V̄ = W =
, and we have
W21 W22
V̂ > Z V̂  0 ⇔ V̂ > Z V̂ = V̂ > V̄ W V̄ > V̂ = W11  0.
Hence,
PZ (Z̃) = arg min kZ − Z̃k2F
Z∈Z

= arg min kV̄ W V̄ > − Z̃k2F
W11 0

= arg min kW − V̄ > Z̃ V̄ k2F
W11 0


=

PS− (W̃11 ) W̃12
W̃21
W̃22


,

where S− denotes the negative semidefinite cone, and we have assumed V̄ > Z̃ V̄ =



W̃11
W̃21

W̃12
W̃22


. Note that

PS− (W11 ) = −PS+ (−W11 ).

3.2

Feasible solution of QAP

Let (Rout , Y out , Z out ) be the output of the ADMM for (3.6). Assume the largest eigenvalue and the corresponding eigenvector of Y are λ and v. We let X out be the matrix reshaped from the second through the
last elements of the first column of λvv > . Then we solve the linear program
maxhX out , Xi, s.t. Xe = e, X > e = e, X ≥ 0
X

by simplex method that gives a basic optimal solution, i.e., a permutation matrix.
7

(3.8)

3.3

Low-rank solution

Instead of finding a feasible solution through (3.8), we can directly get one by restricting R to a rank-one
matrix, i.e., rank(R) = 1 and R ∈ S+ . With this constraint, the R-update can be modified to


Z
>
V̂ ,
(3.9)
R+ = PS+ ∩R1 V̂ Y +
β
where R1 = {R : rank(R) = 1} denotes the set of rank-one matrices. For a symmetric matrix W with largest
eigenvalue λ > 0 and corresponding eigenvector w, we have
PS+ ∩R1 = λww> .

3.4

Different choices for V, Vb

The matrix Vb is essential in the steps of the algorithm, see e.g., (3.4). A sparse Vb helps in the projection if
one is using a sparse eigenvalue code. We have compared several. One is based on applying a QR algorithm
to the original simple V from the definition of V̂ in (2.8). The other two are based on the approach in [10]
and we present the most successful here. The orthogonal V we use is



 
1








1
 I n ⊗ √1
I n ⊗ 1 
h i
 1 


 b 2 c



2




c
b
2 −1
4
−1  . . . Vb 
V = 






−1 
0(n−2b n c),b n c



2
2
0(n−4b n c),b n c
4

n×n−1

4

i.e., the block matrix consisting of t blocks formed from Kronecker products along with one block Vb to
complete the appropriate size so that V > V = In−1 , V > e = 0. We take advantage of the 0, 1 structure of the
Kronecker blocks and delay the scaling for the normalization till the end. The main work in the low rank
projection part of the algorithm is to evaluate one (or a few) eigenvalues of W = Vb > (Y + β1 Z)V̂ to obtain
the update R+ .


1
ρ w>
Y + Z=
.
w W̄
β
We let
K := V ⊗ V,

√
α = 1/ 2,

1
v = √ e,
2n

x=

 
x1
.
x̄

The structure for Vb in (2.8) means that we can evaluate the product for W x as

α
v

0
K

> 

ρ
w

w>
W̄


α
v


0
x
K

=
=
=
=


> 


α 0
ρ w>
αx1
v K
w W̄
x1 v + K x̄



>
>
ραx1 + w (x1 v + K x̄)
α v
0 K>
αx1 w + W̄ (x1 v + K x̄)

 2
ρα x1 + αw> (x1 v + K x̄) + v > αx1 w + W̄ (x1 v + K x̄)

K > αx1 w + W̄ (x1 v + K x̄)

 2

ρα x1 + αw> + v > W̄ (x1 v + K x̄) + v > (αx1 w)

.
K > αx1 w + W̄ (x1 v + K x̄)

We emphasize that V ⊗ V = (V̄ ⊗ V̄ )/(D ⊗ D), where V̄ denotes the unscaled V , D is the diagonal
matrix of scale factors to obtain the orthogonality in V , and / denotes the MATLAB division on the right,
multiplication by the inverse on the right. Therefore, we can evaluate


K > W̄ K = (V ⊗ V )> W̄ (V ⊗ V ) = (V̄ ⊗ V̄ )> (D ⊗ D)\W̄ /(D ⊗ D) (V̄ ⊗ V̄ ).
8

4

Numerical experiments

We illustrate our results in Table 1 on the forty five QAP instances I and II, see [5, 6, 9]. The optimal
solutions are in column 1 and current best known lower bounds from [9] are in column 3 marked bundle.
The p-d i-p lower bound is given in the column marked HKM-FR. (The code failed to find a lower bound on
several problems marked −1111.) These bounds were obtained using the facially reduced SDP relaxation
and exploiting the low rank (one and two) of the constraints. We used SDPT3 [11].1
Our ADMM lower bound follows in column 4. We see that it is at least as good as the current best
known bounds in every instance. The percent improvement is given in column 7. We then present the best
upper bounds from our heuristics in column 5. This allows us to calculate the percentage gap in column 6.
The CPU seconds are then given in the last columns 8 − 9 for the high and low rank approaches, respectively.
The last two columns are the ratios of CPU times. Column 10 is the ratio of CPU times for the 5 decimal and
12 decimal tolerance for the high rank approach. All the ratios for the low rank approach are approximately
1 and not included. The quality of the bounds did not change for these two tolerances. However, we consider
it of interest to show that the higher tolerance can be obtained.
The last column 11 is the ratio of CPU times for the 12 decimal tolerance of the high rank approach in
column 8 with the CPU times for 9 decimal tolerance for the HKM approach. We emphasize that the lower
bounds for the HKM approach were significantly weaker.
We used MATLAB version 8.6.0.267246 (R2015b) on a PC Dell Optiplex 9020 64-bit, with 16 Gig,
running Windows 7.
We heuristically set γ = 1.618 and β = n3 in ADMM. We used two different tolerances 1e − 12, 1e − 5.
Solving the SDP to the higher accuracy did not improve the bounds. However, it is interesting that the
ADMM approach was able to solve the SDP relaxations to such high accuracy, something the p-d i-p
approach has great difficulty with. We provide the CPU times for both accuracies. Our times are significantly
lower than those reported in [4, 9], e.g., from 10 hours to less than an hour.
We emphasize that we have improved bounds for all the SDP instances and have provably found exact
solutions six of the instances Had12,14,16,18, Rou12, Tai12a. This is due to the ability to add all the
nonnegativity constraints and rounding numbers to 0, 1 with essentially zero extra computational cost. In
addition, the rounding appears to improve the upper bounds as well. This was the case for both using
tolerance of 12 or only 5 decimals in the ADMM algorithm.
1 We

do not include the times as they were much greater than for the ADMM approach, e.g., hours instead of minutes and
a day instead of an hour.

9

Esc16a
Esc16b
Esc16c
Esc16d
Esc16e
Esc16g
Esc16h
Esc16i
Esc16j
Had12
Had14
Had16
Had18
Had20
Kra30a
Kra30b
Kra32
Nug12
Nug14
Nug15
Nug16a
Nug16b
Nug17
Nug18
Nug20
Nug21
Nug22
Nug24
Nug25
Nug27
Nug28
Nug30
Rou12
Rou15
Rou20
Scr12
Scr15
Scr20
Tai12a
Tai15a
Tai17a
Tai20a
Tai25a
Tai30a
Tho30

1.
opt
value

2.
Bundle [9]
LowBnd

3.
HKM-FR
LowBnd

4.
ADMM
LowBnd

5.
feas
UpBnd

6.
ADMM
%gap

7. ADMM
vs Bundle
%Impr LowBnd

8 Tol5
cpusec
HighRk

9 Tol5
cpusec
LowRk

10 Tol12/5
cpuratio
HighRk

11 HKM
cpuratio
Tol 9

68
292
160
16
28
26
996
14
8
1652
2724
3720
5358
6922
149936
91420
88700
578
1014
1150
1610
1240
1732
1930
2570
2438
3596
3488
3744
5234
5166
6124
235528
354210
725522
31410
51140
110030
224416
388214
491812
703482
1167256
1818146
88900

59
288
142
8
23
20
970
9
7
1643
2715
3699
5317
6885
136059
81156
79659
557
992
1122
1570
1188
1669
1852
2451
2323
3440
3310
3535
4965
4901
5803
223680
333287
663833
29321
48836
94998
222784
364761
451317
637300
1041337
1652186
77647

50
276
132
-12
13
11
909
-21
-4
1641
2709
3678
5287
6848
-1111
-1111
-1111
530
960
1071
1528
1139
1622
1802
2386
2386
3396
-1111
-1111
-1111
-1111
-1111
221161
323235
642856
23973
42204
83302
215637
349586
441294
619092
1096657
-1111
-1111

64
290
154
13
27
25
977
12
8
1652
2724
3720
5358
6922
143576
87858
85775
568
1011
1141
1600
1219
1708
1894
2507
2382
3529
3402
3626
5130
5026
5950
235528
350217
695181
31410
51140
106803
224416
377101
476525
671675
1096657
1706871
86838

72
300
188
18
32
28
996
14
14
1652
2724
3720
5358
6930
169708
105740
103790
632
1022
1306
1610
1356
1756
2160
2784
2706
3940
3794
4060
5822
5730
6676
235528
367782
765390
38806
58304
138474
224416
412760
546366
750450
1271696
1942086
102760

11.76
3.42
21.25
31.25
17.86
11.54
1.91
14.29
75.00
0.00
0.00
0.00
0.00
0.12
17.43
19.56
20.31
11.07
1.08
14.35
0.62
11.05
2.77
13.78
10.78
13.29
11.43
11.24
11.59
13.22
13.63
11.85
0.00
4.96
9.68
23.55
14.01
28.78
0.00
9.19
14.20
11.20
15.00
12.94
17.91

7.35
0.68
7.50
31.25
14.29
19.23
0.70
21.43
12.50
0.54
0.33
0.56
0.77
0.53
5.01
7.33
6.90
1.90
1.87
1.65
1.86
2.50
2.25
2.18
2.18
2.42
2.47
2.64
2.43
3.15
2.42
2.40
5.03
4.78
4.32
6.65
4.51
10.73
0.73
3.18
5.13
4.89
4.74
3.01
10.34

2.30e+01
3.87e+00
1.09e+01
2.14e+01
3.02e+01
4.24e+01
4.91e+00
1.37e+02
8.95e+01
1.02e+01
3.23e+01
1.75e+02
4.49e+02
3.85e+02
5.88e+03
4.36e+03
3.57e+03
2.60e+01
7.15e+01
9.10e+01
1.81e+02
9.35e+01
2.31e+02
4.16e+02
4.76e+02
1.41e+03
2.07e+03
1.20e+03
3.12e+03
5.11e+03
4.11e+03
7.36e+03
2.76e+01
3.12e+01
1.67e+02
4.40e+00
1.38e+01
1.53e+03
1.79e+00
2.74e+01
6.50e+01
1.28e+02
3.09e+02
1.25e+03
2.83e+03

4.02
4.55
8.09
3.69
4.29
4.27
3.53
4.30
4.80
1.08
1.69
3.15
6.00
12.15
149.32
170.57
200.26
1.04
1.87
3.31
3.06
3.19
4.34
5.47
11.56
15.32
21.82
29.64
39.23
78.18
83.38
133.38
0.93
2.70
10.31
1.17
2.41
9.61
0.90
2.35
4.52
10.10
38.48
142.55
164.86

4.14
2.15
4.53
4.87
4.80
2.72
2.33
2.39
3.83
1.06
1.19
1.04
2.22
4.20
2.22
3.01
4.28
6.61
5.06
5.90
3.28
6.23
3.63
2.43
3.75
1.68
1.39
3.29
1.65
1.58
2.17
1.76
0.98
8.68
10.90
2.40
1.84
1.15
1.04
14.69
7.31
14.32
5.58
10.51
4.74

9.37
8.08
4.88
10.22
8.79
8.63
10.60
8.76
7.93
5.91
10.46
12.51
13.28
14.53
1111.11
1111.11
1111.11
5.93
8.43
7.79
12.24
11.83
13.13
15.23
14.35
14.95
13.90
1111.11
1111.11
1111.11
1111.11
1111.11
6.90
9.46
16.08
5.79
10.75
17.96
6.70
10.34
12.04
15.85
1111.11
1111.11
1111.11

Table 1: QAP Instances I and II. Requested tolerance 1e − 5.

5

Concluding Remarks

In this paper we have shown the efficiency of using the ADMM approach in solving the SDP relaxation
of the QAP problem. In particular, we have shown that we can obtain high accuracy solutions of the
SDP relaxation in less significantly less cost than current approaches. In addition, the SDP relaxation
includes the nonnegativity constraints at essentially no extra cost. This results in both a fast solution and
improved lower and upper bounds for the QAP.
In a forthcoming study we propose to include this in a branch and bound framework and implement it in
a parallel programming approach, see e.g., [8]. In addition, we propose to test the possibility of using warm
starts in the branching/bounding process and test it on the larger test sets such as used in e.g., [7].

10

Index
J, gangster index set, 5
GJ , gangster operator, 4
vec, 3
R̂, 4
Ŷ , 5
Ẑ, 5
Srn
+ , 5
∗
dZ , 7
d∗Y , 5
e, ones vector, 4
g(Z), 7
p∗R , 4
p∗X , 2
p∗RY , 6
2
P1 = {Y ∈ Sn +1 : GJ (Y ) = E00 }, 6
P2 = P1 ∩ {0 ≤ Y ≤ 1}, 6
R := {R  0}, 7
Y := {Y : GJ (Y ) = E00 , 0 ≤ Y ≤ 1}, 7
Z := {Z : V̂ > Z V̂  0}, 7
QAP, quadratic assignment problem, 2
SDP, semidefinite programmming, 2
alternating direction method of multipliers, 2, 5
augmented Lagrange, 5
dual function, 7
facial reduction, 4
gangster constraints, 2, 6
gangster index set, J, 5
gangster operator, GJ , 4
lifting, 4
linear cuts, 4
ones, 4
optimal
optimal
optimal
optimal

value
value
value
value

ADMM relaxation, p∗RY , 6
QAP, p∗X , 2
dual SDP relaxation, d∗Y , 5
primal SDP relaxation, p∗R , 4

primal-dual interior-point, p-d i-p, 2
quadratic assignment problem, 2
semidefinite programmming, 2
strictly feasible pair, (R̂, Ŷ , Ẑ), 4, 5
trace inner product, hY, Xi = trace Y X > , 2
11

References
[1] K.M. Anstreicher. Recent advances in the solution of quadratic assignment problems. Math. Program.,
97(1-2, Ser. B):27–42, 2003. ISMP, 2003 (Copenhagen).
[2] K.M. Anstreicher and N.W. Brixius. A new bound for the quadratic assignment problem based on
convex quadratic programming. Math. Program., 89(3, Ser. A):341–357, 2001.
[3] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning
via the alternating direction method of multipliers. Found. Trends Machine Learning, 3(1):1–122, 2011.
[4] S. Burer and R.D.C. Monteiro. Local minima and convergence in low-rank semidefinite programming.
Math. Program., 103(3, Ser. A):427–444, 2005.
[5] R.E. Burkard, S. Karisch, and F. Rendl. QAPLIB – a quadratic assignment problem library. European
J. Oper. Res., 55:115–119, 1991. www.opt.math.tu-graz.ac.at/qaplib/.
[6] R.E. Burkard, S.E. Karisch, and F. Rendl. QAPLIB—a quadratic assignment problem library. J. Global
Optim., 10(4):391–403, 1997.
[7] E. de Klerk and R. Sotirov. Exploiting group symmetry in semidefinite programming relaxations of the
quadratic assignment problem. Math. Program., 122(2, Ser. A):225–246, 2010.
[8] R. Jain and P. Yao. A parallel approximation algorithm for positive semidefinite programming. In 2011
IEEE 52nd Annual Symposium on Foundations of Computer Science—FOCS 2011, pages 463–471.
IEEE Computer Soc., Los Alamitos, CA, 2011.
[9] F. Rendl and R. Sotirov. Bounds for the quadratic assignment problem using the bundle method. Math.
Program., 109(2-3, Ser. B):505–524, 2007.
[10] H. Sun, N. Wang, T.K. Pong, and H. Wolkowicz. Eigenvalue, quadratic programming, and semidefinite
programming bounds for a cut minimization problem. Comput. Optim. Appl., (accepted Aug. 2015):to
appear, 2015. 32 pages, submitted.
[11] K.C. Toh, M.J. Todd, and R.H. Tütüncü. SDPT3—a MATLAB software package for semidefinite
programming, version 1.3. Optim. Methods Softw., 11/12(1-4):545–581, 1999. Interior point methods.
[12] Z. Wen, D. Goldfarb, and W. Yin. Alternating direction augmented Lagrangian methods for semidefinite
programming. Math. Program. Comput., 2(3-4):203–230, 2010.
[13] Q. Zhao, S.E. Karisch, F. Rendl, and H. Wolkowicz. Semidefinite programming relaxations for the
quadratic assignment problem. J. Comb. Optim., 2(1):71–109, 1998. Semidefinite programming and
interior-point approaches for combinatorial optimization problems (Fields Institute, Toronto, ON, 1996).

12

