Minimum energy configurations on a toric lattice as a quadratic assignment problem. Part I: comparison of bounds

We consider three known bounds for the quadratic assignment problem (QAP): an eigenvalue, a convex quadratic programming (CQP), and a semidefinite programming (SDP) bound. Since the last two bounds were not compared directly before, we prove that the SDP bound is stronger than the CQP bound. We then apply these to improve known bounds on a discrete energy minimization problem, reformulated as a QAP, which aims to minimize the energy between repulsive particles on a toric grid. Thus we are able to prove optimality for several configurations of particles and grid sizes, complementing earlier results by Bouman, Draisma and Van Leeuwaarden [ {\em SIAM Journal on Discrete Mathematics}, 27(3):1295--1312, 2013].


Introduction
A quadratic assignment problem in Koopmans-Beckmann form is given by three matrices A = (a ij ), B = (b ij ), C = (c ij ) ∈ R n×n , and can be written as where S n is the set of all permutations of n elements. If C = 0, then we shorten the notation to QAP (A, B). This is a quadratic optimization problem, which can be seen if we write the objective using permutation matrices: where X, Y = tr(X T Y ) is the trace inner product, and Π n the set of n × n permutation matrices.
Because of the very general form of the problem, it is not surprising that it is NP-complete (see for example §7.1.7 in [5]), which motivates the search for good approximations and bounds; see, e.g., the survey [12] and the book [5] for an overview. In Section 2 we describe three such bounds, in both increasing complexity and strength. The first is a projected eigenvalue bound, which was first introduced in [10], which, similar to the eigenvalue bound of [8], is based on the eigenvalues of the data matrices. The second bound, a convex quadratic bound, then improves this bound by adding a convex quadratic term to the objective, as introduced in [1] (see also [2,3]). The third bound, which was introduced in [20] and later reformulated in [16], is a semidefinite programming relaxation of the quadratic assignment Minimum energy configurations on a toric lattice as a quadratic assignment problem.
Part I: comparison of bounds problem. As it is the most complex computationally, it is natural to expect it to be stronger than the two other bounds, which we prove in our first main result Theorem 2.3.
In Section 3 we then apply the three bounds to a discrete energy minimization problem. It was first described in [17] as the problem of printing a particular shade of grade, by repeating the same tile of black and white squares in all directions. Other applications from physics is the search for ground states of a two-dimensional repulsive lattice gas at zero temperature ( [19]), and more generally the Falicov-Kimball model ( [7,11]), which is relevant for modelling valence fluctuations in transition metal oxides, binary alloys and high-temperature super-conductors ( [19]).
To get a distribution as equal as possible, it is natural to view this as a problem of minimizing the energy between repulsive particles on a grid. This problem can then be reformulated as a quadratic assignment problem, which allows us to apply the three bounds of Section 2 to this problem. We will see in Proposition 3.4 and Proposition 3.5 that both the projected eigenvalue bound, as well as the convex quadratic programming bound, coincide with an eigenvalue bound for this problem introduced in [4]. Finally we do calculate the bounds for instances on different grid sizes, even the semidefinite programming bound, despite its size, and prove this way optimality of certain grid arrangements. This way we expand the results by Bouman, Draisma and Van Leeuwaarden [4].
In Part II of this paper we describe the technique we used to calculate the semidefinite programming bound, which involves a symmetry reduction of the problem to a more manageable size.

Bounds for quadratic assignment problems
In this section we will consider three different bounds for QAPs, of increasing computational complexity. These are then compared to each other in Section 2.4, and applied to an energy minimization problem in Section 3.

Projected eigenvalue bound
The first bound relevant for this paper is the projected eigenvalue bound, which was introduced in [10], a stronger variant of the eigenvalue bound for QAP (see [8]), which is based on projecting the matrices into a space the same dimension as the span of the permutation matrices. [10], cf. Prop. 7.23 in [5]). Let V be the n × (n − 1) matrix, of which the columns form an orthonormal basis of the orthogonal complement of the all-ones vector e. DefineÃ = V T AV ,B = V T BV , and collect their eigenvalues in the vectors λÃ and µB respectively. Set D = 2 n Aee T B. The projection bound for the symmetric QAP(A, B) is given by where x, y − = min ϕ∈Sn n i=1 x ϕ(i) y i . One then has P B(A, B) ≤ QAP (A, B).
One may calculate PB(A, B) by sorting λÃ and µB to compute λÃ, µB − (see Proposition 5.8 in [5]) and solving one linear assignment problem min ϕ∈Sn n i=1 d iϕ(i) .

QP bound
The second bound we consider is a convex quadratic programming (CQP) bound, introduced in [1], which is based on the same projection as the bound in Proposition 2.1. We will see that it is at least as good as the projected eigenvalue bound. Here we relax X ∈ Π n to Xe = X T e = e and X ≥ 0, i.e. we optimize over doubly stochastic matrices instead of permutation matrices.
Hadley, Rendl and Wolkowicz ( [9], [10]) observed that every doubly stochastic matrix can be written as where V is the n × (n − 1) matrix of which the columns form an orthonormal basis of the orthogonal complement of e, as before. We have V T V = I n−1 , V V T = I n − 1 n ee T and Y = V T XV . As before, we setÃ = V T AV and B = V T BV , and collect their eigenvalues in the vectors λÃ and µB.
In Section 3 of [1], Anstreicher and Brixius introduce the following CQP bound for quadratic assignment problems.
Minimum energy configurations on a toric lattice as a quadratic assignment problem.
Part I: comparison of bounds In other words, one always has P B(A, B) ≤ QP B(A, B).
One may compute QPB(A, B) by solving a linear assignment problem to obtainQ, and then solving a CQP in O(n 2 ) variables. For details, see Section 4 in [1].

SDP bound
The following semidefinite relaxation for QAP(A, B, C) was studied by Povh and Rendl [16], which is equivalent to an earlier bound by Zhao, Karisch, Rendl and Wolkowicz [20]: where A, B, C ∈ R n×n , A and B are symmetric, I n and J n are the identity and all-ones matrices respectively of size n × n, E ij denotes the n × n matrix with a single one at position (i, j), and ⊗ the Kronecker product. We write The bound SDP QAP (A, B, C) is expensive to compute, as it involves an SDP with doubly nonnegative matrices of order n 2 × n 2 .

Bound comparison
The three bounds P B(A, B), QP B(A, B) and SDP QAP (A, B) increase in computational complexity, hence we would expect that the bounds do get better accordingly. As it turns out, we can show the expected order of the bound quality. Proof. The only inequality we have to show is QP B(A, B) ≤ SDP QAP (A, B), since the leftmost inequality was shown in [1]. To do this, we introduce another SDP bound, which lies between QPB(A,B) and SDPQAP(A,B). Again we start with the observation that for every doubly stochastic X we can always find an Y with where V is the n × (n − 1) matrix of which the columns form an orthonormal basis of the orthogonal complement of the all one vector e. Set y = vec(Y ). The idea of this bound is now to relax U = yy T to be a semidefinite matrix with certain constraints. We also add a vector variable u with U − uu T 0, which relaxes y = vec(Y ) itself.
Minimum energy configurations on a toric lattice as a quadratic assignment problem. Part I: comparison of bounds We can rewrite the objective function of the QAP at X = 1 n ee T + V Y V T as: Thus we can write it as linear function of yy T and y, which we relax to U and u respectively. To make this bound at least as good as the convex quadratic bound QP B(A, B), we have to add more conditions, which follow from and analogously we can show that E ij ⊗ I, U = δ ij as well. Finally, the property that X = 1 Thus we get a semidefinite programming relaxation of QAP by We now want to show that SDP P B(A, B) ≥ QP B(A, B). For this let (U, u) be an optimal solution of SDP P B(A, B). Then we construct a feasible solution for QP B(A, B) by setting y = u, since (V ⊗V )u+ 1 n ee T ≥ 0, and thus X = 1 By the Schur complement theorem we know that U − uu T 0, hence we have that Minimum energy configurations on a toric lattice as a quadratic assignment problem. Part I: comparison of bounds Thus we can compare the two bounds by Finally we want to show that SDP P B(A, B) is no better than the main SDP -bound SDP QAP (A, B) (2.1). For this we split the n × n-matrix-variable Y of (2.1) into n 2 blocks of size n × n, which we call Y (ij) . We write 1≤i,j≤n , and use similar notation for other block-matrices. By [16], With these we can show that we get a feasible solution for which is the transformation to a Slater-feasible variant of SDP QAP (see e.g. the thesis of Uwe Truetsch [18]). Similarly to Y , we can split U into (n − 1) 2 blocks of size (n − 1) × (n − 1), which we call U (ij) . We get an explicit formula for these blocks in terms of the Y (ij) , if we see V ⊗ V as n(n − 1) blocks of size n × (n − 1), since then all block sizes are compatible with multiplication.
We can now use (i)-(iv) to derive some properties of U . First note that by (i) and (ii) we know that tr(Y (ij) ) = δ ij , and by (ii) and (iii) that tr(Y (ij) J) = 1. Hence we see that Minimum energy configurations on a toric lattice as a quadratic assignment problem. Part I: comparison of bounds Similarly we can use To construct a feasible u with the objective value we need, we use that we can add Y − yy T 0, y = diag(Y ) to (2.1) without changing the optimal value of SDP QAP . With an optimal solution (Y, y) we thus set Since Y and Y − yy T are positive semidefinite, the matrices are positive semidefinite as well, and thus feasible for SDP P B (A, B). What remains to be seen is that the objective values of the two programs are the same.
, 3 Energy minimization on a toric grid as QAP We generalize a problem described by Taillard in [17], which models the problem of printing a certain shade of grey with only black and white squares ("pixels"). An example of these problems is included in the QAPLIB dataset [6], namely Tai64c.
The goal is to print a particular shade of grey with a given density m/n (ratio of black to total squares), which is done by repeating a grid of n = n 1 ×n 2 square cases, exactly m of which are black. We want the cases to be as regular as possible, so it is natural to model this as an energy minimization problem, with repulsive particles corresponding to the black squares. If we have two cases at locations (x 1 , y 1 ) and (x 2 , y 2 ) ∈ [n 1 ] × [n 2 ], the force between them is inverse to their distance squared, where the distance is given by the shortest path metric on the toric grid. We use this grid, since we wish to tile the plane with the n 1 × n 2 -rectangles, as can be seen in Figure 1.
if the coordinates are different, where d Lee ((x 1 , y 1 ), (x 2 , y 2 )) = min(|x 1 − x 2 |, n 1 − |x 1 − x 2 |) + min(|y 1 − y 2 |, n 2 − |y 1 − y 2 |) is the Lee distance, given by the shortest path metric on the toric grid. We set f i,i = 0. We can then formulate this problem as QAP with matrices A, B ∈ R n×n , indexed by grid points i = (x i , y i ), j = (x j , y j ) ∈ [n 1 ] × [n 2 ], Minimum energy configurations on a toric lattice as a quadratic assignment problem.
Note that the QAP has dimension n = n 1 × n 2 , and its semidefinite relaxation has dimension n 2 , which is already 4096 on an 8 × 8 grid.
To reduce the amount of cases one has to look at, we can show a that selecting the complement of an optimal solution leads to another optimal solution. Proof. We can rewrite the objective function as a,b∈T

Eigenvalue bound of Bouman, Draisma and Leeuwaarden
Problem (3.1) was considered before in [4], specifically for the case of the Lee-metric (shortest path length on the toric grid), and m = n 2 (which we can see as special case of our variant). They took a look at a different relaxation of the problem, which they call fractional total energy where B i,j is the potential energy between coordinates i and j, as defined before. If x is the characteristic vector of a subset of m coordinates, then it is exactly the potential energy of the set of particles. They identified the optimal solutions of the relaxation in terms of the eigenvalues of B. Proposition 3.2 (Proposition 2.5. in [4]). Let λ min be the smallest eigenvalue of B. Then the set of optimal solutions of the fractional energy minimization problem (3.2) consists of all vectors of the form m n e + y, where y belongs to the eigenspace of B with eigenvalue λ min , is perpendicular to e, and satisfies y T y = m − m 2 n . This proposition can be used to find lower bounds for the potential energy of a subset of coordinates, if one knows the eigenvalues of B. Corollary 3.3. Let λ min be the smallest eigenvalue of B, and λ 1 the eigenvalue of B corresponding to e. Then the following is a lower bound for the potential energy of m particles on a grid with n nodes: Proof. Since B is a circulant block-matrix of circulant matrices, e is an eigenvector of it. Let λ 1 be its corresponding eigenvalue. Furthermore, let x = m n e + y be an optimal solution of the fractional total energy as in Proposition 3.2. Then its objective value is

Bound comparison
We now want to compare the eigenvalue bound with the QAP relaxations described last section, as well as the different QAP relaxations for this specific case. The matrix A has rank one, soÃ has rank one as well. Since e is an eigenvector of B,B has the same eigenvalues as B, except for λ 1 , thus we get λÃ, µB − = λ max (Ã)λ min (B) = tr(Ã)λ min .
The eigenvalue ofÃ is exactly Combining these, we see that the projection bound is the same as the eigenvalue bound: So the eigenvalue bound is the same as the first of the QAP bounds, namely the bound P B (A, B). Furthermore, even the convex quadratic bound cannot give us better bounds here, as we show now.
Minimum energy configurations on a toric lattice as a quadratic assignment problem.
Part I: comparison of bounds Proposition 3.5. If A and B are of the forms as defined for the energy minimization problem, then we have that where D = (d ij ) = 2 n AJB andQ is positive semidefinite, as defined before (see (2.1), (2.2)).
Proof. If we eliminate the terms which appear in both programs, we see that we want to show By definition of D the two linear terms are equal, except that on the left we minimize over permutations, and on the right over doubly stochastic matrices. Because the terms are linear, and doubly stochastic matrices are convex combinations of permutations, the minima of the two linear terms are equal. SinceQ is positive semidefinite, we thus want to find a doubly stochastic X = 1 n J + V Y V T with y TQ y = 0, which minimizes the linear term. Earlier in (3.4) we have seen that Thus it makes sense to take a look at the SDP-bound for this kind of problem, if one wants to find stronger bounds.

Numerical comparison of the bounds
Here we compare the eigenvalue-bound with the SDP-bound for the energy minimization problems. The upper bounds were found using simulated annealing. Calculating the SDP-bound directly is prohibitively slow, which is why we symmetry reduced the problems first, as described in the second part of this paper. In [4] Bouman, Draisma and van Leeuwaarden prove optimality for the checkerboard arrangement in the cases that n 1 = n 2 are even, and m = n1n2 2 . This can be seen for the grid sizes we checked as well, but we do get some more proofs of optimality.
If one of the bounds is sharp, then we get a proof of optimality for these parameters. Furthermore, we can even prove optimality in some cases even if the bound is not completely sharp, as explained in the following result.
Let S be a subset of grid points with objective value T , and let p be the SDP-bound for the problem. If T − p < r, then S is an optimal solution.
Minimum energy configurations on a toric lattice as a quadratic assignment problem.
Part I: comparison of bounds Proof. Both the matrix A and the optimization variable X ∈ Π n of this QAP are symmetric 0/1-matrices, and B is symmetric as well, with zeros on the main diagonal, which means that the objective value is of the form Hence different objective values have to at least differ by r.
Note that we could find several new optimal configurations by computing the SDP bound SDP QAP (A, B), for example the case m = 20 on a 10 × 10 grid. We do not include results for the weaker SDP bound SDP P B(A, B) in the tables, since these turned out to equal the projected eigenvalue bound for small instances. We do not know if these bounds coincide in general, though.  Minimum energy configurations on a toric lattice as a quadratic assignment problem. Part I: comparison of bounds

Conclusion
We have compared three different bounds for the quadratic assignment problem in Theorem 2.3, showing that the bounds get better with their complexity. We have then applied and compared the bounds on an energy minimization problem, proving new cases of optimality by computing the SDP bound from [20].
In the second part of this paper we go into detail on the symmetry reduction of the SDP bound [20] for the energy minimization problem. It is based on Permenter and Parrilo's Jordan reduction method (see [13,14]), which we extended to include optimization over the (non-symmetric) doubly nonnegative cone, and allows for polynomial time numeric symmetry reduction in the size of the original matrices.