Department of Policy and Planning Sciences
Discussion Paper Series

No.1368

Centering ADMM for the Semidefinite Relaxation of the
QAPbyShin-ichi KANOH , Akiko YOSHISE

January 2020
Revised April 2020

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANCentering ADMM for the Semideﬁnite Relaxation of the QAP∗
Shin-ichi Kanoh†and Akiko Yoshise‡
January 2020
Revised April 2020

Abstract
We propose a new method for solving the semideﬁnite (SD) relaxation of the quadratic assignment
problem (QAP), called Centering ADMM. Centering ADMM is an alternating direction method ofmultipliers (ADMM) combining the centering steps used in the interior-point method. The ﬁrst stageof Centering ADMM updates the iterate so that it approaches the central path by incorporatinga barrier function term into the objective function, as in the interior-point method. If the current
iterate is suﬃciently close to the central path with a suﬃciently small value of the barrier parameter,
the method switches to the standard version of ADMM. We show that Centering ADMM (notemploying a dynamic update of the penalty parameter) has global convergence properties. To
observe the eﬀect of the centering steps, we conducted numerical experiments with SD relaxation
problems of instances in QAPLIB. The results demonstrate that the centering steps are quite eﬃcientfor some classes of instances.

Key words: Quadratic assignment problem; Semideﬁnite relaxation; Alternating direction method ofmultipliers (ADMM); Interior-point method; Centering step; Barrier function
AMS subject classiﬁcations: 90C05, 90C22, 90C25

1Introduction

The quadratic assignment problem (QAP) in the trace formulation [7] is given by
QAP minimize
subject to⟨F XD − C, X⟩,
X ∈ Πn ,(1)

where F, D ∈ Sn are n × n real symmetric matrices, C is a real n × n matrix, ⟨·, ·⟩ denotes the trace innerproduct ⟨A, B⟩ = Tr(AT B), and Πn denotes the set of all n × n permutation matrices. The QAP was
initially introduced to describe a location problem where the task is to assign n facilities to n locations
in a way that minimizes the total cost [10, 11]. It has many applications in areas as divergent as network
design, VSLI design, and image processing (see, e.g., [15, 16, 5]).
∗ This research was supported by the Japan Society for the Promotion of Science through a Grant-in-Aid for Challenging
Exploratory Research (17K18946) and a Grant-in-Aid for Scientiﬁc Research ((B)19H02373) of the Ministry of Education,
Culture, Sports, Science and Technology of Japan.
† Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan.
email: s1930138@s.tsukuba.ac.jp‡ Corresponding author. Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba, Ibaraki
305-8573, Japan. email: yoshise@sk.tsukuba.ac.jp1

It is known that the QAP is NP-hard (see, e.g., [20]) and it has remained diﬃcult to solve even if thesize of the problem is moderate, e.g., n = 30 [1]. This fact implies that ﬁnding better lower and upperbounds of the optimal value is quite important to solve it. An eﬃcient tool for ﬁnding such bounds
is semideﬁnite (SD) relaxation (see, e.g., [21, 18, 19, 6, 13, 17, 14]). SD relaxation in [21] uses facial
reduction to guarantee strict feasibility for both the relaxed problem and its dual and simpliﬁes theconstraints by making many of them redundant.
However, SD relaxation still often forces us to solve a large-scale semideﬁnite programs. Recently,
Oliveira, Workowicz, and Xu [14] performed computational experiments showing that their alternating
direction method of multipliers (ADMM) is promising for solving the SD relaxation of the QAP. The
authors derived the update formula for solving the SD relaxation of the QAP proposed in [21]. They also
derived upper and lower bounds of the QAP from the solution obtained by their ADMM and compared
their bounds with existing bounds.
The ADMM is a ﬁrst-order method, which requires less computation per iteration and is highly scalable.
However, its computation often becomes unstable and requires a signiﬁcant number of iterations tocompute an accurate solution. On the other hand, the interior-point method is a second-order method
and computationally expensive. Still, it is known to be stable and to require only a small number of
iterations due to the centering steps, which forces the current iterate closer to the central path. Our
motivation is to accelerate the ADMM proposed in [14] by combining the beneﬁts of the ADMM andthe interior-point method. In this paper, we devise a new algorithm, called Centering ADMM, which isan ADMM combining the centering steps of the interior-point method in its ﬁrst stage.
A similar approach was taken by Lin et al. [12], where the authors proposed an ADMM-Based Interior
Point Method (ABIP) for solving large-scale linear programs. Their ABIP and our Centering ADMM
are similar in the sense that both methods add a centering eﬀect in the search direction by using a barrier
function. However, Centering ADMM is diﬀerent from ABIP for the following reasons: Centering ADMM
is limited to solving the SD relaxation of the QAP. Also, it performs centering steps only in the primal
problem and in its ﬁrst stage. It reverts to (the standard) ADMM if the current iterate is suﬃciently
close to the central path with a suﬃciently small value of the barrier parameter. On the other hand,
ABIP employs centering steps in both primal and dual problems using the homogeneous self-dual formof linear programs throughout its iterations.
To observe the eﬀect of the centering steps, we conducted numerical experiments using instances in
QAPLIB [3] and compared the solutions obtained with the ADMM in [14] and with Centering ADMM.
The results demonstrate that the centering steps are quite eﬃcient for some classes of instances.
The organization of the paper is as follows: After giving a brief introduction of (the standard) ADMM in
section 2, we describe its in details for solving the SD relaxation of the QAP proposed in [14] in section
3. We provide our new method, Centering ADMM, in section 4. Then, we show that Centering ADMM
(not employing a dynamic update of the penalty parameter) has global convergence properties in section
5. In section 6, we numerically compare these two methods in terms of their lower bounds of the QAP
for instances in QAPLIB [3].

2Standard ADMM

Here, we give a brief introduction of the (standard) alternating direction method of multipliers. To
contrast with Centering ADMM, we will refer to the method as the Standard ADMM in the following.
The Standard ADMM was proposed by Glowinski and Marrocco [9] and Gabay and Mercier [8] forsolving the following type of optimization problem:2

minimizesubject tof (x) + g(z)Ax + Bz = c,(2)where x ∈ Rn , z ∈ Rm , A ∈ Rk×n , B ∈ Rk×m , c ∈ Rk and f : Rn → R ∪ {+∞} and g : Rm → R ∪ {+∞}are closed proper convex functions.
By introducing a penalty parameter ρ > 0, the augmented Lagrangian function for problem (2) is givenbyρLρ (x, z, y) := f (x) + g(z) + ⟨y, Ax + Bz − c⟩ + ∥Ax + Bz − c∥2 ,2where y ∈ Rk is the dual variable or Lagrange multiplier. Using the augmented Lagrangian function,
ADMM updates the variables (xk , y k , z k ) as follows:
xk+1:= argminx Lρ (x, z k , y k ),z k+1:= argminz Lρ (xk+1 , z, y k ),y:= y k + ρ(Axk+1 + Bz k+1 − c).

k+1

The third update formula is a feature of the Standard ADMM. It updates the dual variable y in its
gradient direction, i.e., in the direction that increases the objective function value of the dual problem.
For this reason, the Standard ADMM is sometimes considered to be a method that solves the dual
problem.

3Standard ADMM for the SD relaxation of the QAPAs shown in, e.g., [16], the set Πn of all permutation matrices can be represented as
Πn = On ∩ En ∩ Nn = On ∩ En ∩ Zn ,where On := {X ∈ Rn×n |XX T = X T X = I}, En := {X ∈ Rn×n |Xe = X T e = e}, Nn := {X ∈Rn×n |X ≥ 0}, Zn := {X ∈ Rn×n |X ◦ X − X = O}, I ∈ Rn×n is the identity matrix, e ∈ Rn is the vectorof ones, A ◦ B is the Hadamard product of A ∈ Rn×n and B ∈ Rn×n . Using this fact, Zhao et al. [21]showed that the QAP (1) is equivalent to the following problem QAPO :QAPOminimizeX
subject to⟨F XD − C, X⟩
XX T = I,
X T X = I,
∥Xe − e∥2 + ∥X T e − e∥2 = 0,
X ◦ X − X = O.

We also deﬁne S n := {X ∈ Rn×n | X = X T }, and for any A, B ∈ S n , we say A ⪰ B if A − B is
positive semideﬁnite. By considering the dual problem of the Lagrange dual of QAPO and projecting
the dual problem onto the minimal face, they also showed that the following problem QAPR1 gives an
SDP relaxation problem for the QAP (1):

2

2

QAPR1

minimizeR
subject to22

tr(V̂ T LQ V̂ R)
GJ (V̂ RV̂ T ) = E00 ,R ⪰ 0,

where LQ ∈ S n +1 , V̂ ∈ Rn +1×(n−1) +1 , E00 ∈ S n +1 is the matrix whose (1, 1)-element is one and all2other elements are zero, J ⊆ {(i, j) | 1 ≤ i, j ≤ n2 } are given, R ∈ S (n−1) +1 is the variable matrix, and2
2
GJ : S n +1 → S n +1 is the gangster operator deﬁned by{
Yij (i, j) ∈ J or (j, i) ∈ J,
(GJ (Y ))ij :=(3)0otherwise.
3

More precisely, the matrix LQ is as follows:()
0− 12 vec(C)TLQ :=,− 12 vec(C)
D⊗F

(4)

2

where vec(C) ∈ Rn is the vector formed by stacking the columns of C on top of one another, andA ⊗ B ∈ Rmp×nq is the Kronecker product of A ∈ Rm×n and B ∈ Rp×q . The matrix V̂ is the normalized22
matrix of V̄ ∈ Rn +1×(n−1) +1 given by
)
(
1
0T
,
V̄ := 1n (e ⊗ e) V ⊗ V
where V ∈ Rn×n−1 is a full-rank matrix

(
V :=

)
In−1
,
−eT

and In−1 is the (n − 1) × (n − 1) element matrix. As a result, the matrices V and V̂ satisfy V T e = 0
and V̂ T V̂ = I.
Note that QAPR1 has a relative interior feasible solution since it is obtained by projecting an SD
relaxation problem to the minimal face [21].
In [14], the authors succeeded in deriving upper bounds and lower bounds of the QAP by applying theStandard ADMM to QAPR1 . In what follows, we show how we can solve QAPR1 by using ADMM
according to the descriptions in [14].

3.1

Representation of the QAP for which the Standard ADMM is applicable

In [14], the authors represent QAPR1 as a problem QAPR2 having two variables R and Y by introducing
a new variable Y and adding the equation Y = V̂ RV̂ T :QAPR2

minimizeR,Ysubject to⟨LQ , Y ⟩GJ (Y ) = E00 ,Y = V̂ RV̂ T ,R ⪰ 0.

(5)

Next, they consider the following optimization problem for which the Standard ADMM is applicable:
QAPR3⟨LQ , Y ⟩ + I(R) + I(Y )Y = V̂ RV̂ T ,minimizeR,Ysubject towhere the second and third constraints of QAPR2 are combined in the objective function as the corresponding indicator functions below:{{
0 if R ⪰ 0,
0 if GJ (Y ) = E00 ,I(R) :=
I(Y ) :=
(6)
∞ otherwise,∞ otherwise.

3.2

Update formula of the variables in the Standard ADMM

The augmented Lagrangian function for QAPR3 is given byρLρ (R, Y, Z) := ⟨LQ , Y ⟩ + I(R) + I(Y ) + ⟨Z, Y − V̂ RV̂ T ⟩ + ∥Y − V̂ RV̂ T ∥2F ,2
4

where Z is the dual variable matrix. Using this function, the variables are updated as follows:
Rk+1
Y

k+1

Z

k+1

:= argminR Lρ (R, Y k , Z k ),:= argminY Lρ (Rk=: Z + ρ(Y

k+1

k+1

(7)

k

, Y, Z ),− V̂ R

k+1

T

V̂ ).

(8)
(9)

As shown in [14], the above updates can be explicitly calculated by
Rk+1

=

1
PS+ (V̂ T (Y k + Z k )V̂ ),
ρ

Y k+1

=

1
E00 + GJ C (V̂ Rk+1 V̂ T − (LQ + Z k )),
ρ

(10)

Z k+1

=

Z k + ρ(Y k+1 − V̂ Rk+1 V̂ T ),(11)

where PS+ is the orthogonal projection onto S+ and J C is a set given by
J C := {(i, j) | 1 ≤ i, j, ≤ n2 + 1)} \ J.
In [14], the authors added the constraints 0 ≤ Yij ≤ 1(∀i, j) if 0 ≤ V̂ RV̂ T ≤ 1 are satisﬁed and showedthat adding these constraints has a profound eﬀect in accelerating the convergence of the Standard
ADMM.

3.3

Stopping conditions of the Standard ADMM

The stopping conditions are given by the optimal conditions of the Lagrangian function,
L(R, Y, Z) := ⟨LQ , Y ⟩ + I(R) + I(Y ) + ⟨Z, Y − V̂ RV̂ T ⟩.
Let us deﬁne f (R) := I(R) and g(Y ) := ⟨LQ , Y ⟩ + I(Y ). Then, the Karush-Khun-Tucker conditions of
QAPR3 are given by
Y − V̂ RV̂ T

= O,

(12)

T

∂f (R) − V̂ Z V̂

∋

O,

(13)

∂g(Y ) + Z

∋

O.

(14)

Here, as in Section 3.3 of [4], we call (12) the primal feasibility constraint, and (13) and (14) the dual
feasibility constraints, respectively. We see that condition (14) is always satisﬁed at each iteration(Rk+1 , Y k+1 , Z k+1 ). The update Y k+1 = argminY Lρ (Rk+1 , Y, Z k ) implies that∂g(Y k+1 ) + Z k + ρ(Y k+1 − V̂ Rk+1 V̂ T ) ∋ O,
and the update Z k+1 = Z k + ρ(Y k+1 − V̂ Rk+1 V̂ T ) implies that∂g(Y k+1 ) + Z k+1 ∋ O.
Thus, condition (14), i.e., the dual feasibility of Y , is always satisﬁed. This implies that we only needto consider the primal feasibility (12) of Z and the dual feasibility (13) of R as stopping conditions.
The primal feasibility (12) of Z can be measured by the Frobenius norm of the residual vector rp =∥Y k+1 − V̂ Rk+1 V̂ T ∥F . To measure the dual feasibility (13) of R, we focus on the update formula (7) of
Rk+1 . Since we set Rk+1 := argminR Lρ (R, Y k , Z k ), Rk+1 satisﬁes∂f (Rk+1 ) − V̂ T Z k V̂ + ρ(Rk+1 − V̂ T Y k V̂ ) ∋ O,
5and this implies that∂f (Rk+1 ) − V̂ T Z k+1 V̂ + V̂ T Z k+1 V̂ − V̂ T Z k V̂ + ρ(Rk+1 − V̂ T Y k V̂ ) ∋ O,
and hence,∂f (Rk+1 ) − V̂ T Z k+1 V̂ ∋ ρV̂ T (Y k − Y k+1 )V̂ .
This implies that if ρV̂ T (Y k − Y k+1 )V̂ = O holds, the dual feasibility (13) of R is guaranteed, andhence, the Frobenius norm of the matrix rd = ∥ρV̂ T (Y k − Y k+1 )V̂ ∥F can be considered as the residual
value of the dual problem.
If the values of rp and rd at the iterate (Rk+1 , Y k+1 , Z k+1 ) are suﬃciently small, we stop the Standard
ADMM. In fact, in [14], the authors chose 10−5 or 10−12 as tolerances, for rp and rd , respectively,
whereby if rp and rd become smaller than these tolerances, the update is stopped.

4

Centering ADMM

We propose a new algorithm, called Centering ADMM, to solve the SDP relaxation problem of QAP bycombining the path-following scheme employed by the interior point methods with the Standard ADMM
described in the previous section.
As in the interior point method, we incorporate a barrier function term with a barrier parameter µ > 0
in the objective function of problem QAPR2 (5) as follows:
BQAPR2⟨LQ , Y ⟩ − µlog(det(R))
GJ (Y ) = E00 ,Y = V̂ RV̂ T ,R ≻ O.

minimizeR,Ysubject to(15)

We consider the following problem for which ADMM is applicable:
BQAPR3minimizeR,Ysubject to⟨LQ , Y ⟩ − µlog(det(R)) + I(Y )Y = V̂ RV̂ T .

(16)

In what follows, we derive the update formula for solving the BQAPR3 with ADMM. First, let us consider
the following augmented Lagrangian function:
ρ
LBQAP(R, Y, Z) := ⟨LQ , Y ⟩ − µlog(det(R)) + I(Y ) + ⟨Z, (Y − V̂ RV̂ T )⟩ + ∥Y − V̂ RV̂ T ∥2F .
ρ
2

(17)

The updating formulas for Y and Z are the same as in the Standard ADMM, deﬁned by (8) and (9),and have the explicit forms (10) and (11). On the other hand, the updating formula for R is deﬁned by
Rk+1

:=
=

argminR LBQAP(R, Y k , Z k )ρ{}ρargminR −µlog(det(R)) + ⟨Z, (Y − V̂ RV̂ T )⟩ + ∥Y − V̂ RV̂ T ∥2F .
2(18)

We can easily check that the function F (R) = −µlog(det(R)) + ⟨Z, (Y − V̂ RV̂ T )⟩ + ρ2 ∥Y − V̂ RV̂ T ∥2F isnnstrictly convex for any R ∈ S++, and hence, F (R) has a unique minimum solution in S++. Using thisresult, we obtain the following proposition.
Proposition 4.1. Suppose that we obtain a spectral decomposition of the matrix V̂ T Z V̂ + ρV̂ T Y V̂ into
an orthogonal matrix P and a diagonal matrix D, asV̂ T Z V̂ + ρV̂ T Y V̂ = P DP T .
6(19)

Then, the new iterate Rk+1 in (18) is given by
Rk+1 = P R̄P T ,where the matrix R̄ is a diagonal matrix whose elements are√
2 + 4ρµDii + Dii
R̄ii =(i = 1, 2, . . . , n).
2ρ

(20)

(21)

Proof. Rk+1 can be calculated in a similar way to what is proposed in [4]. Since the function LBQAP(R, Y, Z)
ρ
BQAP
k+1
in (17) is strictly convex, the gradient of Lρ
(R, Y, Z) at Rshould be O, and hence, we have−µR−1 − V̂ T Z V̂ + ρ(−V̂ T Y V̂ + R) = O.
The above equation and the decomposition (19) imply thatρR − µR−1= V̂ T Z V̂ + ρV̂ T Y V̂
= P DP T ,and, by setting R̄ = P T RP , we haveρR̄ − µR̄−1 = D.

Thus, it turns out that R̄ is a diagonal matrix, and R̄ is given by (21), since R and R̄ should be positive
semideﬁnite.
The update formulas of Centering ADMM consist of (20) for R, (10) for Y and (11) for Z.
For the Standard ADMM, it was shown in [4] that the following dynamic update of the penalty parameterρ, depending on the residual values rp and rd , τ incr > 1, τ decr > 1, and θ > 1, is an eﬃcient way toaccelerate convergence:
incr krp > θrd ,
τ ρ
k+1
k
decr
(22)
ρ
= ρ /τ
rd > θrp , kρotherwise.
We use this update with τ incr = τ decr = 2 and θ = 10 for both Centering ADMM and the Standard
ADMM.
If the residual values of the primal and dual problems become smaller than 0.1 at iteration k, we consider
that the point (X k , Y k , Z k ) is suﬃciently close to the central path, and update the barrier parameterµk by µk+1 = 0.75µk , where the ratio 0.75 was determined from experience.
If the barrier parameter µk is suﬃciently small, e.g., µk < 10−3 , we consider that the point (X k , Y k , Z k )is suﬃciently close to the set of optimal solutions. In that case, the centering eﬀect of the barrier function
is not needed, and we can switch to the Standard ADMM instead of Centering ADMM.
The full description of Centering ADMM is in Algorithm 1.

5

Global convergence of Centering ADMM

In section 3.2 of [4], it has been shown that ADMM (with a ﬁxed penalty parameter ρ > 0) in section 2
has global convergence properties if the following assumptions hold (see also Appendix A of [4]):7

Algorithm 1 Centering ADMM
1: initialization

Y 0 , Z 0 , µ0 , ρ02: while µk < 10−3 do3:Compute Rk+1 = argminR LBQAP(R, Y k , Z k ) by (20)
ρk
k+1
4:
Compute Yby (10)5:Compute Z k+1 by (11)6:rp = ∥Y k+1 − V̂ Rk+1 V̂ T ∥F
7:
rd = ∥ρk V̂ T (Y k − Y k+1 )V̂ ∥F
8:
if rp > 10rd then9:
ρk+1 = 2ρk
10:
else
11:
if rd > 10rp then12:
ρk+1 = ρk /2
13:
else
14:
ρk+1 = ρk15:
end if16:
end if17:
r = max (rp , rd )18:
if r < 0.1 then19:
µk+1 = 0.75µk
20:
else
21:
µk+1 = µk
22:
end if
23: end while
24: Start the Standard ADMM with the initial point (Y k , Z k )8Assumption 1. The (extended-real-valued) functions f : Rn → R ∪ {+∞} and g : Rn → R ∪ {+∞}are closed, proper, and convex.
Assumption 2. The unaugmented Lagrangian L0 (the augmented Lagrangian Lρ with ρ = 0) has asaddle point.
Centering ADMM is an ADMM for problem BQAPR3 (16), and the functions f and g are given by
f (R) := −µlog(det(R)), g(Y ) := ⟨LQ , Y ⟩ + I(Y ),where µ > 0 and LQ is given by (4), and the augmented Lagrangian function LBQAP(R, Y, Z) is deﬁnedρby (17). From the deﬁnition (6) of I(Y ), I(Y ) is the indicator function for an aﬃne space, and hence,
we can see that Assumption 1 holds for Centering ADMM. Thus, if Assumption 2 holds, i.e., the(R, Y, Z) has a saddle point, then Centering ADMM (with nounaugmented Lagrangian function LBQAP
0update of ρ , i.e., ρk = ρ > 0 for every k = 0, 1, . . .) has global convergence properties. In what follows,
we show that the function LBQAP(R, Y, Z) has a saddle point (Proposition 5.1).
0
The Karush-Khun-Tucker conditions of BQAPR2 (15) turn out to beV̂ T (LQ + A)V̂ R = µI,
GJ (V̂ RV̂ T ) = E00 ,2

A ∈ SJn +1 ,R ⪰ O,V̂ T (LQ + A)V̂ ⪰ O,
where

2

2

SJn +1 = {X ∈ S n +1 | GJ (X) = X},(23)

and by setting A = −(LQ + Z), we obtain the following system:(−V̂ T Z V̂ )R = µI,(24)

GJ (V̂ RV̂ T ) = E00 ,(25)

2−(LQ + Z) ∈ SJn +1 ,(26)

R ⪰ 0,(27)−V̂ Z V̂ ⪰ O.

(28)

T

For any ﬁxed µ > 0, the objective function of problem BQAPR2 is strictly convex, and hence the above
system has a unique solution. Proposition 5.1 guarantees that the solution is a saddle point of thefunction LBQAP(R, Y, Z).
0
Proposition 5.1. For any ﬁxed µ > 0, the unique solution (R∗ , Z ∗ ) of the system (24)-(28) satisﬁes
max LBQAP(R∗ , Y ∗ , Z) = LBQAP(R∗ , Y ∗ , Z ∗ ) = min LBQAP(R, Y, Z ∗ ),
0
0
0
Z

R,Y

where Y ∗ = V̂ R∗ V̂ T . That is, (R∗ , Y ∗ , Z ∗ ) is a saddle point of the function LBQAP(R, Y, Z).
0
Proof. From the deﬁnition (6) of I(Y ), (25) and setting Y ∗ = V̂ R∗ V̂ T , we can easily see that{}max LBQAP(R∗ , Y ∗ , Z) = max ⟨LQ , Y ∗ ⟩ + I(Y ∗ ) − µlog(det(R∗ )) + ⟨Z, (Y ∗ − V̂ R∗ V̂ T )⟩0Z

Z

= max {⟨LQ , Y ∗ ⟩ − µlog(det(R∗ ))}Z= LBQAP(R∗ , Y ∗ , Z ∗ ).
0
9Let us show that LBQAP(R∗ , Y ∗ , Z ∗ ) = minR,Y LBQAP(R, Y, Z ∗ ) holds. Before doing so, we deﬁne the00
set Y as follows:2
Y := {Y ∈ S n +1 | GJ (Y ) = E00 }.
(29)
Then, from the deﬁnition (6) of I(Y ), we see that{}min LBQAP(R, Y, Z ∗ ) = min ⟨LQ , Y ⟩ + I(Y ) − µlog(det(R)) + ⟨Z ∗ , (Y − V̂ RV̂ T )⟩
0
R,Y
R,Y
{
}
= min ⟨LQ + Z ∗ , Y ⟩ + I(Y ) − µlog(det(R)) + ⟨Z ∗ , −V̂ RV̂ T ⟩R,Y{}= min {⟨LQ + Z ∗ , Y ⟩ + I(Y )} + min −µlog(det(R)) + ⟨Z ∗ , −V̂ RV̂ T ⟩
Y
R
{
}
∗
(30)
= min {⟨LQ + Z , Y ⟩} + min −µlog(det(R)) + ⟨Z ∗ , −V̂ RV̂ T ⟩ .
Y ∈Y

R

2

Since Z ∗ satisﬁers (26), the deﬁnition (23) of SJn +1 , the deﬁnition (3) of GJ (·), and the deﬁnition (29)of Y imply thatmin {⟨LQ + Z ∗ , Y ⟩} =

Y ∈Y

=
=

min {⟨GJ (LQ + Z ∗ ), Y ⟩}Y ∈Y

min {⟨LQ + Z ∗ , GJ (Y )⟩}Y ∈Y

min {⟨LQ + Z ∗ , E00 ⟩}Y ∈Y= ⟨LQ + Z ∗ , E00 ⟩.

(31)

For any ﬁxed µ > 0, the function −µlog(det(R)) + ⟨Z ∗ , −V̂ RV̂ T ⟩ is strictly convex at any R ≻ O, andthe second term of (30) has a unique minimum solution R satisfying−V̂ T Z ∗ V̂ − µR−1 = 0.
Thus, the fact that R∗ satisﬁes (24) implies that R∗ is the minimum solution of the second term of (30)
and we have{}min −µlog(det(R)) + ⟨Z ∗ , −V̂ RV̂ T ⟩ = −µlog(det(R∗ )) + ⟨Z ∗ , −V̂ R∗ V̂ T ⟩.
(32)
R

Equations (30), (31), and (32) imply thatmin LBQAP(R, Y, Z ∗ ) =
0
R,Y

{
}
min {⟨LQ + Z ∗ , Y ⟩} + min −µlog(det(R)) + ⟨Z ∗ , −V̂ RV̂ T ⟩

Y ∈Y

R

∗

= ⟨LQ + Z , E00 ⟩ − µlog(det(R∗ )) + ⟨Z ∗ , −V̂ R∗ V̂ T ⟩.

(33)

By a discussion similar to derive (31), we also see that
LBQAP
(R∗ , Y ∗ , Z ∗ ) =0⟨LQ , Y ∗ ⟩ + I(Y ∗ ) − µlog(det(R∗ )) + ⟨Z ∗ , (Y ∗ − V̂ R∗ V̂ T )⟩= ⟨LQ + Z ∗ , Y ∗ ⟩ − µlog(det(R∗ )) + ⟨Z ∗ , −V̂ R∗ V̂ T ⟩= ⟨LQ + Z ∗ , E00 ⟩ − µlog(det(R∗ )) + ⟨Z ∗ , −V̂ R∗ V̂ T ⟩.

(34)

(R, Y, Z ∗ ) holds.
(R∗ , Y ∗ , Z ∗ ) = minR,Y LBQAP
Therefore, (33) and (34) guarantee that LBQAP006

Numerical experiments

We conducted numerical experiments to examine the performance of Centering ADMM in comparison
with the Standard ADMM on the QAPLIB instances with symmetric matrices in [2, 3]. We used
10

MATLAB R2018b on an Intel (R) Core (TM) i7-6700 CPU @ 3.40GHz 3.41GHz machine. For the sakeof limiting the computational time, we only dealt with instances of size n ≤ 40.
We set the initial points and the accuracy parameters of Centering ADMM and the Standard ADMM,as
Y 0 = I, Z 0 = −I, µ0 = 1, ρ0 = n, ϵr = 0.1
for all instances.
For both methods, we limited the number of iterations to 10000 and outputted the obtained lower bounds
every 100 iterations. Figures 1 – 7 are plots of the diﬀerence between the lower bounds obtained by thetwo methods, i.e., (the value of the lower bound obtained by Centering ADMM) - (the value of the lower
bound obtained by the Standard ADMM) every 100 iterations for each class of instances.
The horizontal axis shows the number of iterations, and the vertical axis shows the diﬀerence between
the obtained lower bounds. If the diﬀerence is positive (negative), it means that Centering ADMM (theStandard ADMM) computes a better lower bound.
The results allow us to make the following observations for each class of instances.

6.1

Observations for the class of “chr” instances (Figure 1)

Except for instance “hr25a,” the diﬀerence is positive when the number of iterations becomes larger than
1000, which implies that Centering ADMM obtains a better lower bound than the Standard ADMM ateach iteration for most instances of this class.

6.2

Observations for the class of “Had” instances (Figure 2)

For “Had” instances, a signiﬁcant increase in the lower bound occurs during the ﬁrst few iterations, andthe Standard ADMM obtains better results than Centering ADMM for every instance.
Figure 9 shows the results of solving “Had12” by using the Standard ADMM as the diﬀerence between thelower bound and the optimal value every 100 iterations. We can see that the lower bound is suﬃciently
close to the optimal value at 300 iterations.
Table 1 compares the results obtained by the Standard ADMM and those obtained by Centering ADMM.
The table lists the problem name (Prob.), its optimal value (Opt.), and for each ADMM, the pair of thelower bound (LB) and the number of iterations (#Iter) for which the diﬀerence from the optimum value
becomes less than or equal to 0.5.
In every case, the lower bound is suﬃciently close to the optimal value within 2000 iterations, and thissuggests that the centering eﬀect is not required.

6.3

Observations for the class of “Kra” instances (Figure 3)

The diﬀerence is positive once the number of iterations becomes larger than 1100, which implies that
Centering ADMM obtains a better lower bound than the Standard ADMM at each iteration for all
instances of this class.

11

6.4

Observations for the classes of “Rou” and “Scr” instances (Figure 4)

The diﬀerence is always positive, which implies that Centering ADMM obtains a better lower boundthan the Standard ADMM at each iteration for all instances of these classes.

6.5

Observations for the class of “Nug” instances (Figure 6)

Similary to the results for the “Had” instances, a signiﬁcant increase in the lower bound occurs at the
beginning of the iterations for the “Nug” instances, and the Standard ADMM obtains better results thanCentering ADMM for every instance. A diﬀerence from the results for the “Had” instances is that theupper bound is not attained at 10,000 iterations. At an early stage, the lower bound increases rapidly
to a certain value, but after that, the increase becomes quite small.
Table 2 lists the diﬀerences between the lower bounds obtained by the Standard ADMM and by Centering
ADMM. The table shows the problem name (Prob.), the value at which the increase in the ratio of thediﬀerences starts to slow down (Slow Down LB), the number of iterations at which the diﬀerence becomesless than or equal to 1 (Small Diﬀ. #Iter.), and the lower bound obtained by the Standard ADMM at
10000 iterations (LB at 10000 Iter.). We omit the lower bound obtained by Centering ADMM since itis quite close to the value of “LB at 10000 Iter.”Table 3 compares the results obtained by the Standard ADMM and by Centering ADMM. The tableshows the problem name (Prob.), its optimal value (Opt.), and for each ADMM, the pair of the lower
bound (LB) and the number of iterations for which the increase in the ratio of the lower bound values
starts to slow down (#Iter.).
Similarly to the results for the the “Had” instances, the lower bound is suﬃciently close to the optimalvalue within 3000 iterations, and this this suggests that the centering eﬀect is not required in any of thecases.

6.6

Observations for the class of “Tai-a” instances (Figure 7)

Except for instance “Tai35a,” Centering ADMM obtains a better lower bound than the Standard ADMM
at almost every iteration for all instances of this class.

6.7

Observations for the classes of “Els19” and “Tho30” instances (Figure
8)

For instance “Els19,” Centering ADMM obtains better lower bounds than the Standard ADMM. For
instance “Tho30,” the diﬀerence is positive at almost every iterations. Thus, Centering ADMM is better
than the Standard ADMM for this instance as well.

7

Concluding remarks

We devised a new method for solving a semideﬁnite (SD) relaxation of the quadratic assignment problem(QAP), called Centering ADMM. Centering ADMM is an alternating direction method of multipliers(ADMM) combining the centering steps used in the interior-point method. The ﬁrst stage of Centering

12

ADMM updates the iterate such that it approaches the central path by incorporating a barrier function
term in the objective function, as in the interior-point method. If the current iterate is suﬃciently close
to the central path with a suﬃciently small value of the barrier parameter, the method then proceeds
to the Standard ADMM. We showed that Centering ADMM (not employing a dynamic update of thepenalty parameter) has global convergence properties. To observe the eﬀect of the centering steps, we
conducted numerical experiments with SD relaxation problems of the instances in QAPLIB [3]. The
results demonstrate that Centering ADMM is quite eﬃcient for some instances, e.g., all instances in“chr,” “Kra,” “Rou,” and “Scr, ” and instances “Els19” and “Tho30.”
Our future research will include further discussions on convergence of Centering ADMM, providing a
way of determining valid initial parameters, and an extension of the method to general semideﬁnite
programs.

References[1] Anstreicher, K.: Recent advances in the solution of quadratic assignment problem. Mathematical
Programming Series B 97(1-2), 27-42 (2003).
[2] Burkard, R. E., Karisch, S. E., Rendl, F.: QAPLIB–a quadratic assignment problem library. Journalof Global optimization 10(4), 391-403 (1997).
[3] R.E. Burkard, S.E. Karisch and F. Rendl. QAPLIB: A Quadratic Assignment Problem Library.
http://www.seas.upenn.edu/qaplib/[4] Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J.: Distributed optimization and statistical
R in Machinelearning via the alternating direction method of multipliers. Foundations and Trends⃝
learning 3(1), 1-122 (2011).
[5] Cela, E.: The Quadratic Assignment Problem: Theory and Algorithms. Kluwer (1998).
[6] De Klerk, E., Sotirov, R..: Exploiting group symmetry in semideﬁnite programming relaxations of
the quadratic assignment problem. Mathematical Programming 122(2), 225–246 (2010).
[7] Edwards, C.S.: A branch and bound algorithm for the Koopmans-Beckmann quadratic assignment
problem. In: Rayward-Smith V.J. (eds) Combinatorial Optimization II. Mathematical Programming
Studies, 13, 35–52 (1980), Springer, Berlin, Heidelberg.
[8] Gabay, D., Mercier, B.: A dual algorithm for the solution of nonlinear variational problems via ﬁnite
element approximation. Computers & Mathematics with Applications 2, 17-40 (1976).
[9] Glowinski, R., Marroco, A.: Surl’approximation, par éléments ﬁnis d’ordre un, et la résolution, parpénalisation- dualité d’une classe de problémes de Dirichlet non linéaires. ESAIM: Mathematical
Modelling and Numerical Analysis - Modélisation Matématique et Analyse Numérique, 9(R2), 41-76(1975).
[10] Koopmans, T. C., Beckmann, M.: Assignment problems and the location of economic activities.
Econometrica 25, 53-76 (1957).
[11] Lawler, E. L.: The quadratic assignment problem. Management Science 9, 586-599 (1963).
[12] Lin. T., Ma. S., Ye. Y., Zhang, S.: An ADMM-based interior-point method for large-scale linear
programming. arXiv:1805.12344 [math.OC] (2018).
[13] Luo, Z.Q., Ma, W.K., So, A.M., Ye, Y., Zhang, S.: Semideﬁnite relaxation of quadratic optimization
problems and applications. IEEE Signal Processing Magazine 27(3), 20-34 (2010).

13

[14] Oliveira, D. E., Wolkowicz, H., Xu, Y.: ADMM for the SDP relaxation of the QAP. Mathematical
Programming Computation 10(4), 631-658 (2018).
[15] Pardalos, P., Wolkowicz, H. (eds.).: Quadratic assignment and related problems. American Mathematical Society, Providence, RI, 1994. Papers from the workshop held at Rutgers University, New
Brunswick, New Jersey, May 20–21 (1993).
[16] Pardalos, P., Rendl, F., Wolkowicz, H. : The quadratic assignment problem: A survey and recent developments. In Quadratic Assignment and Related Problems, Pardalos, P., Wolkowicz, H.
(eds) DIMACS Series in Discrete Mathematics and Theoretical Computer Science 16, 1-42, AMS,
Providence, RI (1994).
[17] Pong, T.K., Sun, H., Wang, N., Wolkowicz, H.: Eigenvalue, quadratic programming, and semidefinite programming relaxations for a cut minimization problem. Computational Optimization andApplication 63(2), 333–364 (2016).
[18] Povh, J., Rendl, F.: Copositive and semideﬁnite relaxations of the quadratic assignment problem.
Discrete Optimization 6(3), 231–241 (2009).
[19] Rendl, F., Sotirov, R.: Bounds for the quadratic assignment problem using the bundle method.
Mathematical Programming Series B 109, 505–524 (2007). DOI 10.1007/s10107-006-0038-8[20] Sahni, S., Gonzalez, T.: P-complete approximation problems. Journal of the ACM (JACM) 23(3),555-565 (1976).
[21] Zhao, Q., Karisch, S. E., Rendl, F., Wolkowicz, H. : Semideﬁnite programming relaxations for thequadratic assignment problem. Journal of Combinatorial Optimization 2(1), 71-109 (1998).

14




































(a) chr12a




(b) chr15a(c) chr18a







































(d) chr20a(e) chr22a
(f) chr25a





































(g) chr12b




(h) chr15b
(i) chr18b



















































(j) chr20b




(k) chr22b






























(m) chr15c(n) chr20c

Figure 1: Results for “chr” instances15(l) chr12c



















































(a) Had12(b) Had14(c) Had16























(d) Had18(e) Had20

Figure 2: Results for “Had” instances













































(a) Kra30a

(b) Kra30b(c) Kra32

Figure 3: Results for “Kra” instances








































(a) Rou12

(b) Rou15(c) Rou20

Figure 4: Results for “Rou” instances













































(a) Scr12(b) Scr15

Figure 5: Results for “Scr” instances
16(c) Scr20














































(a) Nug12(b) Nug14
(c) Nug15













































(d) Nug16a
(e) Nug16b(f) Nug17





















































(g) Nug18
(h) Nug20(i) Nug21




















































(j) Nug22






(k) Nug24
(l) Nug25































(m) Nug27

(n) Nug28

Figure 6: Results for “Nug” instances

17(o) Nug30













































(a) Tai12a



(b) Tai15a

(c) Tai17a











































(d) Tai20a

(e) Tai25a(f) Tai30a















(g) Tai35a

Figure 7: Results for “Tai-a” instances























(a) Els19(b) Tho30

Figure 8: Results for “Els19” and “Tho30” instances























Figure 9: Diﬀerence between the lower bounds for “Had12”

18

Prob.
Had12
Had14
Had16
Had18
Had20

Opt.
1652
2724
3720
5358
6922

Centering LB
1651.666717
2723.602388
3719.512521
5357.509444
6921.514997

Centering #Iter.
1400
2200
2000
5400
3800

Standard LB
1651.918164
2723.683821
3719.322802
5357.523159
6921.599327

Standard #Iter.
300
500
400
1800
1400

Table 1: Comparison of results for “Had” instances

Prob.
Nug12
Nug14
Nug15
Nug16aNug16b
Nug17
Nug18
Nug20
Nug21
Nug22
Nug24
Nug25
Nug27
Nug28
Nug30

Slow Down LB
567
1009
1140
1598
1217
1706
1892
2505
2380
3525
3399
3624
5125
5023
5947

Small Diﬀ. #Iter.
2000
2700
3000
3200
2900
3400
3300
4200
5100
8300
5100
5500
9300
7300
8700

LB at 10000 Iter.
567.991010.111140.56
1599.27
1218.25
1707.10
1893.532506.33
2381.893528.53
3401.053625.855129.365025.565949.43

Table 2: Comparison of results for “Nug” instances (1)

Prob.
Nug12
Nug14
Nug15
Nug16aNug16b
Nug17
Nug18
Nug20
Nug21
Nug22
Nug24
Nug25
Nug27
Nug28
Nug30

Opt.
578
1014
1150
1610
1240
1732
1930
2570
2438
3596
3488
3744
5234
5166
6124

Centering LB
567.0016543
1009.018699
1140.052527
1598.013623
1217.009463
1706.06183
1892.09605
2505.029389
2380.015468
3525.102227
3399.109253
3624.029764
5125.1364
5023.01236
5947.069915

Centering #Iter.
2100
2600
3500
3000
2700
3400
3000
3900
4200
5000
4300
4700
5800
5700
6600

Standard LB
567.0206964
1009.047161
1140.04133
1598.175632
1217.175222
1706.139849
1892.124008
2505.033302
2380.053905
3525.137032
3399.234062
3624.147168
5125.054708
5023.162624
5947.159729

Table 3: Comparison of results for “Nug” instances (2)19Standard #Iter.
500
600
1100
1200
1000
1400
1200
1400
1800
2400
1900
1900
2500
2900
1700