A semidefinite programming based branch-and-bound framework for the quadratic assignment problem

(1)

Tilburg University

A semidefinite programming based branch-and-bound framework for the quadratic

assignment problem

Truetsch, U.

Publication date:

2014

Document Version

Publisher's PDF, also known as Version of record

Link to publication in Tilburg University Research Portal

Citation for published version (APA):

Truetsch, U. (2014). A semidefinite programming based branch-and-bound framework for the quadratic assignment problem. CentER, Center for Economic Research.

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal

Take down policy

(2)

branch-and-bound framework for the

quadratic assignment problem

by

Uwe Truetsch

(3)

(4)

branch-and-bound framework for the

quadratic assignment problem

Proefschrift

ter verkrijging van de graad van doctor aan Tilburg University op gezag van de rector magnicus prof. dr. Ph. Eijlander in het openbaar te verdedigen ten overstaan van een door het college voor promoties aangewezen commissie in de aula van de Universiteit op vrijdag 31 oktober 2014

om 10.15 uur door Uwe Truetsch,

(5)

Overige Leden prof. dr. M.F. Anjos prof. dr. ir. E.R. van Dam prof. dr. M. Laurent prof. dr. F. Liers

(6)

Let me explore the motivation of this thesis with a short personal anecdote concerning the start of my academic life as a Ph.D. student at Tilburg University. At my very rst ocial working day, I headed for the main building of Tilburg University, the K building, in order to search my new oce. Entering the building, there was no doubt anymore whence the short cut K origins: A highly visible commemorative plaque of Tjalling C. Koopmans with his achievements welcomes me in the entree hall; my new academic home for the upcoming years, the Koopmans building.

Born in 1910 in 's-Graveland, North Holland, the mathematician and economist Tjalling C. Koopmans crowned his career in 1975 when he was announced to be the jointly winner of the Nobel Memorial Prize in economic sciences. Amongst other seminal academic work in his life, he introduced together with his colleague Martin Beckmann [85] the quadratic assignment problem in 1957. Since that time, researchers in the operation research community tried to ght the hardness and complexity of these facility layout problems in order to solve them to optimality. Plenty of (at that time) unsolved problems applied to the quadratic assignment problem by Koopmans and Beckmann were posed to the community and various new methods and techniques were introduced in the last decades.

Still today, more than 50 years later, some of these instances have still remained unsolved. Therefore, inspired by Tjalling C. Koopmans and his commemorative plaque at the main entry of Tilburg University that I have been passing every morning ever since, I became part of the community that tries to improve and invent new exact methods in order to attack the quadratic assignment problem.

Almost 4 years later, this dissertation is the result of hard but passionate work in dealing with the quadratic assignment problem. Working towards the completion of my thesis has been extremely challenging and would be impossible without the support of many important people to whom I would like to take the opportunity to express my gratitude.

First of all, I would like to express my sincere thanks to my family. Growing up, surrounded by such lovely people who guided me through my life to this milestone is priceless. Above all, I would like to thank my parents who supported me in every situation and always believed in me. They taught me never to despair, always to think positive and to ght for one's dreams as they have exemplied it to me through their own life. I am grateful for this endless support, their generous love and their good education.

(7)

who has always been a role model. Therefore, it is not surprising that I followed her way and enrolled for mathematical studies at University of Cologne.

In addition to a lot of other great friends who I have met and who have made my life enjoyable and special ever since, I would like to thank one very special person in particular: Diana, my personal source of energy. You came into my life, relieved me from a hole of listlessness and instilled enthusiasm in me again. It is not easy to express how much I feel blessed to have such wonderful friends and such a great family surrounding me.

However, it is not only the circle of friends and family that guided me to this point. I belong to the people who are lucky to meet incredible personalities and excellent teachers throughout my academic life. In particular, I can honestly say that I was blessed with two ideal supervisors and great persons, Renata Sotirov and Etienne de Klerk. I would like to express to both of you my most sincere thanks for all your patience and support. Next to the endless help and hours in front of the computer in order to debug my source code, Renata also oered me support during some crisis that probably rise up in most Ph.D.'s lives. It was a pleasure to work with you and to be supervised by you.

Further, I also would like to thank the members of my Ph.D. committee: prof. dr. M.F. Anjos, prof. dr. ir. E.R. van Dam, prof. dr. M. Laurent, prof. dr. F. Liers and prof. dr. ir. D. den Hertog.

Finally, I would like to thank Etienne de Klerk. Being supervised by Etienne is the best that could ever happen to me. Next to his extraordinary knowledge and his capability to explain in a comprehensible way, he was always patient, seless and supportive - in other words, the best boss I could have imagined. In addition, thanks to Etienne, I obtained the possibility of an once-in-a live time experience by living at NTU campus in Singapore.

I am so thankful to be blessed with such wonderful friends, teachers and family. Uwe Truetsch

(8)

1 Introduction 1

1.1 The quadratic assignment problem (QAP) . . . 1

1.1.1 Introduction to the QAP . . . 1

1.1.2 Historical example by Krarup . . . 1

1.1.3 Reformulations of the QAP . . . 2

1.1.4 Special case: the traveling salesman problem . . . 4

1.1.5 Complexity of the QAP . . . 6

1.1.6 Symmetry within QAP instances . . . 6

1.2 Existing methods to solve the QAP . . . 7

1.2.1 General overview . . . 7

1.2.2 Branch-and-bound method . . . 8

1.2.3 Existing approaches via branch-and-bound methods . . . 11

1.3 Semidenite programming . . . 12

1.3.1 Symmetry in SDPs . . . 15

1.4 Overview of the thesis . . . 16

I Lower bounds for the QAP 18 2 Semidenite programming relaxations of order n2 ₁₉ 2.1 The QAPR2 and the QAPR3 relaxation . . . 19

2.1.1 The (QAPR2)relaxation . . . 20

2.1.2 The (QAPR3)relaxation . . . 21

2.2 Slater feasible versions . . . 23

2.2.1 The Slater feasible version of the QAPR2 relaxation . . . 29

2.2.2 The Slater feasible version of the QAPR3 relaxation . . . 35

2.3 Symmetry reduction for the QAPR3 relaxation . . . 37

2.4 Numerical results . . . 42

3 Sherali-Adams-type relaxations 45 3.1 The reformulation-linearization technique (RLT) . . . 45

3.2 RLT cuts for the quadratic assignment problem . . . 49

3.2.1 First level RLT relaxation . . . 49

3.2.2 Second level RLT relaxation . . . 50

3.3 QAPR3 as an RLT relaxation . . . 52

3.4 Symmetry reduction for the QAPSDP+RLT2 relaxation . . . 55

(9)

4.3.1 Further variants of the eigenspace relaxation . . . 72

4.4 Symmetry reduction for the eigenspace relaxation . . . 73

4.4.1 Example: the traveling salesman problem . . . 76

4.4.2 Symmetry reduction for selected QAPLIB instances . . . 79

4.5 Relations between SDP bounds . . . 80

4.6.1 Relaxations of order n . . . 89

4.6.2 Symmetry reduced eigenspace relaxation of the TSP . . . 93

4.6.3 All introduced QAP relaxations . . . 94

II Upper bounds for the QAP 100 5 Old versus new heuristics for the QAP 101 5.1 Motivation . . . 101

5.2 Upper bounds for the QAP via heuristics from literature . . . 101

5.2.1 Local search . . . 102

5.2.2 Iterated local search (ILS) . . . 104

5.2.3 Simulated annealing (SA) . . . 104

5.3 A new QP-based heuristic . . . 105

5.3.1 KKT point heuristic . . . 105

5.3.2 A polynomial-time potential reduction algorithm . . . 107

5.3.3 An algorithm for the new QP-based KKT point heuristic . . . 109

5.4.1 KKT point heuristic versus SA and ILS . . . 115

III Complete Branch-and-Bound framework 117 6 An SDP based branch-and-bound framework 118 6.1 Decision tree . . . 118

6.1.1 Subproblem at a node of the decision tree . . . 119

6.1.2 Symmetry at any subnode in the branching tree . . . 121

6.2 Choice of branching decisions . . . 122

6.2.1 Branching rule . . . 122

6.2.2 Node selection strategy . . . 124

6.3 Choice of relaxation . . . 125

6.4 Choice of heuristic . . . 126

6.5 Summary of the branch-and-bound framework . . . 127

(10)

B Details for QAPLIB instances 139

(11)

2.4.1 Numerical (in)stability concerning Example 2.4.1 . . . 43

3.5.1 Dierent bounds on the min k-section of the Higman-Sims graph . . . 65

3.5.2 Dierent bounds on the max k-section of the Higman-Sims graph . . . 65

3.5.3 Dierent bounds on the min/max 11-section of the Cameron graph . . . . 65

4.4.1 Results for selected TSPLIB instances with symmetry. . . 80

4.6.1 Performance of QAPEig , QAPEig+SDRMS and (QAPSDRMS) - part 1 . . 89

4.6.4 New best known lower bounds for several QAPLIB instances . . . 91

4.6.5 Lower bounds on some small TSPLIB instances from various convex relax-ations . . . 93

4.6.6 Results for instances on n = 8 cities, constructed from the 24 facet dening inequalities . . . 94

4.6.7 Comparison of the lower bounds for various dierent QAP relaxations . . 95

4.6.9 Detailed comparison (QAPSDRMS) versus QAPEig - part 1 . . . 96

4.6.10Detailed comparison (QAPSDRMS) versus QAPEig - part 2 . . . 97

4.6.11Detailed comparison (QAPSDRMS) versus QAPEig - part 3 . . . 98

5.4.1 Performance of various rounding approaches . . . 113

5.4.2 Detailed look at performance of Barvinok's rounding approach [12] using dierent starting points . . . 113

5.4.3 Variation in starting points for KKT point heuristic . . . 114

5.4.4 Local search (LS) versus KKT point heuristic . . . 114

5.4.5 Comparison of heuristics: KKT versus SA versus ILS . . . 115

5.4.6 Comparison of heuristics: KKT versus SA versus ILS . . . 116

6.6.1 Detailed look at (BaBcombined) , BaBQAP_R3 ,BaBQAP_Eig and BaBQAP_SDRMS 128 6.6.2 Summary of all non-symmetric QAPLIB instances with n ≤ 16 . . . 130

6.6.3 Detailed comparison for Scr15 . . . 131

B.0.1Details for QAPLIB instances - part 1 . . . 139

(12)

1.2.1 Example of a branch-and-bound tree . . . 9

1.2.2 A general branch-and-bound framework for the QAP . . . 10

4.6.1 Comparison of (QAPSDRMS) and QAPEig on QAPLIB instances of n = 16 92 4.6.2 Comparison of (QAPSDRMS) and QAPEig on QAPLIB instances of 30 ≤ n ≤ 35 . . . 93

4.6.3 Comparison of three SDP bounds on Nugent instances: the relaxation gap 99 4.6.4 Comparison of three SDP bounds on Nugent instances: the CPU time . . 99

6.1.1 Example of a generic branch-and-bound tree for n = 3 . . . 119

6.2.1 Decision tree generated by the depth-rst search algorithm . . . 125

6.2.2 Decision tree generated by the breadth-rst search algorithm . . . 126

(13)

Rn : ndimensional Euclidian vector space; Rn+ Rn++

: nonnegative (strictly prositive) orthant of Rn; Rm×n : space of (m × n) real matrices;

Mij : (i, j)th entry of a matrix M ∈ Rm×n;

M(ij) : the block in the ith row-block and jth column-block of M; MT : the transpose of a matrix M ∈ Rm×n;

Sn×n = M ∈ Rn×nM = MT

: set of all symmetric matrices; M 0 (M 0) : M ∈ Sn×n is a positive semidenite (positive denite) matrix; M 0 (M ≺ 0) : M ∈ Sn×n is a negative semidenite (negative denite) matrix;

Sn×n+ = M ∈ Sn×n|M 0

: set of all positive semidenite matrices; Sn : the symmetric group on {1, ..., n};

Πn : set of all (n × n) permutation matrices;

coli(M ) ∈ Rm : consists of all entries of the ith column of M ∈ Rm×n; In : the (n × n) dimensional identity matrix;

Jn : the (n × n) dimensional all-ones matrix;

e : the all-ones vector; ei : ith standard unit vector;

Eij = eieTj : all-zeros matrix except the (i, j)th entry which equals one;

0n : the all-zeros vector of dimension n;

vec (M) ∈ Rmn _{: refers to the column-wise vectorized matrix M ∈ R}m×n_;

diag (M ) ∈ _Rn : a vector that contains all diagonal elements of M ∈ Rn×n; Diag (v) ∈ _Rn×n : diagonal matrix with iths diagonal entry equals vi for v ∈ Rn;

λi(M ) : ith largest eigenvalue of a matrix M ∈ Sn×n;

det(M ) = Y i λi(M ) : determinant of M ∈ Rn×n; tr(M ) = X i Mii= X i

λi(M ) : the trace of a matrix M ∈ Rn×n;

hM1, M2i = tr(M1M2T) : the trace inner product of M1, M2 ∈ Rm×n;

M1◦ M2 : the component-wise (Hadamard) product of M1, M2 ∈ Rm×n;

M1⊗ M2 : Kronecker product : block matrix with (M1⊗ M2)(ij)= (M1)_ijM2;

(14)

Introduction

1.1 The quadratic assignment problem (QAP)

1.1.1 Introduction to the QAP

The quadratic assignment problem (QAP) was introduced in 1957 by Koopmans and Beckmann [85] to model facility layout problems where we need to assign departments to locations including material ow between the departments. Thus, let us assume in the following to have n facilities and n locations. Let Aij with i, j ∈ {1, ..., n} denote the ow

between every pair (i, j) of facilities and let Bklrepresent the distance between every pair

(k, l)of locations for k, l ∈ {1, ..., n}. We view the product AijBkl as the cost incurred if

facilities i and j are respectively assigned to locations k and l. Further cidkwith c, d ∈ Rn

denotes the cost for allocating facility i to location k. (For example, dk may represent

the cost per worker at location k, and ci the number of workers needed to run facility i.)

The objective is to nd a permutation π that assigns facilities to locations such that the total costs of the assignment is minimal. Note that π (i) = j holds if and only if facility i is assigned to location j. Thus, the QAP can be formulated in its Koopmans-Beckmann formulation: min π∈Sn n X i,j=1 AijBπ(i)π(j)+ n X i=1 cidπ(i), (1.1.1)

with given symmetric data matrices A, B ∈ Sn×n_{, c, d ∈ R}n _{and where S}

n denotes the

symmetric group (of all permutations) on {1, ..., n}.

Since the introduction of the quadratic assignment problem by Koopmans and Beck-mann [85], the QAP has attracted a great deal of attention in literature in terms of theory, applications, and solution methodologies; see e.g. [24, 25, 104].

1.1.2 Historical example by Krarup

(15)

For a hospital, facilities include surgery, radiology, etc, and locations denote the physical locations of the hospital buildings.

Among the criteria that had to be considered in order to build the university hospital, the announcement explicitly asked for the usual QAP objective: nd a layout minimizing the sum CDIST of Communication x DISTance taken over all pairs of facilities to be located, see [73]. Krarup and some colleagues were commissioned to nd a best bound on CDIST or if possible to nd provably optimal solutions.

At that time, no algorithms were capable (not even close) to solve QAP instances of dimension 30 or higher to optimality. Thus, Krarup decided to generate upper bounds based on a randomized heuristic [86] in the hope that he could beat the objective function values reached by the architects. Indeed, no architect managed to come up with a layout beating Krarup's bound minCDIST = 91730. The contractor was satised and Krarup succeeded in his commission.

However, Krarup himself and many other researchers in the eld of combinatorial opti-mization felt challenged to improve that bound or to prove optimality. So, the Krarup30a instance became a standard benchmark problem. In the subsequent years, a number of heuristic approaches were developed for the QAP in order to improve Krarup's bound. In 1976, Burkard and Stratmann [29] managed to reduce CDIST to 90420 and Krarup himself proposed subsequently a cyclic rotation on the latter solution that yields a new improved bound CDIST = 89100, see Chapter 3 in [73] for details. In the 1980s, not much happened regarding the improvements for the Krarup30a bound until 1990 when Skorin-Kapov [116] found an upper bound of 88900 using a strict tabu search. Since the heuristic did not guarantee optimality, the question of a better solution still remained.

Finally, after 27 years, two approaches provided multiple optimal solutions in 1999. The rst one was an iterated local search algorithm due to Stützle [118] (published only in 2006) that discovered 256 global optimal solutions. The second one was a branch-and-bound enumeration by Hahn et al. [72] (published in 2001) that proved that Stützle's 256 solutions were indeed optimal. According to [73] the branch-and-bound enumeration of the Krarup30a was begun in December of 1998 and ended on April 12, 1999, thus a total run time of 98.6 days. The number of nodes examined in the tree was 29, 764, 589. And the result was as expected: the optimal objective function value of the Krarup30a instance was indeed 88900.

1.1.3 Reformulations of the QAP

Since the introduction of the QAP by Koopmans and Beckmann [85], dierent versions of the quadratic assignment problem have been formulated that are equivalent to (1.1.1). In this context, the so-called trace formulation of the QAP was introduced in order to work with permutation matrices, say Pπ ∈ Rn×n, instead of using permutations π ∈ Sn. Note

that (Pπ)ij = 1 if facility i is assigned to location j, and thus that

(Pπ)ij =

(

1 if π (i) = j

0 otherwise (1.1.2)

(16)

holds. Since the mapping is bijective, it is easy to see that for any permutation matrix Pπ dened by (1.1.2) there exists only one entry equal to one in each row and column of

Pπ, respectively. The set of all permutation matrices, say Πn, can be dened as:

Πn:=    Pπ ∈ {0, 1}n×n n X i=1 (Pπ)ij = 1 ∀j = 1, ..., n; n X j=1 (Pπ)ij = 1 ∀i = 1, ..., n    . (1.1.3) Denition 1.1.1. We dene a doubly stochastic matrix as a square matrix of non-negative real numbers where each of its rows and columns sums to one.

In the following, we state the well-known theorem attributed to the work of Birkho [16] in 1946 and the work of von Neumann [124] in 1953, respectively. The proof is omitted but can be found throughout the literature, amongst others in [80].

Theorem 1.1.2. (Birkho-von Neumann theorem) Every doubly stochastic (n×n) matrix is a convex combination of permutation matrices.

Note that the doubly stochastic (n × n) matrices are precisely the convex hull of Πn. Further, recall that there is a one-to-one correspondence between the set of all

permutations Sn and the set of all (n × n) permutation matrices Πn such that if π ∈ Sn

corresponds to Pπ ∈ Πn according to (1.1.2) then Bπ(i)π(j) = PπTBPπ

ij. Hence, we can

express the quadratic assignment problem (1.1.1) in its equivalent trace formulation: min Pπ∈Πn tr APT πBPπ +tr Diag (c) PπTDiag (d) Pπ (1.1.4) where Diag (c) ∈ Rn×n _{denotes the diagonal matrix with the components of the vector}

c ∈ Rn on its diagonal.

Under the assumption that the input matrices are symmetric, say A, B ∈ Sn×n_{, we}

can derive a third equivalent reformulation of the QAP. To introduce this formulation, we rst recall the denition of the Kronecker product of matrices and some of its properties. Denition 1.1.3. Let A ∈ Rk×l _{and B ∈ R}m×n_{. The Kronecker Product of these two}

matrices, say B ⊗ A, is dened in the following way: B ⊗ A :=    [B11A] · · · [B1nA] ... ... ... [Bm1A] · · · [BmnA]   ∈ R (mk)×(nl)_.

The following properties of the Kronecker product are well-known:

Let A ∈ Rk×l_{, B ∈ R}l×m_{, C ∈ R}n×p_{, D ∈ R}p×q _{be some arbitrarily matrices then}

(A ⊗ C) (B ⊗ D) = AB ⊗ CD . (1.1.5)

For two square matrices A ∈ Rk×k _{and B ∈ R}l×l_{the following property is satised:}

(17)

[37].

Lemma 1.1.4. Let A, B ∈ Sn×n_{, F ∈ R}n×n _{and let P ∈ R}n×n_{. The following equivalence}

illustrates the relationship between the Kronecker product and the vec-operator:

AP B = F ⇔ (B ⊗ A) vec (P ) = vec (F ) (1.1.7) Furthermore, we itemize two operator rules without any proof that are well known: Let M1, M2 ∈ Rn×n be an arbitrarily matrix. Then it follows directly from the

denition of the trace operator that

tr (M1M2) = tr (M2M1) (1.1.8)

is satised.

Let M1, M2 ∈ Rn×n be two arbitrarily matrices and let h·, ·i represents the trace

inner product on Rn×n_{. Therefore,}

hM₁, M2i := tr M1TM2 = vec (M1)T vec (M2) (1.1.9)

is satised.

The upcoming Kronecker product formulation of the quadratic assignment problem follows directly from Lemma 1.1.4 and from the operator rules (1.1.8) and (1.1.9).

Theorem 1.1.5. Let A, B ∈ Sn×n _{and c, d ∈ R}n_{. The quadratic assignment problem in}

its formulation (1.1.4) is equivalent to the problem min

Pπ∈Πn

vec (Pπ)T [B ⊗ A + Diag (d) ⊗ Diag (c)] vec (Pπ) . (1.1.10)

Note that the Kronecker product formulation (1.1.10) of the QAP assumes the input matrices A, B ∈ Sn×n _{to be symmetric. Provable by direct verication, the following}

statement shows that it is enough to assume that at least one of the input matrices A and B has to be symmetric.

Proposition 1.1.6. Let w.l.o.g. ¯A /_{∈ S}n×n _{be a non-symmetric matrix and B ∈ S}n×n

be symmetric. Then we can replace the non-symmetric matrix by its symmetric part A := 1₂ A + ¯¯ AT_{∈ S}n×n without aecting the optimization problem (1.1.10).

In this thesis, we may assume that both input matrices A, B ∈ Sn×n _{are symmetric.}

In the case that one input matrix might not be symmetric, we can apply Proposition 1.1.6 to convert the non-symmetric data matrix.

1.1.4 Special case: the traveling salesman problem

In this section, we want to present another application of the QAP. Therefore, we will refer to the QAP in its trace formulation (1.1.4) without a linear term:

min

X∈Πn

(18)

where A, B ∈ Sn×n_{. One of the most famous applications for the QAP in combinatorial}

optimization is the traveling salesman problem.

Given a list of cities and the respective distances between each pair of cities, the traveling salesman problem searches for the shortest route that visits each city exactly once and returns to the city where we started. Let us consider in the following a complete graph Kn(D)with n vertices where D denotes the matrix of weights which represents the

distances Dij = Dji > 0 between two locations i 6= j for all i, j = 1, ..., n. It is asked to

nd a Hamiltonian cycle of a minimum length in Kn(D). To formulate the TSP, we use

the adjacency matrix Cn of a canonical circuit on n nodes:

Cn=             0 1 0 · · · 0 1 1 0 1 ... 0 0 1 ... ... ... ... ... ... ... ... 1 0 0 ... 1 0 1 1 0 · · · 0 1 0             ∈ Sn×n_.

Now, the TSP is to nd an optimal permutation X ∈ Πn, such that a feasible cycle with

the least distances is found. This, we can model as a special case of the QAP in (1.1.11) in the following way:

min X∈Πn tr 1 2DX T_C nX (1.1.12) where Cn, D ∈ Sn×n are dened as above.

One of the most famous lower bounds for the TSP is the Held-Karp bound [75] whose optimal value coincides with the LP relaxation that was introduced in 1954 by Dantzig, Fulkerson and Johnson [36]. Held and Karp [75] use the fact that a Hamiltonian cycle is 2-connected and describe it by the following subtour elimination inequalities:

X

i∈I,j /∈I

Xij ≥ 2 ∀∅ 6= I ⊂ {1, ..., n} . (1.1.13)

Furthermore, for the all-ones vector e ∈ Rn _{and the all-ones matrix J ∈ R}n×n_{, the}

LP relaxation which equals in its optimal value the Held-Karp bound can be written as follows: (TSPHK)                          min X 1 2tr(DX) s.t. Xe = 2e diag(X) = 0n 0 ≤ X ≤ J X i∈I,j /∈I Xij ≥ 2 ∀∅ 6= I ⊂ {1, ..., n} . (1.1.14)

(19)

Karisch, Rendl, Wolkowicz type [128]. In the upcoming second chapter, we will introduce this relaxation in detail, see the so-called TSPQAPR3

relaxation in (2.3.8) on page 42. Comparing the TSPQAPR3

relaxation to the (TSPHK) relaxation, de Klerk et al. [44]

show that neither their semidenite programming bound dominates the Held-Karp bound nor vice versa (see [44], Theorem 5.1). Nevertheless, in their numerical comparisons, see [39], it turns out that the SDP relaxation seems to oer a more powerful lower bound for some special TSP instances.

1.1.5 Complexity of the QAP

The quadratic assignment problem is well-known to be among the most dicult combi-natorial optimization problems. Speaking in complexity terms, the QAP belongs to the class of computationally hard problems known as NP-complete (since the TSP is a special case of the QAP). The proof that the QAP indeed belongs to the class of NP-complete problems was rst shown in 1976 by Sahni and Gonzalez [113]. Furthermore, Sahni and Gonzalez [113] also proved that nding an approximate solution is also NP-hard. This makes the QAP among the hardest of all combinatorial optimization problems. To fur-ther illustrate this fact, Pardalos, Rendl and Wolkowicz [104] explain how ofur-ther famous NP-hard problems are special cases of the QAP.

However, it takes more than the theoretical complexity to explain why it is so dicult to solve the QAP to optimality. For most of the applications, we can nd many solutions whose value is close to the optimum. So, even when the best solution is obtained, it is very hard to prove its optimality, see [24] for more details. Remark that a straightforward enumeration method for all n! feasible solutions leads to an overwhelming number of permutations if n ≥ 12, say (see Povh [108], Table 1).

Hahn and Krarup [73] give a short overview over the computational development in the last fty years concerning the QAP: In the late 1960s, it was an achievement to nd the optimal solution to a dicult instance of size n = 8. Later, in the 1970s and 80s, one could only expect to solve dicult instances for n < 16. It was not until the mid-1990s that Clausen and Perregaard [32] were able to enumerate the very dicult Nugent20 instance from the QAPLIB. Clearly, much progress has been made since then amongst others by Hahn et al. [70], Anstreicher and Brixius [9, 20] and, recently Fischetti et al. [54] who were able to solve dierent instances of the QAP of size 30 ≤ n ≤ 32 to optimality for the rst time. In general, solving QAPs of size n > 30 is still computationally impractical.

1.1.6 Symmetry within QAP instances

(20)

Let us rst introduce the concept of an automorphism group of a matrix.

Denition 1.1.7. The automorphism group, aut (M), of a matrix M ∈ Rn×n_{is the group}

of permutation matrices Pπ ∈ Πnthat satises

aut (M ) :=Pπ ∈ Πn

M = P_πTM Pπ .

Surprisingly often, QAP instances have data matrices with large automorphism groups (see [45]) such that symmetry reduction is applicable. For further applications in combi-natorial optimization and for more details, we refer to [38].

In the following, we give a crucial lemma that shows how to gain information on optimal assignments based on the symmetry in the input matrices. Lemma 1.1.8 is proven in Appendix A.1.

Lemma 1.1.8. Consider the QAP in its form (1.1.1) without linear term. There are at least as many optimal assignments π ∈ Sn as elements in the automorphism group

aut (A).

The following example illustrates the statement given by Lemma 1.1.8.

Example 1.1.9. For a better understanding of the previously introduced denitions, let us consider the following example of the traveling salesman problem. Therefore, let us consider a 5-cycle where

C5=       0 1 0 0 1 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 1 0 0 1 0      

denotes the adjacency matrix of the standard circuit on 5 vertices. For the 5-cycle denoted by C5, its automorphism group is the so-called dihedral group on 5 elements of the order

10 = |aut (C5)|. Thus, we have at least 10 optimal assignments.

1.2 Existing methods to solve the QAP

1.2.1 General overview

In the eld of combinatorial optimization, and in particular for the quadratic assignment problem, it is an ambitious challenge to prove optimality. Therefore, dierent types of exact algorithms were introduced to provide the global optimal solution to the QAP. Finke, Burkard and Rendl [53] classied them in three dierent types of enumeration procedures that are applied to the QAP:

single assignment algorithms, pair assignment algorithms, cutting plane methods.

(21)

combinatorial optimization problems to optimality, in particular the quadratic assignment problem. As described by Burkard et al. [24], the branch-and-bound algorithm appears to be the most ecient exact algorithm for solving QAPs.

Pair assignment algorithms are branch-and-bound methods in which always pairs of facilities are assigned to pairs of locations. Among others, pair assignment algorithms were developed by Land [89], Gavett and Plyter [58] and Nugent et al. [101] but, according to [53], did not yield successful numerical results, so far.

The class of cutting plane methods can get subdivided into two types; traditional cut-ting plane methods and polyhedral cutcut-ting plane methods, see [24]. The traditional cutcut-ting plane method for the QAP has been developed by Bazaraa and Sherali [14], Kaufmann and Broeckx [83] and others. These algorithms are based on mixed integer linear pro-gramming formulations for the QAP which are suitable for Benders' decomposition, see [15]. In spite of the insightful mathematical approach the computational results were not successful according to [53]. Whereas the polyhedral cutting plane method or also known as branch-and-cut method provides more promising numerical results, in particular for the QAP. Recently, Fischetti, Monaci and Salvagnin [54] achieved impressive numerical results for some specic QAPLIB instances (the highly symmetric ones by Eschermann and Wunderlich [51]) using a branch-and-cut method. The branch-and-cut approach is in-spired by two prior methods; it combines elements of the traditional cutting plane method and the branch-and-bound method. For more details and applications on branch-and-cut methods in combinatorial optimization where semidenite programming is used, see for example Ghaddar, Anjos and Liers [60].

1.2.2 Branch-and-bound method

As mentioned before, branch-and-bound algorithms are the most commonly used exact solution methods for the QAP. Beginning in 1963 with the branch-and-bound approach due to Lawler [91], several other implementations for the QAP have been followed, for example by Graves and Whinston [66], Mautor and Roucairol [99], Clausen and Perregaard [32], Hahn et al. [70], Brixius and Anstreicher [20] and others.

They all use a branch-and-bound algorithm that repeatedly divides a QAP into sim-pler subproblems of the same type which are then either solved to optimality or further subdivided. Note that the root node solution to the relaxation constitutes a lower bound to the original subproblem. If a lower bound at a given node exceeds the value of a previ-ously known discrete solution (dened in the following as the incumbent value, see [19]), then continued search in this subtree cannot lead to an improved solution, and thus the node gets pruned.

In Figure 1.2.1 the nodes Ni with i ∈ N are labeled with their lower bound values,

say opt [R (Ni)], given by an appropriate relaxation R (Ni) of the QAP that corresponds

to the respective node Ni. Let F be any pending (father) node at the kth level of the

tree with its m = n − k child nodes C1, ..., Cm. If a father node is not pruned, then the

(22)

(23)

to a quadratic assignment problem.

Figure 1.2.2: A general branch-and-bound framework for the QAP

Both, the correctness and the eectiveness of a branch-and-bound algorithm critically depends on its two main phases:

The branching phase:

This phase denes how child problems are created. For QAPs, the child problems are generated by assigning an unassigned location to an unassigned facility in each branching step. An eective branching strategy aims to x assignments that result in child problems which prune as quickly as possible.

(24)

Lemma 1.2.1. Let F be any pending node at the kth level (0 ≤ k ≤ n) of the tree with its m = n − k child nodes C1, ..., Cm. Then the following two properties have

to hold for every correct branching approach:

Every feasible solution of F is a feasible solution of at least one of C1, ..., Cm.

Every feasible solution of C1, ..., Cm is a feasible solution of F .

The bounding phase:

The bounding phase deals with lower bounds for the QAP and, in particular, how to compute them. Again, with respect to the notation in Figure 1.2.1, we can describe a valid lower bounding procedure according to [59] if it produces relaxations at any node in the tree, say F , with the following properties:

If R(F ) has no feasible solutions, neither does F .

The minimum value of F is no less than the minimum value of R(F ).

If an optimal solution of R(F ) is feasible in F , then it is an optimal solution of F .

In the following section, we give a detailed overview of existing methods including their explicit choice of a valid lower bounding procedure for the QAP.

1.2.3 Existing approaches via branch-and-bound methods

Although many bounding techniques have been developed for the QAP in the last decades, the Gilmore-Lawler bound (GLB) still represents a benchmark. The reason is the simplic-ity of the linear programming bound due to Gilmore [62] and Lawler [91] which implies its computational eciency. In the subsequent decades, various lower bounds for the QAP have been obtained based on linear programming, amongst others the bound of Carraresi and Malucelli [30], the bound of Assad and Xu [11], the bound of Adams and Johnson [3], and the bound of Hahn and Grant [69]. In 2001, a breakthrough in branch-and-bound methods for more complex classes of relaxations was achieved by Anstreicher and Brixius [9, 20] who successfully used a convex quadratic programming bound (QPB0) for their branch-and-bound approach. This approach yields proven optimality to various QAP in-stances of larger size, among others to the famous Nugent30 instance. In the following, we will give a short introduction for these two lower bounding procedures.

A linear programming bound due to Gilmore [62] and Lawler [91]:

The GLB is well known to be one of the rst to compute bounds for the quadratic assignment problem and can be obtained by relaxing an equivalent integer linear programming formulation of the QAP. The Gilmore-Lawler bound can be computed by solving a linear assignment problem of the following form, see [93]:

(25)

Rn×n excluding the value Aii ∈ R. In an analogous way, we also dene ˆbj ∈ Rn−1

for the given QAP data matrix B ∈ Rn×n_.

A convex quadratic programming bound due to Anstreicher and Brixius [9, 20]: In 2001, Anstreicher and Brixius [9, 20] applied the following convex quadratic programming relaxations within a branch-and-bound framework:

(QPB0)            min X vec(X) T_{Qvec(X) + hC, Xi +}_λ(VT_{AV ), λ(V}T_{BV )} − s.t. Xe = XTe = e X ≥ 0, (1.2.2) where Q := (B ⊗ A) − (I ⊗ S) − (T ⊗ I) ∈ Rn2_×n2 , S = VSVb T, T = VT Vb T and V ∈ Rn×(n−1) represents any matrix with orthonormal columns such that eTV = 0, see [20]. The matrices S, bb T ∈ Rn−1×n−1 are optimal solutions of a semidenite programming problem          max b S, bT tr bS + tr bT s.t. I ⊗ bS + b T ⊗ I VTBV ⊗ VTAV ,

but they can be eciently obtained from the spectral decomposition of the terms VTAV and VTBV, see [9], Section 4. Furthermore,

hx, yi₋:= min π∈Sn n X i=1 xiyπ(i)

denotes the minimal product and λ(M) ∈ Rn_{is dened as the vector of eigenvalues}

of a matrix M ∈ Rn×n_.

In this work, we want to enlarge the portfolio of relaxations used within a branch-and-bound framework to the class of semidenite programming (SDP) relaxations. Thus, we will generate a rst semidenite programming based branch-and-bound algorithm that solves QAPs to optimality. This approach is motivated by the recent development in semidenite programming in particular concerning ecient available solvers as MOSEK2_,

see [6]. Furthermore, SDP relaxations are known to oer qualitatively strong lower bounds, in particular compared to (GLB) and (QPB0).

1.3 Semidenite programming

Previously, we introduced linear and quadratic programming relaxations. In order to improve the approximation quality, stronger relaxation methods have attracted the focus

(26)

of recent research, amongst others semidenite programming (SDP). The idea to construct semidenite relaxations for a combinatorial optimization problem dates back to a seminal work in 1979 when Lovasz [97] introduced an SDP bound for the stability number of a graph. For an overview of SDP relaxations for some combinatorial optimization problems, see for example the work by Sotirov [117].

Semidenite programming can be seen as linear programming over the cone of positive semidenite matrices, denoted by Sn×n

+ .

Denition 1.3.1. A symmetric matrix M ∈ Sn×n _{is called positive semidenite, say}

M 0, if xT_{M x ≥ 0} _{holds for all x ∈ R}n_.

The following characterization of a positive semidenite matrix can be found for ex-ample in Helmberg's survey on semidenite programming [76].

Theorem 1.3.2. ([76], Theorem 2.4) For M ∈ Sn×n _{the following statements are}

equiv-alent:

1. M is positive semidenite. 2. λi(M ) ≥ 0 ∀i = 1, ..., n.

3. There exists a matrix C ∈ Rm×n _{so that M = C}T_C. _{For any such C, rank (C) =}

rank (M ).

4. hM, Di ≥ 0 for all positive semidenite matrices D ∈ Sn×n + .

In the following, we explain the relevant mathematical background for the semidenite programming by presenting its duality theory. Therefore, let us in the following consider the primal version of a standard semidenite programming problem:

(PSDP)        min X0 tr (A0X) s.t. tr (AiX) = b ∀i = 1, ..., k (1.3.1) with Ai ∈ Sn×nfor all i = 0, ..., k are given. By introducing Lagrange multipliers yi ∈ R for

each of the i = 1, ..., k equality constraints in (PSDP), we can derive the dual formulation

of (1.3.1) as follows: (DSDP)              max y∈Rk b T_y s.t. A0− k X i=1 (yiAi) 0. (1.3.2)

For more details on the duality of semidenite programming problems, we may refer the reader to [92]. Analogously to the duality theory in linear programming, the construction of (DSDP)implies that the dual objective value can never exceed the objective value the

(27)

tr (A0X) − bTy ≥ 0. (1.3.3)

If equality holds in (1.3.3) then the the primal dual pair (X, y) represents an optimal solution for (PSDP) and (DSDP), respectively. This status is well known under the name

of strong duality. In contrast to linear programming, remark that it is no longer true that optimality implies equality in (1.3.3), see Example 1.3.5.

Denition 1.3.3. A semidenite programming problem (1.3.1) is Slater feasible if there is a strictly feasible point X 0.

Analogously, the denition also holds for the dual problem (DSDP) with a strictly

feasible solution y. The following strong duality theorem (see e.g. Theorem 2 in [90]) shows that the Slater feasibility condition insures the strong duality. This statement has its origin in the work by Dun [49].

Theorem 1.3.4. Let us denote the optimal solution value of (PSDP)by p∗:= opt (PSDP)

and d∗ _{respectively for (D}

SDP). If (1.3.1) satises the Slater constraint qualication and

p∗ is nite, then p∗= d∗, and the dual inmum is attained.

The following example adapted from [123] demonstrates diculties that may arise for SDP problems that do not meet the Slater constraint qualication.

Example 1.3.5. Let us describe the primal problem by

(PEx)                  min y y1 s.t.    0 y1 0 y1 y2 0 0 0 y1+ 1    0.

By construction, the value y1 is forced to be zero since the top-left entry in the

pos-itive semidenite matrix equals zero. Thus, y1 = 0 is part of any feasible solution

(y1 = 0, y2≥ 0) with a primal optimal value opt (PEx) = 0. Its dual formulation can

be rewritten as (DEx)                  max X∈S3 −X33 s.t. X12+ X21+ X33= 1 X22= 0 X 0.

Fixing X22 = 0 together with the positive semideniteness constraint, it follows that

X12= X21= 0. Thus, X33= 1 yields an optimal solution value of opt (DEx) = −1which

does not coincide with opt (PEx). The reason therefore is the fact that both, primal and

dual, have no strictly feasible solution. Thus, i.e. Slater condition is not satised and Theorem 1.3.4 does not hold.

(28)

Rendl [90], in the work by Todd [120] or by Vandenberghe and Boyd [123]. Among others, the latter authors oer in a subsequent work a proof for the Schur complement theorem which is needed later in this thesis.

Theorem 1.3.6. (Schur complement theorem, e.g. [18], p. 650f.) Consider a matrix M ∈ Sn×n to be symmetric partitioned as M = MA MB M_BT MC , where MA∈ Sk×k. If MA is invertible, the matrix

S = MC− MBTM −1 A MB

is called the Schur complement of MA in M. Furthermore, if MA 0, then M 0 if and

only if S 0.

1.3.1 Symmetry in SDPs

The idea of symmetry reduction for instances that consist of group symmetric data, say with large automorphism groups, is especially useful when solving semidenite programs. Exploiting structure in the SDP data is much more dicult than in the LP case. While i.e. sparsity in LP can be easily exploited this does not hold in the case of SDP. Thus, the methodology of symmetry reduction in SDP gained more and more interest, especially from the 1990s onwards (due to availability of interior point methods). Recently, two detailed surveys about symmetry reduction in SDP were published by Vallentin [121] and de Klerk [38]. For further reading, we might refer to the Appendix A.2 and also to the work by Gatermann and Parillo [57].

To exploit the group symmetry for SDPs, we will consider in the following a semidef-inite programming problem in its standard form (1.3.1) on page 13 for which an optimal solution does exist. Further, we will denote its respective dual problem in the notation including a slack variable S ∈ Sn×n

+ such that the set of dual feasible solutions can be

described by ( (y, S) m X i=1 (yiAi) + S = A0; S 0 ) .

Furthermore, we assume within this section that the data matrices, denoted by Ai with

i = 0, ..., m, belong to the commutant of a nontrivial multiplicative group of permutation matrices GSDP ⊆ Πn, say A0, ..., Am∈ CGSDP. In particular, let GSDP ⊆ Πnsuch that

AiPπ = PπAi ∀Pπ ∈ GSDP; i = 0, ..., m.

Theorem 1.3.7. ([38], Theorem 2) If the primal SDP problem (1.3.1) on page 13 and its dual problem meet the Slater condition, then there exists an optimal primal-dual pair (X, S) ∈ CGSDP × CGSDP.

(29)

RGSDP (X ∗_{) :=} 1 |G_SDP| X Pπ∈GSDP PπX∗PπT is also optimal.

The proofs of the previous two theorems can be found in [38] and in Appendix A.2 on page 137, respectively. Thus, we have shown that every SDP of the form (1.3.1) has an optimal solution in the commutant of GSDP. In general, this solution can be constructed

via RGSDP, the so-called Reynolds operator of GSDP.

To illustrate the benet gained by exploiting symmetry reduction in semidenite pro-gramming, let us consider the following small example.

Example 1.3.9. Consider the 2-dimensional semidenite programming problem min X0tr (A0X) = minX0X11+ X22 (1.3.4) with A0 = 1 0 0 1 and X ∈ S2×2

+ . This problem is invariant under the group

GSDP = 1 0 0 1 , 0 1 1 0 . Thus, CGSDP is given by:

C_G_SDP :=_{X ∈ S}2×2₊ |X₁₁= X22 .

The following linear program results after applying symmetry reduction: min X∈C_GSDP∩S2×2 + tr (A0X) = min X11≥0 2X11.

In the thesis, we will apply this technique of symmetry reduction on more complex SDPs. However, Example 1.3.9 illustrates how to reduce the complexity beforehand via symmetry reduction in order to simplify a problem. Thus, with this approach, higher-dimensional problems might be solved in practice. In [40], de Klerk, Dobre and Pasechnik propose a pre-processing technique for QAP instances that exhibit algebraic symmetry.

1.4 Overview of the thesis

Since we want to introduce a new SDP-based branch-and-bound framework throughout the thesis, we subdivide our work in the following three dierent parts:

1. The lower bound of the QAP

(30)

thesis dierent SDP relaxations for the QAP that can be classied via a dierent trade-o between eciency and strength of the resulting bound.

We start in Chapter 2 with a benchmark SDP relaxation for the QAP, the so called Zhao, Karisch, Rendl, Wolkowicz QAPR3

relaxation that represents a milestone in the eld of semidenite programming applied to the QAP due to its strength. We will modify this strong SDP relaxation in order to make it Slater feasible. In Chapter 3, we will introduce a Sherali-Adams type relaxation based on the refor-mulation linearization technique (RLT) and combine it with a semidenite program-ming relaxation. The resulting new relaxation leads to some best known bounds (at that time) that we achieved to obtain for certain maximum k-section instances on strongly regular graphs.

Moreover, in Chapter 4, we will expand our portfolio of SDP relaxations with some new smaller dimensional SDP relaxations of order O (n) based on spectral decom-position that are applicable for higher dimensional QAPs. There, we will introduce amongst others our new eigenspace relaxations called (QAPEig)and (QAPEig+SDRMS).

In particular with (QAPEig), we complement the existing SDP relaxations from

lit-erature for a smart size-dependent choice within our branch-and-bound approach. Further, via the stronger (QAPEig+SDRMS)relaxation, we will obtain improved lower

bounds for some QAPLIB instances. Therefore, we will be able to compute new best known lower bounds for the Tail35b and Tail40b instances.

2. The upper bound of the QAP

In the second part of the thesis, we search for a good approach to obtain an upper bound. Therefore, in Chapter 5, we work out a comprehensive comparison of QAP heuristics from literature and introduce a new insightful heuristic that is based on Karush-Kuhn-Tucker [82, 88] points (KKT points) of a quadratic programming reformulation of the QAP.

3. The complete branch-and-bound framework

(31)

(32)

Semidenite programming

relaxations of order n

2 Preamble

Within this chapter, we introduce the Zhao-Karisch-Rendl-Wolkowicz-type SDP relax-ations based on [128]. Further, we will develop Slater feasible versions of these type of relaxation and show how to exploit group symmetry for the QAPR3

relaxation. There-fore, it is not surprising that the results presented in the course of Chapter 2 are not new. Instead they form a collection of results spread over the literature, for example [109, 45]. Herein, the intention is to state a rst established SDP bound that is well known to oer strong bounds. Further, we include a symmetry reduced version of QAPR3

and modify the Zhao-Karisch-Rendl-Wolkowicz-type SDP relaxations to be Slater feasible and, thus, to be solvable in a numerical stable way in practice.

2.1 The QAP

R2

and the QAP

R3

relaxation

Semidenite programming has turned out to be a powerful tool in the sense of providing strong lower bounds to the QAP. In 1998, Zhao, Karisch, Rendl and Wolkowicz [128] introduced a seminal type of semidenite programming relaxations for the QAP, in par-ticular the QAPR2

and the QAPR3

relaxation. This type of relaxations is based on the idea to lift the problem to a higher-dimensional space, Sn2_×n2

, of symmetric matrices. Therefore, let us introduce the matrix

Y = vec (X) vec (X)T ∈ Sn+2×n2 (2.1.1)

for X ∈ Πn. Applying the property of the Kronecker product (see (1.1.7) on page 4), we

(33)

In the following, we rst state a Zhao-Karisch-Rendl-Wolkowicz-type SDP relaxation called (QAPR2). Therefore, we will dene the matrix

G := I ⊗ (J − I) + (J − I) ⊗ I ∈ Rn2_×n2

related to the so called Gangster operator ˆG : Rn2×n2 → Rn2_×n2

with ˆ

G (Y ) = (Jn2 − G) ◦ Y ∀Y ∈ Rn 2_×n2

that was introduced in [128] and [111]. The Gangster operator xes certain entries of Y to zero as implied by (2.1.1). (QAPR2)                            min Y tr [B ⊗ A + Diag (vec (C))] Y s.t. tr (J ⊗ J) Y = n2 Yij = 0 ∀ (i, j) : Gij = 1 tr (I ⊗ Eii) Y = 1 ∀i = 1, ..., n tr (Eii⊗ I) Y = 1 ∀i = 1, ..., n Y 0. (2.1.2)

The (QAPR2) relaxation was originally introduced in 1998 by Zhao et al. [128] and later

reformulated to the form (2.1.2) by Povh and Rendl [109].

To verify that (QAPR2) is indeed a relaxation of the quadratic assignment problem

(1.1.4), note that Y = vec (X) vec (X)T _{is a feasible point of (2.1.2) for X ∈ Π}

n. The

(QAPR2)relaxation is obtained by relaxing the condition (2.1.1) to Y 0. Further, note

that

Y(ij)= coli(X) colj(X)T (2.1.3)

with X ∈ Πn. Thus, the sum of all block matrices for the special case i = j is the identity

matrix: _n

X

i=1

Y(ii)= I. The rst constraint of (QAPR2), tr (J ⊗ J) Y = n

2_{, follows from the fact that the}

per-mutation matrix X ∈ {0, 1}n×n _{consists of exactly n entries equal to one. Furthermore,}

forcing all entries Yij to zero if Gij = 1 represents a relaxation of the orthogonality

condi-tion

XXT = XTX = I,

which also has to be satised for all permutation matrices X ∈ Πn. The nal two

con-straints, tr (I ⊗ Eii) Y = 1 and tr (Eii⊗ I) Y = 1, are relaxations of the assignment

constraints

Xe = XTe = e.

(34)

Lemma 2.1.1. ([109], Lemma 6) A matrix Y =    Y(11) · · · Y(1n) ... ... ... Y(n1) _{· · ·} _Y(nn)   ∈ S n2_×n2 + , Y(ij)∈ Rn×n, i, j = 1, ..., n,

is feasible for (QAPR2) if and only if Y satises

1. Yij = 0 for all i, j = 1, ..., n2 with Gij = 1,

2. tr Y(ii)_{= 1} _{for 1 ≤ i ≤ n;} n X i=1 diag Y(ii) = e, 3. eT_Y(ij)_{= diag Y}(jj)T for 1 ≤ i, j ≤ n, 4. n X i=1

Y(ij)= e diagY(jj)T for 1 ≤ j ≤ n.

2.1.2 The (QAPR3) relaxation

Restricting the variables in (2.1.2) to the nonnegative orthant strengthens the relaxation. The resulting relaxation is called (QAPR3):

(QAPR3)                                min Y tr [B ⊗ A + Diag (vec (C))] Y s.t. tr (J ⊗ J) Y = n2 tr (GY ) = 0 tr (I ⊗ Eii) Y = 1 ∀i = 1, ..., n tr (Eii⊗ I) Y = 1 ∀i = 1, ..., n Y 0 Y ≥ 0. (2.1.4)

For the (QAPR3) relaxation, the following lemma gives us an explicit description of the

feasible set with Y ∈ Rn2_×n2

in its block matrix structure Y(ij) _{∈ R}n×n _{accordingly to}

(2.1.3).

Lemma 2.1.2. ([41], Theorem 3) A nonnegative matrix Y =    Y(11) _{· · ·} _Y(1n) ... ... ... Y(n1) · · · Y(nn)   ∈ S n2_×n2 + , Y(ij)∈ Rn×n, i, j = 1, ..., n,

(35)

3. eT_Y(ij) _{= diag Y}(jj)T for 1 ≤ i, j ≤ n, 4. n X i=1 Y(ij) = e diag Y(jj) T for 1 ≤ j ≤ n.

The resulting (QAPR3) relaxation oers strong bounds in practice but due to its

complexity in the constraints and in the dimension of the matrix variable Y ∈ Rn2×n2

, this SDP relaxation turns out to be too time-consuming for solving larger QAP instances, say if n ≥ 14. For numerical results see Chapter 2.4 or the more detailed tables worked out by Rendl and Sotirov (Table 3 and 4 in [111]).

Further, it is not obvious that we can replace the positive semideniteness con-straint Y 0 for both, (QAPR2) and (QAPR3), by its seemingly stronger constraint

Y − vec (X) vec (X)T 0.

Lemma 2.1.3. (see e.g. [61], Proposition 7) Let M ∈ Sn×n _{such that diag (M) = cMe}

for some c ∈ R. Then the following are equivalent: M diag (M ) diag (M )T 1 0,

M is positive semidenite and eT_{M e ≥ (trM )}2_.

Indeed, in our particular case of a feasible Y for (QAPR2)or (QAPR3), respectively, the

assumptions in Lemma 2.1.3 hold for c = n since n2_{= hJ}

n2, Y i = eTY eand (trY )2 = (n)2

are satised and Y 0. Thus, using the essence of Lemma 2.1.3 and applying it to the Schur complement Theorem 1.3.6 on page 15, we get the following necessary condition for all feasible solutions for each of the Zhao, Karisch, Rendl, Wolkowicz relaxations.

Corollary 2.1.4. Every feasible solution Y of (QAPR2) and (QAPR3), respectively, also

satises:

Y − vec (diag (Y )) vec (diag (Y ))T 0. Lemma 2.1.5. (QAPR2) and (QAPR3) are not strictly feasible.

Proof. We will show that if Y is feasible for (QAPR2) and (QAPR3), respectively then

Y v = 0for v = eT, 0T, ..., 0T, −eTT ∈ Rn2

.

Therefore, we use the third property of Lemma 2.1.1 and 2.1.2, respectively, Y(ji)e = diag

Y(jj)

,

and we exploit the block-structure (2.1.3) of Y such that it follows: Y(ji)v = Y(j1)e − Y(jn)e = diag Y(jj) − diagY(jj) = 0.

(36)

2.2 Slater feasible versions

In theory, it can be shown that the Zhao-Karisch-Rendl-Wolkowicz-type SDP relaxations, (QAPR2) and (QAPR3), have a feasible point and the existence of an optimal solution,

therefore, follows from the compactness of the feasible sets. Alternatively, it follows from the conic duality theorem (see [112] or Theorem 1.1 in [42] respectively) that if a conic optimization problem is primal feasible and dual strictly feasible than there exists a primal optimal solution.

For the duals of (QAPR2) and (QAPR3), strictly feasibility can be shown. But there

exist instances where the dual feasible set is unbounded and there is no dual optimal solution. As one can see in Section 2.4, already for small instances, say n = 6, inaccuracies of 10−3 _{appear while solving (QAP}

R2) and (QAPR3); this is due to the lack of strict

feasibility.

Thus, in this section, it is our aim to analyze the problem of this numerical instability and to nd a modied reformulation of (QAPR2) and (QAPR3) which can guarantee

numerical stability when using interior-point methods.

To avoid these numerical problems we need to generate a valid relaxation for which it can be shown that there exists a strictly feasible solution (i.e. where the Slater conditions are satised) for the primal problem. To do so, let us keep the following well-known fact in mind.

Lemma 2.2.1. For any symmetric matrix M ∈ Sn×n _{it follows}

R(M ) = (Null(M ))⊥.

Since we already know from Lemma 2.2.1 that the range space is orthogonal to the null space of a symmetric matrix, we can use this knowledge to initially generate a basis of the smaller dimensioned null space. Then we can nd its orthogonal complement to nally receive a basis for the range.

Therefore, let us consider the following convex combination of outer products of all permutation matrices as a feasible solution for QAPR2

and QAPR3 , respectively: ˆ Y = 1 n! X XΠn vec (X) vec (X)T ∈ Rn2×n2. (2.2.1)

The matrix ˆY is structured in blocks of size (n × n) where the diagonal elements equals

1

n. Furthermore, in the non-diagonal blocks all elements are equal to 1

n(n−1) except their

diagonal elements that equal zero. Therefore, we have the following structure for the matrix ˆY: ˆ Y =          ₁ nIn h ₁ n(n−1)(J − I) i · · · h 1 n(n−1)(J − I) i h 1 n(n−1)(J − I) i ... ... ... ... ... ... h 1 n(n−1)(J − I) i h 1 n(n−1)(J − I) i · · · h 1 n(n−1)(J − I) i ₁ nIn          .

(37)

The-orem 3.1): ˆ Y = 1 n2J ⊗ J + 1 n2_{(n − 1)}(nI − J ) ⊗ (nI − J ) . (2.2.2)

Lemma 2.2.2. Let ˆY ∈ Rn2×n2 be dened as in (2.2.2), then its rank equals

rank( ˆY ) = (n − 1)2+ 1. (2.2.3)

Proof. We will regard in this proof the two terms of (2.2.2) separately. Thus, we dene M1, M2 ∈ Rn 2_×n2 so that M1 := 1 n2J ⊗ J, M2 := 1 n2_{(n − 1)}(nI − J ) ⊗ (nI − J ) .

Claim. Rank (M1) = 1with nonzero eigenvalue equals to 1 and eigenvector e ∈ Rn

2

. Note that J ⊗ J is the all-ones matrix and, thus, this is trivial.

Claim. M2 has eigenvalues equal to _(n−1)1 (with multiplicity (n − 1)2) and the eigenvalue

0(with multiplicity 2n − 1).

This can be proved by using the fact that the eigenvalues of the Kronecker product of two arbitrarily matrices A ∈ Rn×n _{and B ∈ R}m×m_{equals the product of each eigenvalues.}

Thus, let λi with i = 1, ..., n and µj with j = 1, ..., m denote the eigenvalues of A and B,

respectively. Then the nm eigenvalues of A ⊗ B are the following: λ1µ1, ..., λ1µm, λ2µ1, ..., λnµm.

Using this fact yields the eigenvalues of M2 , namely n

2

n2_(n−1) =

1

(n−1) with multiplicity

(n − 1)2 and 0 with multiplicity 2n − 1. To prove that e ∈ Rn2

is indeed an eigenvector of M2 with eigenvalue zero, we use the

properties (1.1.5) and (1.1.6) on page 3 of the Kronecker product so that: tr (J ⊗ J ) (nI − J ) ⊗ (nI − J ) = tr [J (nI − J )] tr [J (nI − J )] = 0. So, we can conclude that the rank of ˆY is the sum of the ranks of M1 and M2:

rank( ˆY ) = rank(M1) + rank(M2)

= 1 + (n − 1)2.

Thus, we proved that ˆY ∈ Rn2×n2 does not have full rank: rank( ˆY ) = (n − 1)2+ 1 < n2.

This means that ˆY is not positive denite. In general, ˆY denes a solution in the relative interior with the biggest rank due to its construction. We will show in the following how to nd a valid formulation of lower dimension, say (n−1)2₊₁_{, which has a strictly feasible}

(38)

Lemma 2.2.3. Let x, y ∈ Rn _{be any vectors chosen such that e}T_{x = e}T_{y = 0}_{. Then}

ˆ

Y ∈ Rn2×n2 as dened in (2.2.1) satises ˆ

Y (x ⊗ e + e ⊗ y) = 0_n2. (2.2.4)

Proof. To prove (2.2.4), we divide the formula into two parts and concentrate on each one separately.

For any X ∈ Πn it follows that

vec (X) vec (X)T [x ⊗ e] = vec (X)    col1(X) ... coln(X)    T   x1e ... xne    = vec (X) (x1+ ... + xn) = vec (X) eTx = 0. Analogously, we receive the same result for the second part:

vec (X) vec (X)T

[e ⊗ y] = 0.

So both, x ⊗ e and e ⊗ y, lie in the null space of ˆY. Thus, (x ⊗ e + e ⊗ y) is also contained in the null space.

In fact, all vectors in the null space of ˆY are of the form (x ⊗ e + e ⊗ y) for some x, y ∈ Rnwith eTx = eTy = 0.

Lemma 2.2.4. The null space of ˆY is dened by

Null ˆY= span {e ⊗ êi, êj⊗ e i, j = 1, ..., n − 1} (2.2.5) where ˆ e1, ê2, ..., ên−1 :=       1 0 0 −1       ,        0 1 0 ... −1        , ...,        0 ... 0 1 −1        .

Proof. Applying Lemma 2.2.3 implies that ˆ

L := span {e ⊗ ˆei, ˆej⊗ e i, j = 1, ..., n − 1} ⊆ Null ˆY

.

We still need to show that equality holds. Therefore, let us consider Γ := z ∈ Rn

eTz = 0 that has dimension n − 1 and, thus,

(39)

From (2.2.3) we know that rank ˆY = 1 + (n − 1)2. This implies that dim Null ˆY = n2− rank ˆY = n2−h1 + (n − 1)2i = 2n − 2.

Since dim ˆL= dim

Null ˆY with 2n−2 linear independent basis vectors of the form e ⊗ êi, êj⊗ e i, j = 1, ..., n − 1, it follows that ˆ L = span {e ⊗ êi, êj⊗ e i, j = 1, ..., n − 1} = Null ˆY . This completes the proof.

The next step is to nd a basis for the range of ˆY by constructing (n − 1)2+ 1 linear independent vectors that are orthogonal to the basis vectors of the null space (2.2.5). Lemma 2.2.5. The range of ˆY is given by

R ˆY

= span {e ⊗ e, ê1⊗ ê1, ê1⊗ ê2, ..., ê1⊗ ên, ê2⊗ ê1, ..., ên−1⊗ ên−1}

= span {e ⊗ e, ˆei⊗ ˆej ∀i, j = 1, ..., n − 1} . (2.2.6)

Proof. We check element-wise that each element of the basis of Null ˆY is orthogonal to all elements of (2.2.6). Remark that eT_e_ˆ

i = 0such that

(e ⊗ e)T (e ⊗ êi) = eT ⊗ eT (e ⊗ êi) = eTe ⊗ eTeî = 0.

Further, the fact that ˆejTe = 0leads to:

( êi⊗ êj)T ( êk⊗ e) = eîTeˆk ⊗ êjTe = 0.

Analogously, we can continue for the remaining combinations. Accordingly to the rank-nullity theorem,

n2 = dimR ˆY+ dimNull ˆY, it follows directly that

dimR ˆY= n2− (2n − 2) = (n − 1)2+ 1.

These results and the fact that all (n − 1)2_{+ 1} _{vectors in (2.2.6) are linear independent}

prove the lemma.

According to (2.2.6), we dene the basis of the range of ˆY by V := [e ⊗ e, ê1⊗ ê1, ê1⊗ ê2, ..., ên−1⊗ ên−1] ∈ Rn

2_×(n−1)2

(40)

Theorem 2.2.6. Let Y ∈ Sn2_×n2

+ be a feasible solution for (QAPR2) and (QAPR3),

re-spectively then the range of Y is contained in the column space of V ∈ Rn2_×(n−1)2₊₁

. Proof. To show R (Y ) ⊆ R (V ), we will prove that

Null (Y ) ⊇ Null (V ) . This can be subdivided according to (2.2.5) such that

Y (e ⊗ ˆek) = 0 and Y (ˆej ⊗ e) = 0

has to hold.

First, we will prove that Y (e ⊗ ˆek) = 0 is always satised. Therefore, let us consider

Y(ij) in its block-wise notation (2.1.3) such that

Y (e ⊗ ˆek) =    Y(11) · · · Y(1n) ... ... ... Y(n1) _{· · · Y}(nn)       ˆ ek ... ˆ ek    =    Y(11)_e_ˆ k+ . . . + Y(1n)eˆk ... Y(n1)eˆk+ . . . + Y(nn)eˆk    =              n X j=1 Y(1j)  eˆ_k ...   n X j=1 Y(nj)  eˆk            . Since

Y(ij)= coli(X) colj(X)T

holds for X ∈ Πn, we can use the fact that, for a permutation matrix, the columns sum

up to one such that

n

X

j=1

Y(ij)= coli(X) eT.

Since coli(X) ∈ {0, 1}n with just one entry equal to one, it follows in combination with

any vector ˆek∈ Rn that

  n X j=1 Y(ij)  eˆ_k = 0.

We still have to show that Y ( ˆej⊗ e) = 0 is satised. Therefore, we consider Y(ij) again

in its block-structure such that

(41)

As a direct consequence of the third property of Lemma 2.1.1 (respectively the third property of Lemma 2.1.2), it follows that

Y(ij)e − Y(in)e = eTY(ji)T −eTY(ni)T

= diagY(ii)T − diagY(ii)T = 0n,

for i, j = 1, ..., n.

Corollary 2.2.7. For (QAPR2) and (QAPR3), ˆY is a feasible solution of maximal rank.

Proof. This statement follows directly from Theorem 2.2.6 and from the fact that the basis of the range of ˆY is precisely V , see (2.2.3).

The following property origins from linear algebra and is crucial for the correctness of the upcoming statement in Corollary 2.2.9.

Lemma 2.2.8. Let R (Υ) ⊆ R (M) for some M ∈ Rn×k_{and Υ ∈ S}n×n

+ with rank (Υ) = k.

Then there exists an U ∈ Sk×k

++ such that

Υ = M U MT. Proof. Applying the spectral decomposition on Υ yields

Υ = q1 · · · qk    λ1 ... λk       qT₁ ... qT_k   , (2.2.8)

where k = rank (Υ), qi ∈ Rn denotes the ith eigenvector of Υ with i = 1, ..., k and

λi ∈ R++ its corresponding eigenvalue. Since R (Υ) is contained in the range of M, there

exist k vectors si∈ Rnwith i = 1, ..., k such that

M si= qi i = 1, ..., k.

We may rewrite (2.2.8) in the following way: Υ = M s1 · · · M sk    λ1 ... λk       (M s1)T ... (M sk)T    = M SΛST MT_, with S = s1 · · · sk ∈ Rn×k and Λ =    λ1 ... λk    ∈ S k×k

++. Let us dene the

matrix U := SΛST_{. Thus, U ∈ S}k×k

++ is positive denite by construction, and satises

Υ = M U MT_.

(42)

Corollary 2.2.9. For any feasible solution Y ∈ Sn2_×n2

+ of (QAPR2) or (QAPR3), one has

Y = V U VT for some U ∈ S(n−1)2+1×(n−1)2+1

+ and V ∈ Rn

2_×(n−1)2₊₁

as given in (2.2.7). Thus, we can rewrite

ˆ Y = 1 n! X X∈Πn vec (X) vec (X)T = 1 n! X X∈Πn V UXVT = 1 n!V X X∈Πn UX ! VT, where ˆU := P

X∈ΠnUX 0 is positive semidenite as the sum of positive

semide-nite matrices. Furthermore, it can be shown via spectral decomposition that ˆU is, by construction, a full rank matrix, positive semidenite and, additionally, has eigenvalues λ1 = 1 and λi = _n−11 that are strictly positive for n > 1. Thus, ˆU has to be positive

denite. Therefore, ˆU is a Slater feasible point for the upcoming formulations of the Zhao-Karisch-Rendl-Wolkowicz-type relaxations in (2.2.9) and (2.2.13) for Y = V UVT_. 2.2.1 The Slater feasible version of the QAPR2

relaxation

We have shown the existence of a strictly feasible solution of QAPR2

and QAPR3

, respectively within the previous section. Based on those theoretical results, we want to demonstrate how to transform the QAPR2

relaxation to an equivalent h(n − 1)2+ 1 i

dimensional problem with a strictly feasible solution. This idea is related to the one used by Povh and Rendl [109].

Lemma 2.2.10. One may rewrite QAPR2

as:                            min U V T_{A ⊗ BV, U}_{+ U, V}T_{Diag (vec (C)) V} s.t. VT _{(J ⊗ J ) V, U = n}2 VT_E ijV, U = 0 ∀ (i, j) : Gij = 1 VT _(E ii⊗ I) V, U = 1 ∀i = 1, ..., n VT _{(I ⊗ E} ii) V, U = 1 ∀i = 1, ..., n U 0. (2.2.9)

Proof. The proof can be constructed by applying Y = V UVT _{in (2.1.2) on each constraint}

of QAPR2

as well as on the objective. For the objective function, we obtain: tr [B ⊗ A + Diag (vec (C))] Y = VT_{A ⊗ BV + V Diag (vec (C)) V}T_{, U}_.

(43)

tr [(J ⊗ J) Y ] = n2 ⇔ VT_{(J ⊗ J ) V, U} = n2, tr [(Eii⊗ I) Y ] = 1 ⇔ VT_(E ii⊗ I) V, U = 1, ⇔ tr [(I ⊗ E_ii) Y ] = 1 ⇔ VT_{(I ⊗ E} ii) V, U = 1.

Fixing the variable equal to zero for those elements that are forced from the Gangster constraint to be zero, we obtain:

Yij = 0 ∀ (i, j) : Gij = 1

⇔ VT_E ijV, U

= 0 ∀ (i, j) : Gij = 1.

Thus, we end up with an equivalent problem in its form (2.2.9). Finally, we can simplify (2.2.9) as described below.

Lemma 2.2.11. The constraint VT _{(J ⊗ J ) V, U = n}2 _{in (2.2.9) is equivalent to U} 11= 1

n2.

Proof. The denition of V ∈ Rn2_×(n−1)2₊₁

accordingly to (2.2.7) implies that

VT    1 ... 1   =      n2 0 ... 0     

is satised. Further, it is easy to see that the Kronecker product of the all-ones matrix J ∈ Rn×n with itself can be written as follows:

(J ⊗ J ) = (e ⊗ e) (e ⊗ e)T .

Thus, it is possible to simplify the rst term of the inner product in the following way: VT(J ⊗ J ) V = VT (e ⊗ e) (e ⊗ e)T V = VT(e ⊗ e) VT (e ⊗ e)T = n2_{, 0, · · · , 0}T n2_{, 0, · · · , 0} = n4E11, with E11∈ R(n−1) 2_+1×(n−1)2₊₁

A semidefinite programming based branch-and-bound framework for the quadratic assignment problem

Tilburg University