How to start a heuristic? Utilizing lower bounds for solving the quadratic assignment problem

* Corresponding author E-mail: [email protected] (J. Kudela) 2022 Growing Science Ltd. doi: 10.5267/j.ijiec.2021.12.003 International Journal of Industrial Engineering Computations 13 (2022) 151–164 Contents lists available at GrowingScience International Journal of Industrial Engineering Computations home p a g e: www.Growin g Science.co m /i j iec How to start a heuristic? Utilizing lower bounds for solving the quadratic assignment problem Radomil Matouseka, Ladislav Dobrovskya and Jakub Kudelaa* aInstitute of Automation and Computer Science, Faculty of Mechanical Engineering, Brno University of Technology, Brno, Czech Republic C H R O N I C L E A B S T R A C T Article history: Received July 30 2021 Received in Revised Format September 21 2021 Accepted December 31 2021 Available online Decem b e r , 31 2021 The Quadratic Assignment Problem (QAP) is one of the classical combinatorial optimization problems and is known for its diverse applications. The QAP is an NP-hard optimization problem which attracts the use of heuristic or metaheuristic algorithms that can find quality solutions in an acceptable computation time. On the other hand, there is quite a broad spectrum of mathematical programming techniques that were developed for finding the lower bounds for the QAP. This paper presents a fusion of the two approaches whereby the solutions from the computations of the lower bounds are used as the starting points for a metaheuristic, called HC12, which is implemented on a GPU CUDA platform. We perform extensive computational experiments that demonstrate that the use of these lower bounding techniques for the construction of the starting points has a significant impact on the quality of the resulting solutions. © 2022 by the authors; licensee Growing Science, Canada Keywords: Heuristics Lower bounds Metaheuristics Quadratic assignment problem Starting values 1. Introduction The NP-hard quadratic assignment problem (QAP) in its Koopmans and Beckmann form (Koopmans & Beckmann, 1957), which is notoriously difficult in practice, can be described as follows (Cela, 2013): The problem is structured on a complete directed graph 𝐺 = (𝑉,𝐴) with 𝑛 nodes and 𝑛 arcs, together with a set of 𝑛 facilities that have to be assigned to the nodes. The indices 𝑖,𝑗 ∈ 𝑉 correspond to the nodes, the indices 𝑓,𝑔 ∈ 𝑁 = {1,...,𝑛} correspond to the facilities, 𝑏,≥0 is a given (directed) distance from node 𝑖 to node 𝑗, 𝑎,≥0 is a given flow from facility 𝑓 to facility 𝑔. By using binary variables 𝑥,=1 if facility 𝑓 is assigned to node 𝑖, and 0 otherwise, the QAP can be stated as the following quadratic 0-1 optimization problem: min 𝑎,𝑏,𝑥 𝑥,  ∈  ∈  ∈∈ (1) s.t. 𝑥,  ∈  =1 ∀ 𝑓 ∈𝑁 (2) 𝑥,  ∈ =1 ∀𝑖∈𝑉 (3) 𝑥 ∈{0 , 1} ∀𝑖∈𝑉 ∀ 𝑓 ∈𝑁, (4) It is quite well known that the constraint matrix, defined by (2)-(3) is totally unimodular, implying that the optimization of any linear objective function over the QAP feasible set is just a relatively easy linear programming problem, known as the linear assignment problem (LAP) (Burkard et al., 2012). 152 Despite its rather simple definition, QAP is among the most difficult optimization problems that arise in practice – campus planning problem (Dickey & Hopkins, 1972), backboard wiring problem (Steinberg, 1961), hospital layout problem (Helber et al., 2016), airport gate assignment (Haghani & Chen, 1998), turbine runner in electricity generation (Laporte & Mercure, 1988), statistical analysis (Hubert & Schultz, 1976), and optimal placing of letters on touchscreen devices (Dell’Amico et al., 2009) have all been modeled as a QAP. There are several other well-known combinatorial optimization problems which can be formulated as QAPs with specific coefficient matrices (Cela, 2013) – e.g., the travelling salesman problem, graph partitioning and maximum clique, the linear arrangement problem, and packing problems in graphs, to name a few. An intriguing feature of the QAP is that even for some problems of size 𝑛 ≤ 50, such as sko42 or tai30a from the QAPLIB problem library (Burkard et al., 1997), the optimal solution is still not confirmed. Even finding an ε-optimal solution is a difficult problem. There are, however, several QAP instances/structures for which the optimal solution is attainable in polynomial time (Cela et al., 2018) or which were generated in such a way that the optimal solution is known (Li & Pardalos, 1992). Furthermore, several directions for enriching the QAP formulation have been proposed – among the most notable of these are the multi-objective formulation (Samanta et al., 2018; Sanhueza et al., 2017) and stochastic programming formulation (Popela et al., 2016; Matousek et al., 2017). For these reasons, the study of the QAP attracted quite a large amount of research from both mathematical programming and heuristics communities. In this paper, we show that the approaches from these two communities can be successfully combined. We will utilize the lower-bounding techniques for the construction of advantageous starting points for a hill climbing metaheuristic and show on extensive computational experiments that starting the heuristic at these points yields significant improvements over the usual random starting points. The remained of the article is structured as follows: In Section 2, the state-of-the-art in both mathematical programming and heuristics approaches to the QAP is reviewed. In Section 3 a suitable metaheuristic algorithm HC12 is described. Section 4 provides computational comparison of the mathematical programming approaches, and the results obtain from HC12. Conclusions and future research direction are summarized in Section 5. 2. Methods for approaching QAP 2.1. (Meta)Heuristics Because of the computational difficulty of the QAP, myriads of heuristics were proposed to tackle this combinatorial problem. Among the first ones were simulated annealing (Burkard & Rendl, 1984), robust tabu search (Taillard, 1991) and genetic hybrids (Fleurent & Ferland, 1994) – although these are no longer the most efficient methods, they were able to find the best known solutions for some of the QAPLIB instances, that are yet to be beaten or proven optimal. The comparison between tabu search and simulated annealing based on a size of the QAP was conducted in (Hussin & Stützle, 2014). A local search heuristic called breakout local search enhanced by a Levenshtein Distance metric for checking solutions for similarity was described in (Aksan et al., 2017). The state-of-the-art in metaheuristics for the QAP includes population based memetic algorithms (Benlic & Hao, 2015), genetic algorithms (Ahmed, 2015), differential evolution (Hameed et al., 2020) and particle swarm algorithms (Hafiz & Abdennour, 2016). Hybrid algorithms, combining several heuristics and metaheuristics are also very prevalent. A hybrid teaching-learning based algorithm integrating tabu search within a swarm intelligence metaheuristic was described in (Dokeroglu, 2015). A memetic algorithm that uses a ternary tree structure for its population and the tabu search algorithm, which runs simultaneously, for its local search mechanism was proposed in (Harris et al., 2015). A parallel hybrid algorithm with three phases was proposed by (Tosun, 2015) – this algorithm initially benefits from a genetic algorithm to obtain a high-quality initial seed on which a diversification mechanism is run. Finally, this modified solution is used for a robust tabu search to find a near-optimal result. In (Abdel-Basset et al., 2018) the authors describe an algorithm integrating the whale optimization algorithm with a tabu search. A multistart hyper-heuristic algorithm on the grid is proposed in (Dokeroglu & Cosar, 2016) – it makes use of different metaheuristics (simulated annealing, robust tabu search, ant colony optimization, and breakout local search) and reports computations on a high-performance cluster with 368 cores and 736 GB of RAM. Since it offers speed-up opportunity that can outperform current multicore processors, (Tsutsui & Fujimoto, 2009) applied Graphics Processing Unit (GPU) computation with compute unified device architecture (CUDA) to solve the QAP. In (Czapiński, 2013) is proposed a Parallel Multistart Tabu Search (PMTS) algorithm. It is implemented on a highly powerful GPU hardware intended for high-performance computing with the CUDA platform. Therefore, PMTS is shown to perform competitively with a single-core or a parallel CPU implementation on a high-end six-core CPU. Another GPU based algorithm is described in (Mohammadi et al., 2015) – a parallel genetic algorithm, that (as the authors report) can run up to 30 times faster than its serial counterpart. Finally, the bees algorithm implemented on the CUDA platform is proposed in (Chmiel & Szwed, 2016). 2.2. MIP Reformulations One common mathematical programming approach for solving the QAP is to “linearize” it, that is, reformulate it as a pure or mixed integer linear programming problem. This was first done in (Gilmore, 1962) by replacing the terms 𝑥,𝑥, in the R. Matousek et al. / International Journal of Industrial Engineering Computations 13 (2022) 153 objective function by 𝑛 variables 𝑦,,,=𝑥,𝑥,. This reformulation was further improved upon in (Adams & Johnson, 1994), calling it the level-1 reformulation-linearization technique (RLT-1). The reformulated problem than has the following form: min 𝑎,𝑏,𝑦,,, ∈∈∈∈ (5) s.t. 𝑦,,,  ∈  =𝑥, ∀ 𝑗 ∈𝑉 ∀ 𝑓 ,𝑔∈𝑁 (6) 𝑦,,,  ∈ =𝑥  ∀𝑖 𝑗 ∈𝑉 ∀𝑔∈𝑁 (7) 𝑦,,,=𝑦   ≥0 ∀𝑖 𝑗 ∈𝑉 ∀ 𝑓 𝑔∈𝑁 (8) 𝑥,  ∈  =1 ∀ 𝑓 ∈𝑁 (9) 𝑥,  ∈ =1 ∀𝑖∈𝑉 (10) 𝑥 ∈{0,1} ∀𝑖∈𝑉 ∀ 𝑓 ∈𝑁, (11) By relaxing the binary constraint (11) (using a LP relaxation), the above formulation can be used to obtain a valid lower bound on the QAP (1)-(4). The RLT-1 reformulation was further strengthened by introducing additional 𝑛 variables, called RLT-2 in (Adams et al., 2007), and even further with additional n8 variables, called RLT-3 in (Hahn et al., 2012), which for the time being is still too large even for modern day computers – for problems of size 𝑛 = 25, the computations needed to be done on a server with 384 GB of RAM. A different mixed integer linearization scheme, called the Kaufman-Broeckx formulation, was proposed in (Kaufman & Broeckx, 1978) with 𝑂(𝑛) additional variables. Although this is the smallest QAP linearization, its LP relaxation is known to be usually weak. This relaxation was tightened in (Xia & Yuan, 2006) using the Gilmore-Lawler bound (GLB) (Gilmore, 1962; Lawler, 1963) and further enhanced in (Zhang et al., 2013). A formulation based in the Kaufman-Broeckx family was used in (Fischetti et al., 2012) to solve (prove optimality) all the esc instances (Eschermann & Wunderlich, 1990) (including the one of size 𝑛 = 128). 2.3. Lower Bounding Techniques Exact solution of a QAP typically requires the use of a branch-and-bound framework (Anstreicher, 2003). In practice, the lack of efficiently computable, tight lower bounds for the QAP has been the key factor in the problem’s difficulty, as the tighter the bound is, the more difficult it generally is to compute. There are various approaches for obtaining lower bounds. One of the oldest methods, the Gilmore-Lawler bound (GLB), is still widely used. A comparison of older bounds based on linearization of the QAP can be found in (Karisch et al., 1999). A great success in solving previously unsolved QAP instances was achieved using the convex quadratic programming bound introduced in (Anstreicher & Brixius, 2001). A seminal breaking point in combinatorial optimization was the emergence of semidefinite programming (SDP). The SDP bounds for the QAP were first studied in (Zhao et al., 1998). The problem with this relaxation was that it involved a matrix variable of order 𝑛, and can therefore only be solved efficiently by interior point methods for, say, 𝑛 ≤ 20. This limitation has encouraged research into exploiting group symmetry of the QAP data matrices to obtain smaller and more tractable SDP problems (de Klerk & Sotirov, 2012). It has also prompted recent research into SDP relaxations of QAP where the matrix variables are of order 𝑛; see (Peng et al., 2010) and (Peng et al., 2015). In both these lines of research the authors were able to compute the best-known lower bounds for some QAPLIB instances. As we will use the lower bounding techniques for the construction of starting points for our metaheuristic, we will describe these techniques in greater detail. The computation of the GLB can be done in the following way: Denote the row vectors of matrices 𝐴 and 𝐵 by 𝑎 and 𝑏,𝑖=1,2,…,𝑛,. Let 𝑎 be the vector consisting of the (𝑛 − 1) components of 𝑎, without 𝑎,, and let 𝑏 be the vector consisting of the (𝑛−1) components of 𝑏, without 𝑏,. Define a matrix 𝐿=𝑙, as follows: 𝑙=𝑎,𝑏, 𝑖,𝑗=1,2,…,𝑛, where ⟨𝑎,𝑏⟩ is the minimal scalar product, which can be computed by ordering the vector 𝑎 nondecreasingly and 𝑏 nonincreasingly. The GLB is given by the optimal value of the n-dimensional LAP with cost matrix 𝑙,+𝑎,𝑏, : min 𝑙 ,+𝑎, 𝑏, 𝑥       (12) s.t. 𝑥,   ∈  =1 ∀ 𝑓 ∈𝑁 (13) 154 𝑥,   ∈ =1 ∀𝑖∈𝑉 (14) 𝑥 ∈{0, 1} ∀𝑖∈𝑉 ∀ 𝑓 ∈𝑁 (15) which requires only 𝑂(𝑛) computational time. The second type of the QAP lower bounds are the eigenvalue bounds. These use the fact that the set of permutation matrices Π can be characterized as: Π=𝒬𝓃∩ℰ𝓃∩𝒩𝓃, where 𝒬𝓃 is the set of orthogonal matrices, ℰ𝓃 is the set of doubly stochastic matrices and 𝒩𝓃 is the set of matrices with positive elements of size 𝑛×𝑛. The QAP can then be equivalently formulated as: min ∈𝑡𝑟(𝐴𝑋𝐵𝑋), where tr(·) is the trace of a matrix. The first eigenvalue bound that uses this QAP formulation was introduced in (Hoffman & Wielandt, 2003) and is based on the relaxation of the feasible region: min ∈𝒬𝓃tr (𝐴𝑋𝐵𝑋)=⟨λ(𝐴),λ(𝐵)⟩, where 𝜆(·) denotes the vector of eigenvalues of the matrix. This bound can be computed with very little effort but tends to be extremely weak. The improvement of this bound was done in (Hadley et al., 1992) and is called the Hadley-RendlWolkowitz (HRW): Let 𝑢 be a vector of all ones and let 𝑉 be an 𝑛×(𝑛−1) matrix with 𝑉𝑢=0 and rank(𝑉) = 𝑛 − 1. Then {𝑋∈ℛ𝓃×𝓃:𝑋𝑢=𝑋𝑢=𝑢}={ 1 𝑛𝑢𝑢 +𝑉𝑀𝑉:𝑀∈ℛ(𝓃)×(𝓃)} which can be used to reparametrize the trace formulation as: tr(𝐴𝑋𝐵𝑋)=𝑡𝑟󰇡(𝑉𝐴𝑉)𝑋(𝑉𝐵𝑉)𝑋󰇢+2 𝑛tr(𝐴𝐽𝐵)𝑋−const where 𝐽=𝑢𝑢  is a matrix of all ones, and use the eigenvalue bound to obtain the improved HRW bound: ⟨𝜆(𝑉𝐴𝑉) 𝜆(𝑉𝐵𝑉)⟩+LAP  2 𝑛𝐴 𝐽 𝐵  −const (16) The third type of the QAP lower bounds are based on a convex quadratic programming relaxation of the trace reparameterization shown above. (Anstreicher & Brixius, 2001) used the above-mentioned parametrization and showed that the following convex quadratic optimization problem gives a lower bound on the QAP: min vec(𝑋)𝑄vec(𝑋)+⟨λ(𝑉𝐴𝑉) λ(𝑉𝐵𝑉)⟩ (17) s.t. 𝑋𝑢=𝑋𝑢=𝑢 𝑋≥0, (18) where 𝑄=(𝐵⊗𝐴)−𝐼⊗𝑉𝑆 󰆹𝑉−𝑉𝑇 𝑉⊗𝐼 and 𝑆󰆹 and 𝑇 can be obtained from the spectral decomposition of 𝑉𝐴𝑉 and 𝑉𝐵𝑉. The last type of the QAP lower bounds we consider are based on a SDP relaxation developed by (Peng et al., 2015): Let (𝐵,𝐵) be a minimal trace matrix splitting of the matrix 𝐵 and compute a decomposition 𝐵=𝐵𝐵. Let 𝐵=𝐵+𝐵 the SDP relaxation model of QAP based on minimal trace matrix splitting is the following: min tr(A Y ) (19) s.t. 𝑌=𝑌−𝑌 𝑌=𝑌+𝑌 (20) 󰇧 𝐼𝐵 𝑋 𝑋𝐵𝑌 󰇨 ⪰0 󰇧 𝐼𝐵 𝑋 𝑋𝐵 𝑌 󰇨 ⪰0 (21) diag(𝑌)=𝑋diag(𝐵) 𝑌𝑒=𝑋𝐵𝑒 (22) diag(𝑌)=𝑋diag(𝐵) 𝑌𝑒=𝑋𝐵𝑒 (23) 𝑋min󰇟𝐵󰇠  ≤󰇟𝑌󰇠  ∀𝑖≠ 𝑗 (24) 󰇟𝑌󰇠 ≤𝑋max󰇟𝐵󰇠  ∀𝑖≠ 𝑗 (25) 𝑋min󰇟𝐵󰇠  ≤󰇟𝑌󰇠  ∀𝑖≠ 𝑗 (26) R. Matousek et al. / International Journal of Industrial Engineering Computations 13 (2022) 155 󰇟𝑌󰇠 ≤𝑋max󰇟𝐵󰇠  ∀𝑖≠ 𝑗 (27) 𝑋min󰇟𝐵󰇠≤󰇟𝑌󰇠  ∀𝑖≠ 𝑗 (28) 󰇟𝑌󰇠 ≤𝑋max󰇟𝐵󰇠 ∀𝑖≠ 𝑗 (29) 𝑋min󰇟𝐵󰇠  ≤󰇟𝑌󰇠  ∀𝑖≠ 𝑗 (30) 󰇟𝑌󰇠 ≤𝑋max󰇟𝐵󰇠 ∀𝑖≠ 𝑗 (31) ‖󰇟𝑌󰇠 :‖ ≤𝑋 ‖󰇟𝐵󰇠 :‖ ∀𝑖 (32) ‖󰇟𝑌󰇠 :‖ ≤𝑋 ‖󰇟𝐵󰇠 :‖ ∀𝑖 (33) ‖󰇟𝑌󰇠 :‖ ≤𝑋 ‖󰇟𝐵󰇠 :‖ ∀𝑖 (34) ‖󰇟𝑌󰇠 :‖ ≤𝑋 ‖󰇟𝐵󰇠 :‖ ∀𝑖 (35) 𝑋≥0 𝑋𝑒=𝑋𝑒=𝑒 (36) 𝐵 denotes the matrix consisting of all the off-diagonal elements of 𝐵, i.e., 𝐵 =𝐵− diag(𝑏,𝑏,…,𝑏),max (𝐵)( or min(𝐵))denotes the vector whose 𝑖-th component is the maximal element (or minimal element) in the 𝑖-th row (denoted by 𝐵,:) of 𝐵. In order to obtain a starting point (in our case, a starting permutation) for the upcoming heuristic, we use a projection of the matrix obtained from the lower bounding schemes to the space of permutation matrices: Let 𝑋 be a matrix obtained from the computation of the lower bounds. The “closest” permutation matrix 𝑋 to 𝑋 can be computed by solving the following problem: min 𝑋,−𝑋 ,      (37) s.t. 𝑋,     =1 𝑋,    =1 ∀𝑖 ∀ 𝑗 (38) 𝑋 ∈{0 ,1} ∀𝑖 ∀ 𝑗 , (39) The problem above can be reformulated into an equivalent problem using the fact that 𝑋, is binary: min 1−2𝑋, 𝑋,       (40) s.t. 𝑋,     =1 𝑋,   =1 ∀𝑖 ∀ 𝑗 (41) 𝑋 ∈{0,1} ∀𝑖 ∀ 𝑗 , (42) which is a simple LAP. 2.4 HC12 algorithm The binary HC12 algorithm, described in detail in (Matousek & Zampachova, 2011), is a stochastic heuristic searching algorithm which belongs to the class of pseudo global search methods. The basic step of the algorithm is a generation of a neighborhood of the current solution, which serves as a base for the construction of a new (improved) population. The method of generating the neighborhood is the pivotal characteristic of HC12. The paradigm of the algorithm is the search of the optimal solution in the binary (Hamming) space, that encodes the solution. In this context, it is a parallel approach to genetic algorithms, where the solution is encoded as a binary vector. The best individual of the 𝑖th generation (or iteration) is chosen as the base for the following (𝑖 + 1) generation. The approach is depicted in Fig. 1. The binary vector of the current solution is called a kernel a is denoted with an index “ker” (e.g., 𝑎). The newly generated neighborhood creates a set of 𝑐 new binary vectors 𝑎 with the same length as the vector 𝑎. These new vectors can be viewed as a population and represented by a matrix 𝐴=(𝑎,…,𝑎). The degree of locality/globality of the optimization depends on the particular way the new population 𝐴is generated. The goal of the search is to find optimal parameters 𝑥(43) with respect to the define objective function 𝑓(𝑥) on a parametric space 𝐷∈𝑁. Because of the binary representation, the parametric space is defined by a mapping Γ: {0,1} →𝐷. An important implementation detail of the mapping Γ is the translation of the binary vector from the Gray code into direct binary; afterwards, there is a problem-based decoding of the binary vector (0-1 156 problem, integer problem, or mixed integer problem). The following relationship is used 𝑥=Γ(𝑎) to denote the optimal solution as follows: 𝑥=argmin ∈ f (𝑥) (43) 𝑎=argmin ∈{ } 𝑓 Γ(𝑎) (44) Over this binary representation is defined the neighborhood relation, that describes the neighborhood s for each feasible 𝑎 as points 𝑎∈𝑠(𝑎). The choice of the neighborhood function s determines the behavior and character of the HC algorithm (Fig. 1). Fig. 1. An schematic example of the progress of the binary HC algorithm The HC12 algorithm is very effectively parallelizable. Using the neighborhood function 𝑠 (45) on a binary vector 𝑎, the population 𝐴 is generated. The set of possible neighborhood functions is denoted by 𝐻 (46). 𝑠:𝑎→𝐴 i, e, 𝑠: {0 1}→({0 1}) (45) 𝐻={𝑠 𝑠 … 𝑠} (46) The number 𝑐 of newly generated vectors in the population A0 depends on the chosen neighborhood function 𝑠 and on the length 𝑛 of the binary vector 𝑎 – it is computed as 𝑐= . For the realization of the transformations from the set 𝐻 a system of matrices 𝑀 is defined. The matrix 𝑀 corresponding to the function 𝑠 will be called a matrix of the 𝑘-th order and denoted by 𝑀. Matrix of the 𝑘-th order (𝑀) is a matrix whose rows represent all points of the Hamming metric space that are distance k from the origin (zero vector of length 𝑛): 𝑀=(0,0,⋯0,) 𝑀=⎝ ⎛ 1,0,0,⋯0, 0,1,0,⋯0, ⋮ ⋱ 0,0,0,⋯1,⎠ ⎞ 𝑀=⎝ ⎛ 1,1,0,⋯0, 1,0,1,⋯0, ⋮ ⋱ 0,0,⋯1, 1,⎠ ⎞ ⋮ 𝑀=(1,1,⋯1,) R. Matousek et al. / International Journal of Industrial Engineering Computations 13 (2022) 157 Using the 𝑘-th order matrices, the function 𝑠 can be effectively computed as: 𝑠:𝐴=𝑎⊕𝑀 Fig. 2. An example of neighbourhood generation for 4-bit binary string using transformations H = {s0,s1,s2} and matrixes M0,M1,M2, i.e. utilization in algorithms HC1, HC2 and their union HC12 From the practical point of view (because of the combinatorial expansion), only the transformations 𝑀,𝑀, and 𝑀 are used. The algorithm HC12 encodes in the last digit of its name the utilized transformations (upto order 2). The general paradigm of the HC12 algorithm is implemented using several input parameters: fun (indicator of the objective function 𝑓), nRun (the number of runs/restarts of the algorithm). The value nRun depends on the difficulty of the problem. The section of rows 6 to 10 are the HC12 algorithm itself. This part is very well (row 8) and well (row 9) parallelizable. The computations in row 8 contain the conversion from the Gray code into direct binary, implicitly. The main focus of this paper is on row 5 of the algorithm. How does one select a good starting solution? One possibility is to start at a random solution with the hope that after a sufficiently large number of tries, one does get a “good enough” solution. The other possibility is to start from a solution that is obtained by some heuristic. In this case, the heuristic in question entails the computation of the lower bound for the QAP (by one of the methods described earlier) and the projection of the obtained lower bound solution on the space of permutation matrices, by solving (40)-(42). The implementation of the HC12 algorithm for the QAP was described in (Matousek et al., 2019). As noted earlier, the optimization problem in (1)-(4) can be interpreted as a search over the space of permutation matrices 𝑋∈Π. From the problem structure of the QAP it is clear that swapping arbitrary columns of a (feasible) matrix 𝑋 always results in a different feasible matrix (and swapping the rows of the matrix has the same effect). Algorithm 1 The HC12 algorithm (Pseudo code of the general paradigm). 1: fun,nRun ← inputs 2: 𝑀←(𝑀 0 ,𝑀 1 ,𝑀 2 ) 𝑇 3: 𝑓 𝑏𝑒𝑠𝑡 ← ∞ 4: for i = [1 : nRun] do 5: 𝑎 𝑜𝑝𝑡 ← random / heuristic 6: repeat 7: 𝑎 𝑘𝑒𝑟 ←𝑎 𝑜𝑝𝑡 8: 𝐴←𝑎 𝑘𝑒𝑟 ⊕𝑀 9: 𝑎 𝑜𝑝𝑡 ←argmin 𝑎∈{0,1} 𝑛 𝑓 Γ(𝑎) 10: until 𝑎 𝑜𝑝𝑡 =𝑎 𝑘𝑒𝑟 11: 𝑓 𝑏𝑒𝑠𝑡 (𝑖)← 𝑓 󰇡Γ𝑎 𝑜𝑝𝑡 󰇢 12: 𝐴 𝑏𝑒𝑠𝑡 (𝑖,:)←𝑎 𝑜𝑝𝑡 13: endfor 14: 󰇟𝑖, 𝑓 𝑚𝑖𝑛 󰇠←min 𝑖 𝑓 𝑏𝑒𝑠𝑡 (𝑖) 15: 𝑎 𝑚𝑖𝑛 ←𝐴 𝑏𝑒𝑠𝑡 (𝑖,:) 16: 𝑥 𝑚𝑖𝑛 ←Γ(𝑎 𝑚𝑖𝑛 ) 17: return { 𝑓 𝑚𝑖𝑛 ,𝑎 𝑚𝑖𝑛 ,𝑥 𝑚𝑖𝑛 } 158 4. Results and discussion The computational experiments were carried out on 53 symmetrical QAP instances from the QAPLIB. For these instances the GLB (12)-(15), HRW (16), convex quadratic (17)-(18), and semidefinite (19)-(36) bounds, and their projections (40)- (42) were computed. For the computation of the convex quadratic (AB) and semidefinite (PE) bounds, and for the computation of the LAP for the projection, the corresponding optimization problems were implemented in JuMP environment in JULIA language and the MOSEK solver was used to obtain the solutions. The results from these computations are summarized in Table 1. For one of the instances (tho150), the AB and PE formulations were too big to handle. These computations were carried out on Intel Xeon E5530 2.40GHz CPU with 16GB of RAM. The HC12 algorithm was implemented for HPC computations on GPU CUDA 7.x (i.e., NVIDIA RTX 2080, 8GB), where not more than 6GB were used for any of the QAP instances. From Table 1 we can see that, at least in general, the more complicated formulations (convex quadratic and semidefinite) produce better (higher) lower bounds, but not necessarily better (lower) starting point values, when judged solely on the resulting projection. The trade-off is that these more complicated formulations need quite a lot more computational resources (judged by the computational time) and are only feasible for instances up to n = 100. Also, every method produced the best lower bound and best projected value for at least one problem instance. It should be noted that the lower bounds are not only useful for constructing possible starting solutions for heuristics, but also help to judge the closeness of the solution obtained by the heuristic to the true optimum. This is especially important in situation, where the is otherwise no information about what the optimal value of the QAP instance might be. Next, we used the projected values from the lower bounding techniques as the starting points for the HC12 metaheuristic and run it 1,000 times for each problem instance. The best results of these simulations (the solution with the lowest objective value out of the 1,000) are reported in Table 2. We also include the results from simulations that used random permutations as the starting point (again 1,000). Similarly, to the results of the lower bounds, there is not a clear winner, as for each of the methods (even for the random start) there are instances where it produced solutions that were better than the ones from the other methods. However, we can compare each of the bounding methods with the random start to see if there is significant difference. This comparison is summarized in Table 3 – we can see that even the “worst” performing lower bounding technique (HWB) was significantly better that random start, beating it in 30 of the 53 instances. The “best” performing lower bounding technique was the most complicated semidefinite formulation (PE), which was better than random start in 41 of the 52 instances. We can also see that the GLB method performed a bit better than the much more complicated convex quadratic (AB) one. Similar pattern can be observed for the median results reported in Table 4. The main difference is that the random start was never the best scoring method, while each of the lower bounding methods were the best in at least 6 instances. The comparison of the lower bounding method with random start for median results summarized in Table 5 shows even bigger difference than the one for best (minimum) results – the “worst” lower bounding technique (HWB, again) was better than random start in 41 of the 53 instances, and the “best” one (PE, again) was better in every one of the 52 instances. The GLB and AB methods perform similarly well. From these results, it is clear that starting a heuristic from a carefully chosen points leads to an increase in quality of the resulting solutions. The choice of the technique for constructing these starting points mainly depends on the computational resources at our disposal. While for the QAP the semidefinite (PE) formulation produced the best behaving starting points, it was also the most computationally demanding method, requiring the use of advanced convex optimization algorithms or the use of powerful solvers. In contrast to this, the GLB method produce starting points that are almost as good as the PE one, but the computational requirements for GLB are negligible. 5. Conclusion In this paper we have studied the effects of using the lower bounding techniques for the QAP as for the generation of starting points for the HC12 heuristic, that subsequently tried to find the optimal solution for the QAP. We have shown through extensive numerical computations that this utilization of the lower bounding techniques significantly improves the values of the resulting solutions. Out of the four compared lower bounding techniques, the best overall results were obtained by using the semidefinite relaxation method, which was also the most computationally demanding one. When the computational resources, or the access to high quality semidefinite optimization solvers are limited, the GLB bound can serve as an excellent surrogate – although the resulting solutions are not as good, the computational requirements are negligible. Future research will focus on extending the multicriteria and stochastic QAP instances. Also, the evaluation of various other heuristics that can use the starting points could be interesting, as different methods could benefit more (or less) from starting from an already decent point. Lastly, we expect to work on the evaluation of the starting solutions for other NP-hard optimization problems. R. Matousek et al. / International Journal of Industrial Engineering Computations 13 (2022) 159 Table 1 Lower bounds (LB), values of the projections (PV), and computational time (T) of the four considered QAP lower bounding techniques for selected QAPLIB problems (BKW – best known value). Best results (highest for LB and lowest for PV) are emphasized in bold. GLB (Gilmore, 1962) HRW (Hadley et al., 1992) AB (Anstreicher & Brixius, 2001) PE (Peng et al., 2015) Instance BKW LB PV T [s] LB PV T [s] LB PV T [s] LB PV T [s] chr12a 9552 724 44232 0.001 0 25638 0.001 0 23246 0.095 8499 25752 1.136 chr12b 9742 7146 25580 0.001 0 26712 0.001 0 27088 0.107 7340 41170 1.317 chr12c 11156 7976 18784 0.001 0 31608 0.001 0 30876 0.089 9832 40438 1.052 chr15a 9896 5625 50174 0.001 0 25958 0.001 0 25880 0.253 7441 46976 2.515 chr15b 7990 4653 54254 0.001 0 38470 0.001 0 19008 0.256 5166 30798 2.870 chr15c 9504 6165 44602 0.001 0 26144 0.001 0 35936 0.260 8808 43370 2.207 chr18a 11098 6779 89486 0.001 0 38778 0.001 0 51800 0.331 9077 78282 4.383 chr18b 1534 1534 4640 0.001 0 2372 0.001 0 3922 0.338 1534 4156 2.580 chr20a 2192 2150 9778 0.001 0 7742 0.001 0 6904 0.489 2156 10302 5.814 chr20b 2298 2196 10430 0.001 0 6418 0.001 0 7386 0.497 2237 10320 7.654 chr20c 14142 8601 79430 0.013 0 63268 0.001 0 63350 0.511 8825 69124 7.588 chr22a 6156 5924 17622 0.003 0 10178 0.001 0 10588 0.815 5964 10316 8.907 chr22b 6194 5936 13486 0.002 0 9530 0.001 0 11744 0.995 6015 10264 7.664 chr25a 3796 2765 18964 0.002 0 14186 0.001 0 13062 1.444 3244 11708 17.67 had20 6922 6166 7550 0.001 6625 7460 0.001 6671 7184 0.519 6778 7486 4.486 kra30a 88900 68360 120000 0.004 63717 117490 0.001 68467 111970 10.52 73983 118920 34.27 kra30b 91420 69065 118720 0.005 120990 140040 0.001 68876 122440 3.782 68737 127900 20.88 kra32 88700 67390 121620 0.007 59735 119450 0.002 64591 117010 4.780 72297 119540 29.96 nug18 1930 1554 2402 0.002 1663 2250 0.001 1703 2320 0.303 1753 2278 2.937 nug20 2570 2057 3080 0.002 2196 2940 0.001 2253 2884 0.481 2338 2972 4.191 nu g 21 2438 1833 3198 0.002 1979 2914 0.001 2051 3030 0.627 2215 3266 5.777 nug22 3596 2483 4598 0.002 2966 4448 0.001 3074 4346 0.816 3284 4312 9.638 nug24 3488 2676 4222 0.035 2960 4170 0.001 3024 4272 1.378 3178 4066 7.517 nu g 25 3744 2869 4728 0.003 3190 4634 0.001 3267 4534 1.521 3404 4496 8.902 nug27 5234 3701 6246 0.003 4493 6586 0.001 4604 6378 2.495 4820 6548 18.64 nug28 5166 3786 6470 0.003 4433 6310 0.001 4538 6208 2.705 4732 6274 15.07 nu g 30 6124 4539 7706 0.003 5266 7342 0.001 5360 7270 3.702 5608 7664 25.51 scr15 51140 44737 70154 0.002 10355 75950 0.001 12478 77786 0.299 46015 83280 1.339 scr20 110030 86766 203736 0.001 16113 172306 0.001 22714 193250 0.474 92426 193396 5.518 sko42 15812 11311 19522 0.008 13830 19088 0.002 14029 18748 18.23 14612 18492 121.7 sko64 48498 32522 57316 0.030 43890 56376 0.003 44513 56214 126.2 45467 56766 623.7 sko72 66256 44280 75860 0.048 60402 75772 0.004 61069 75632 222.4 61497 76082 1209 sko81 90998 60283 105932 0.061 82277 104844 0.005 83433 103642 466.3 85795 103954 2757 sko90 115534 75531 132996 0.083 105983 131398 0.006 107171 131794 1448 109260 131646 4585 sko100a 152002 98953 171886 0.091 139365 170880 0.009 140946 172554 2847 144091 171294 9946 sko100b 153890 99028 174226 0.100 141251 173878 0.007 143138 174290 3055 145783 174950 10745 sko100c 147862 95979 169888 0.101 135011 167142 0.007 136773 170340 3306 140146 169792 8920 sko100 d 149576 95921 172024 0.091 136979 167642 0.007 138736 169774 2928 140072 171152 10162 sko100e 149150 95551 171360 0.093 136996 168732 0.007 138711 168640 3090 139277 169914 9688 sko100 f 149036 96016 169768 0.093 136860 170008 0.007 138661 169090 3634 140885 169564 10627 ste36a 9526 7124 14866 0.006 0 16112 0.001 0 17454 9.257 7731 20690 60.95 ste36b 15852 8653 47768 0.007 0 42034 0.001 0 37398 8.375 12930 51936 60.06 ste36c 8.239e6 6.393e6 2.191e7 0.006 0 2.119e7 0.001 0 1.952e7 8.193 6.546e6 1.775e7 101.4 tai25a 1.167e6 962417 1.457e6 0.002 956657 1.310e6 0.001 967207 1.366e6 1.758 958027 1.411e6 12.89 tai50a 4.9386e 3.854e6 5.838e6 0.017 3.840e6 5.680e6 0.002 3.870e6 5.613e6 67.49 3.842e6 5.667e6 253.0 tai60a 7.205e6 5.55e6 8.481e6 0.27 5.537e6 8.398e6 0.004 5.575e6 8.464e6 82.67 5.544e6 8.216e6 841.2 tai80a 1.349e7 1.032e7 1.570e7 0.065 1.030e7 1.568e7 0.005 1.035e7 1.573e7 413.5 1.031e7 1.556e7 3213 tai100a 2.105e7 1.582e7 2.403e7 0.101 1.579e7 2.348e7 0.007 1.585e7 2.347e7 1667 1.552e7 2.347e7 6460 tho30 149936 90578 195698 0.003 119255 193756 0.001 124217 200194 4.067 131588 198154 27.20 tho40 240516 143804 298906 0.007 191042 303704 0.002 197661 319004 13.45 210210 313114 109.9 tho150 8.133e6 4.123e6 9.703e6 0.226 7.350e6 9.756e6 0.013 – – – – – – wil50 48816 38069 53942 0.014 45731 53420 0.003 46194 52938 40.64 46901 53754 281.4 wil100 273038 210949 293908 0.095 260827 293206 0.007 262584 291630 2591 264724 294948 11078