Solving Challenging Large Scale QAPs

For a positive integer $n$, we let $N=\{1,\ldots,n\}$ represent a set of locations and also a set of facilities. Given $n\times n$ symmetric matrices $\mbox{\boldmath$A$}=[a_{ik}]$, $\mbox{\boldmath$B$}=[b_{j\ell}]$ and an $n\times n$ matrix $\mbox{\boldmath$C$}=[c_{ij}]$, the quadratic assignment problem (QAP) is stated as

where $a_{ik}$ denotes the flow between facilities $i$ and $k$, $b_{j\ell}$ is the distance between locations $j$ and $\ell$, $c_{ij}$ the fixed cost of assigning facility $i$ to location $j$, and $(\pi(1),\ldots,\pi(n))$ a permutation of $1,\ldots,n$ such that $\pi(i)=j$ if facility $i$ is assigned to location $j$.

The QAP is known as NP-hard in theory, and solving QAP instances of size larger than 35 in practice still remains challenging. Various heuristic methods for the QAP such as tabu search [21, 22], genetic method [5] and simulated annealing [4] have been developed. Those methods frequently attain near-optimal solutions, which often also happen to be the exact optimal solutions. The exactness is, however, not guaranteed in general.

Most existing methods for finding the exact solutions of QAPs are designed with the branch-and-bound (B&B) method [1, 3, 6, 8, 14, 16]. As its name indicates, branching and (lower) bounding are the main procedures of the method. The bounding based on doubly nonnegative (DNN) relaxations has recently attracted a great deal of attention as it provides tight bounds. Among the methods for solving large scale DNN relaxation problems, we mention four types of methods that can be applied to such DNN problems.

First two of the methods are SDPNAL+ [YST2015] (a majorized semismooth Newton-CG augmented Lagrangian method for semidefinite programming with nonnegative constraints) and BBCPOP [9] (a bisection and projection method for Binary, Box and Complementary constrained Polynomial Optimization Problems). Some numerical results on these two methods applied to QAP instances with dimensions $n=15$ to $50$ from QAPLIB [7] were reported in [9]. where BBCPOP attained tighter lower bounds for many of instances with dimensions $30$ to $50$ in less execution time. In addition, new and improved lower bounds were computed by Mittelmann using BBCPOP for the unknown minimum values of larger scale QAP instances including tai100a and tai100b. See QAPLIB [7]. The third method is an alternating direction method of multipliers (ADMM) proposed by Oliveira et al. [13] in combination with the facial reduction for robustness. The fourth method is the Newton-bracketing (NB) method [2, 10], which was developed to further improve the lower bounds obtained from BBCPOP by incorporating the Newton method into BBCPOP.

Numerical results on the NB method and BBCPOP are given in Table 2 of [10], and ones on ADMM in Table 7 of [13]. Lower bounds for the QAP instances sko42, sko49, sko56, sko64, sko72, sko81, tai60a, and tai80a were commonly computed by these methods. While the lower bounds obtained by the NB method for the first three instances are not smaller than those by ADMM and the difference is as big as $1$, the NB method clearly attained the tightest lower bound among the three methods for the last five larger-scale instances with dimension $n\geq 60$.

RLT [17] is known as a powerful method for global optimization of polynomial programming problems. In their paper [6], Goncalves et al. implemented the level 2 RLT in their B&B method run on heterogeneous CPUs and GPUs environment. They solved tai35b and tai40b from QAPLIB for the first time.

The main motivation of our project is to challenge larger scale QAP instances from QAPLIB [7] that have not been solved yet. We implement the NB method [10] combined with the B&B method in the specialized Ubiquity Generator (UG) framework [ShinanoAchterbergBertholdHeinzKoch2012] to find the exact solutions of large scale QAP instances.

The NB method has the following advantages:

(a) The NB method generates a sequence of intervals $[\mbox{lb}^{k},\mbox{ub}^{k}]$ $(k=1,2,\ldots)$ such that each $[\mbox{lb}^{k},\mbox{ub}^{k}]$ contains the lower bound $\underline{\zeta}$ to be computed and both $\mbox{lb}^{k}$ and $\mbox{ub}^{k}$ converge monotonically to $\underline{\zeta}$ as $k\rightarrow\infty$. Hence, if $\mbox{lb}^{k}$ becomes larger than the incumbent objective value $\hat{\zeta}$ at some iteration $k$, then $\hat{\zeta}<\mbox{lb}^{k}\leq\underline{\zeta}$ follows. Consequently, the node can be pruned. If $\mbox{ub}^{k}<\hat{\zeta}$ occurs at some iteration $k$, then branching the node should be determined as $\underline{\zeta}\leq\mbox{ub}^{k}<\hat{\zeta}$ holds. In these two cases, the NB method can terminate at iteration $k$ before the interval attains an accurate lower bound. This feature of the NB method works effectively to reduce a considerable amount of the computational time.

(b) The lower bounds computed by the NB method (and also the lower bounds computed by BBCPOP) are reliable or valid in the sense that they are guaranteed by very accurate dual feasible solutions of the Lagrangian DNN relaxation problem (see [10, Section 2.3].) It would be plausible to incorporate the technique used there for ADMM’s lower bound, but the resulting lower bound would deteriorate.

The UG is a generic framework to parallelize branch-and-bound based solvers and has achieved large-scale MPI parallelism with 80,000 cores [18]. It has also been generalized to handle non-branch-and-bound based solvers. Its experimental implementation for solving Shortest Vector Problem (SVP) handles more than 100,000 cores [23] .

The aim of this article is to present fundamentals of our joint project for solving larger scale QAPs, report the progress, and discuss future plans to further develop the methods. We mention that tai35b, tai40b, tai30a and sko42 from QAPLIB [7] have been successfully solved by our method, where the last two instances could be solved for the first time. QAPs with dimension at least 50 are our next target problems.

In Section 2, we first reformulate QAP (1) to QAP (6), which is a binary quadratic optimization problem form in $(1+n^{2})$-dimensional vector variable $u$. Then we introduce a class of its sub QAPs and describe the enumeration tree of those subproblems used in the branching process. In Section 3, a simple Lagrangian DNN relaxation problem of a sub QAP and the NB method to solve the relaxation problem are presented for the lower bounding procedure. In Section 4, a method to retrieve a feasible solution of (6) from an approximate optimal solution of the Lagrangian DNN relaxation problem is described for the upper bounding procedure. We present three types branching rules in Section 5, and preliminary numerical results in Section 6. Some additional techniques are proposed for future development in Section 7.

Throughout this section, we fix $(F,L)\in\Pi_{r}\times\Pi_{r}$ $(r\in\{0,1,\ldots,n\})$, and present a method to compute a lower bound for the optimal objective value $\zeta(F,L)$ of QAP$(F,L)$. Before deriving a doubly nonnegative (DNN) relaxation of QAP$(F,L)$ in Section 3.1 and a Lagrangian DNN relaxation of QAP$(F,L)$ in Section 3.2, we first transform QAP$(F,L)$ into $\mbox{\mbox{$\overline{\mbox{QAP}}$}}(F,L)$ described below.

We note that each homogeneous linear equality constraint $\displaystyle\sum_{k\in L^{c}}x_{ik}-x_{0}=0$ of QAP$(F,L)$ can be written as $\hat{\mbox{\boldmath$c$}}(i,L)^{T}\mbox{\boldmath$x$}=0$ for some $\hat{\mbox{\boldmath$c$}}(i,L)\in\mathbb{R}^{\{0,F^{c}\times L^{c}\}}$ $(i\in F^{c})$. Let $\hat{\mbox{\boldmath$a$}}(i,L)=\mbox{\bf vec}(\hat{\mbox{\boldmath$c$}}(i,L)\hat{\mbox{\boldmath$c$}}(i,L)^{T})^{T}$. Then, the linear equality is equivalent to $\hat{\mbox{\boldmath$a$}}(i,L)\mbox{\bf vec}(\mbox{\boldmath$x$}\mbox{\boldmath$x$}^{T})=0$. Similarly, we can rewrite each homogeneous linear constraint $\displaystyle\sum_{k\in F^{c}}x_{kj}-x_{0}=0$ as $\tilde{\mbox{\boldmath$a$}}(j,F)\mbox{\bf vec}(\mbox{\boldmath$x$}\mbox{\boldmath$x$}^{T})=0$ for some $(1+|F^{c}||L^{c}|)^{2}$-dimensional row vector $\tilde{\mbox{\boldmath$a$}}(j,F)$ $(j\in L^{c})$. Hence, the entire linear constraints

can be expressed as $\mbox{\boldmath$A$}(F,L)\mbox{\bf vec}(\mbox{\boldmath$x$}\mbox{\boldmath$x$}^{T})=\mbox{\bf 0}\in\mathbb{R}^{|F^{c}|+|L^{c}|}$, where $\mbox{\boldmath$A$}(F,L)$ denotes the $(|F^{c}|+|L^{c}|)\times(1+|F^{c}||L^{c}|)^{2}$ matrix consisting of the row vectors $\hat{\mbox{\boldmath$a$}}(i,L)$ $(i\in F^{c})$ and $\tilde{\mbox{\boldmath$a$}}(j,F)$ $(j\in L^{c})$.

Define

By construction, the constraint $x_{ij}(x_{0}-x_{ij})=0$ $((i,j)\in F^{c}\times L^{c}))$ is equivalent to $X_{0\mbox{\scriptsize\boldmath$\alpha$}}-X_{\mbox{\scriptsize\boldmath$\alpha$}\mbox{\scriptsize\boldmath$\alpha$}}=0\ (\mbox{\boldmath$\alpha$}\in F^{c}\times L^{c})$ if $\mbox{\boldmath$X$}=\mbox{\boldmath$x$}\mbox{\boldmath$x$}^{T}\in\mathbb{S}^{\{0,F^{c}\times L^{c}\}}$ and $\mbox{\boldmath$x$}\in\mathbb{R}^{\{0,F^{c}\times L^{c}\}}$, and $\langle\mbox{\boldmath$H$}(F,L),\,\mbox{\boldmath$X$}\rangle=X_{00}$ for every $\mbox{\boldmath$X$}\in\mathbb{S}^{\{0,F\times L\}}$. Therefore, we can rewrite QAP$(F,L)$ given in Section 2.2 as

Let $(F,L)\in\Pi_{r}\times\Pi_{r}$ for some $r\in\{0,1,\ldots,n\}$. An approximate optimal solution of QAP$(F,L)$ needs to be computed for the upper bounding procedure. A simple rounding of an approximate optimal solution of its Lagrangian relaxation DNN${}^{p}_{\lambda}(F,L)$ and the tabu search method [21, 22] for QAP$(F,L)$ have been used for the computation. We plan to replace this method by the one presented below in this section to effectively utilize information from DNN${}^{p}_{\lambda}(F,L)$. In the remainder of this section, we assume $F=L=\emptyset$ so that QAP$(F,L)$ coincides with the root node QAP (6), for simplicity of notation. All the discussions, however, can be easily adapted to any sub QAP, QAP$(F,L)$.

Suppose that we have applied the NB method to the primal-dual pair of DNN problems DNN${}^{p}_{\lambda}(\emptyset,\emptyset)$ and DNN${}^{d}_{\lambda}(\emptyset,\emptyset)$ and obtained an approximate optimal solution $\widehat{\mbox{\boldmath$X$}}$ of DNN${}^{p}_{\lambda}(\emptyset,\emptyset)$. The set of feasible solutions $X$ of DNN${}^{p}_{\lambda}(\emptyset,\emptyset)$ may be regarded as a doubly nonnegative relaxation of the set of $\mbox{\boldmath$u$}\mbox{\boldmath$u$}^{T}$ with feasible solutions $u$ of QAP (6), and the first column of $\widehat{\mbox{\boldmath$X$}}$ corresponds to a feasible (or approximate optimal) solution $u_{0}\mbox{\boldmath$u$}=\mbox{\boldmath$u$}$ of QAP (6). We know that $\mbox{\boldmath$U$}=[u_{ik}]$ forms a permutation matrix for every feasible solution $u$ of QAP (6). Through preliminary numerical experiment, we have observed that the $n\times n$ matrix

approximately forms a doubly-stochastic matrix. A popular approach for recovering a feasible solution from $\widehat{\mbox{\boldmath$U$}}$ is through appropriate rounding procedures. Here we propose a method to compute a permutation matrix $\mbox{\boldmath$U$}=\overline{\mbox{\boldmath$U$}}$ which minimize the distance $\parallel\mbox{\boldmath$U$}-\widehat{\mbox{\boldmath$U$}}\parallel$ among all permutation matrices $U$:

where $\parallel\cdot\parallel$ denotes the Frobenius norm. This problem forms a linear assignment problem, which can be solved the well-known Hungarian method [12].

Let $(F,L)\in\Pi_{r}\times\Pi_{r}$ for some $r\in\{0,1,\ldots,n\}$. As we have briefly mentioned in Section 2.3, if node QAP$(F,L)$ is not pruned, we select either an $f_{+}\in F^{c}$ to branch QAP$(F,L)$ into $|L^{c}|$ child nodes QAP$((F,f_{+}),(L,\ell))$ $(\ell\in L^{c})$ or an $\ell_{+}\in L^{c}$ to branch QAP$(F,L)$ into $|F^{c}|$ child nodes QAP$((F,f),(L,\ell_{+}))$ $(f\in F^{c})$.

We present three types of rules for selecting an $f_{+}\in F^{c}$ or an $\ell_{+}\in L^{c}$. In each rule, an evaluation function $\varphi(f,\ell)$ is introduced for possible child nodes QAP$((F,f),(L,\ell))$ $((f,\ell)\in F^{c}\times L^{c})$ so that $\varphi(f,\ell)$ can represent the optimal values of QAP$((F,f),(L,\ell))$ or its Lagrangian DNN relaxation problem DNN${}^{p}_{\lambda}((F,f),(L,\ell))$. Then, the following functions are defined:

Our method is to choose $f_{+}=\mbox{arg max}\{\bar{\varphi}(f,\cdot):f\in F^{c}\}$ if $\mbox{max}\{\bar{\varphi}(f,\cdot):f\in F^{c}\}\geq\mbox{max}\{\bar{\varphi}(\cdot,\ell):\ell\in L^{c}\}$, and $\ell_{+}=\mbox{arg max}\{\bar{\varphi}(\cdot,\ell):\ell\in L^{c}\}$ otherwise. Each rule employs a different function $\varphi(f,\ell)$ from others. We describe how $\varphi(f,\ell)$ is defined with a fixed $((f,\ell)\in F^{c}\times L^{c})$ for the three rules in Sections 5.1, 5.2 and 5.2, respectively.