 S. Podolsky, Yu. Zorin, 2015 112 ISSN 1681–6048 System Research & Information Technologies, 2015, № 2 UDC 004.89:004.4 O(1) DELTA PART COMPUTATION TECHNIQUE FOR THE QUADRATIC ASSIGNMENT PROBLEM S.V. PODOLSKY, YU.M. ZORIN The quadratic assignment problem is rightfully considered to be one of the most challenging problems of combinatorial optimization. Since this problem is NP-hard, the use of heuristic algorithms is the only way to find in a reasonable time a solution that is close to optimal. One of the most effective heuristic algorithms is the Robust Tabu Search, which is the basis of many subsequent metaheuristic algorithms. The paper describes a novel approach to scanning the neighborhood of the current solu- tion that allows reducing by half the number of delta values that were required to be computed with complexity )(NO in most of the heuristics for the quadratic assign- ment problem. Using the correlation between the old and new delta values, obtained in this work, a new formula of complexity )1(O is proposed. The results obtained leads up to 25% performance increase as compared to such well-known algorithms as the Robust Tabu Search and others based on it. The formula obtained in this paper may be successfully applied to other heuristics using a full scan of the solution neighborhood. INTRODUCTION The quadratic assignment problem (QAP) was first mentioned by Koopmans and Beckmann in 1957 [1] and still remains one of the most challenging combinato- rial optimization problems. The problem formulation is as following. There are N locations and N facilities. Distances between each pair of locations and flow values, i.e. number of transportations, are provided. The goal is to assign all the facilities to different locations so that to minimize the sum of all distances multi- plied by the corresponding flows. This can be formulated in a form of the objec- tive function    N i N j jiij fd 1 1 )()(min  , (1) where )(x is a facility number assigned to location with number x ; ijd is a dis- tance between locations i and ;j )()( jif  is a flow between facilities )(i and )( j . Since the problem is NP-hard [2] there is no exact algorithm that can solve the QAP in polynomial time. Furthermore, the travelling salesman problem (TSP) may be seen as a special case of the QAP in the case when all the facilities are connected via flows having a constant value into a single ring while the other flows have zero value. Thus the QAP can be considered even more challenging than the TSP. The only methods that allow obtaining feasible solutions for the QAP in- stances of size 30 and higher are heuristics. One of the most efficient heuristic algorithm for the QAP is the Robust Tabu Search (Ro-TS) developed by Eric O(1) delta part computation technique for the quadratic assignment problem Системні дослідження та інформаційні технології, 2015, № 2 113 Taillard in 1991 [3], which produces high quality solutions even for very large instances. Among other successful heuristics are genetic algorithms [4], ant sys- tems [5] and others. Typically, predominant majority of them are based on the same approach of the neighbor solutions representations obtained by pairwise ex- change of elements in the solution vector. The purpose of this study is to propose a technique of the current solution neighborhood scanning that allows obtaining a significant performance increase of the metaheuristic algorithm for the quadratic assignment problem solution. BACKGROUND OF THE NEIGHBORHOOD SCANNING The process of the neighborhood scanning described for the Ro-TS is one of the most representative environment subject to improvement, since the Ro-TS per- forms initial solution construction and objective function (1) evaluation only once at the beginning of the algorithm before main iterations started. During the initial stage, a random solution vector is generated and its cost is computed by the objec- tive function as defined by (1). In each main iteration of the Ro-TS a move which relies on exchange of the two elements from the current solution vector is per- formed. A pair of elements to be exchanged is selected among all 2NC pairs using ij minimization criterion. Actually the selection criterion is more complicated than just delta value minimization, but for now we omit its discussion. After the exchange a new cost can be obtained from the previous cost by adding a ij term, where all 2NC values ij are computed on the initial stage of the Ro-TS, each with the time complexity O(N), before main iterations start  ))(())(( )()()()()()()()( jiijjiijiijjjjiiij ffddffdd  )].)(()()[( )()()()()()()()( ,,1 gigjjgigigjg N jigg gjgi ffddffdd     It is convenient to store all computed values ij in a jagged matrix and ex- tract them by indices i and j once the solution cost must be modified after the elements i and j exchange performed. However, each time after the exchange of two elements, all 2NC values ij become unfeasible and needed to be corrected. It is obvious, that after elements p and q are exchanged, the new pq value is evaluated trivially as pq reflecting the reverse exchange of the same pair. For all those ij which do not involve the two elements p and q that were exchanged just before, i.e. }Ø{},{},{ qpji  , new values *ij can be evaluated exploiting their old values ij with the complexity )1(O as following  ))(( )()()()()()()()(* jqiqipjpqiqjpjpiijij ffffdddd  ))(( )()()()()()()()( qjqipipjiqjqjpip ffffdddd   , ., jig  (2) S. Podolsky, Yu. Zorin ISSN 1681–6048 System Research & Information Technologies, 2015, № 2 114 However, as it was observed by Taillard [3], for those ij which involve one of the indices p or q from the last exchange, i.e. 1},{},{ qpji  all new values *ij must be evaluated anew with complexity )(NO by (2). This technique of complete delta recomputation was widely used since it was introduced by Taillard and proceeded without any improvements in further researches such as Reactive Tabu Search [6]. Therefore it is still used in miscellaneous recent Ro-TS implementations such as efficient heuristic for sparse matrices [7] etc. This paper shows how to bypass this issue and evaluate with the linear time complexity only a half of deltas that involve nodes from the last pair exchange. However, it would be important to prove that this issue is a real performance bottleneck and needed to be resolved. RELEVANCE AND PERFORMANCE As it was mentioned before, the total number of all possible exchanges equals N NNCN )1(2  . Among them there are possible pair exchanges that involve only one of the two just exchanged elements and the total number of them equals )2(2 N . The last number is obtained from the considerations that each of the two elements being swapped can be exchanged with any element from the rest of 2N other distinct elements, as it shown in Fig. 1. Consequently, the total num- ber of those ij which can be corrected with the complexity )1(O is equal to the difference of all possible exchanges number and the total number of exchanges being corrected with complexity :)(NO .22 )1)(4()2(22 )1(  NNNNN So, what is to be done is to compute only )2( N instead of )2(2 N new delta values with complexity )(NO and the rest ones computed with .)1(O How- ever, the main question still remains — whether this will affect iteration perform- ance significantly? The answer was obtained after the Ro-TS open source C++ code profiled with the use of MS Visual Studio Profiling tools on Tai100a QAP instance. … k 21 ji N – 2 Fig. 1. An illustration of further exchange of one of the two exchanged elements O(1) delta part computation technique for the quadratic assignment problem Системні дослідження та інформаційні технології, 2015, № 2 115 As it is shown in Fig. 2, the compute_delta method which computes del- tas anew takes over 50% of total samples, being the most expensive call, while compute_delta_part takes only 29% of total samples to update each of oth- er delta values with the )1(O time complexity. Therefore, we suggest that compu- tational time per iteration can be reduced up to a quarter with the use of proposed technique, especially on large QAP instances. NEW O(1) DELTA PART COMPUTATION TECHNIQUE EXPLORATION Suppose we have three elements with indices kji ,, in the solution vector and the values jkikij  ,, are already known. Our purpose is to exchange a pair of ele- ments i and j and compute new values *** ,, jkikij  afterwards. First of all let’s consider the assignments* of the flows )()()()()()( ,, gkgjgi fff  to the dis- tances kgjgig ddd ,, (Fig. 3), where by g we denote an arbitrary element from the rest of the other 3N elements distinct from .,, kji Though in terms of the QAP an assignment usually means the assignment of facility to specific location, here we imply the assignment of flow to distance caused by a pair of regular QAP assignments. All the values jkikij  ,, include the following terms, which in- dicate the cost change caused by reassignments of flows to distances to (from) element g after corresponding pair exchanges: Fig. 2. A sample profiling report for the Ro-TS on Tai100a instance *Though in terms of the QAP an assignment usually means the assignment of facility to specific loca- tion, here we imply the assignment of flow to distance caused by a pair of regular QAP assignments. S. Podolsky, Yu. Zorin ISSN 1681–6048 System Research & Information Technologies, 2015, № 2 116 ,)()()()()()()()( gijggjiggjjggiigij fdfdfdfd   ,)()()()()()()()( gifkggkfiggkfkggifigik dddd   .)()()()()()()()( ijkggkjggkkggjjgjk fdfdfdfd   Now, let’s consider the assignments of the same three flows to the same three distances after the elements i and j are exchanged (see Fig. 4). The new values *** ,, jkikij  that we need to compute after this exchange will contain the following terms ,* ijij   ,)()()()()()()()(* gjkggkiggkkggjigik fdfdfdfd   .)()()()()()()()(* gikggkjggkkggijgjk fdfdfdfd   igd gif  jgd gjf  kgd gkf  i j k g N–3 Fig. 3. Assignments of the flows to the distances before any exchanges Fig. 4. Assignments of the flows to the distances after i and j exchanges igd gif  jgd gjf  kgd gkf  i j k g N–3 O(1) delta part computation technique for the quadratic assignment problem Системні дослідження та інформаційні технології, 2015, № 2 117 It is important to note, that exactly these 3N terms *ik and *jk for each arbi- trary element g require the complexity )(NO for each new *ik and *jk values evaluation. After a simple analysis of the right-hand side expressions for ,, ikij  ** ,, jkikjk  the following dependency between them has been discovered .** ijjkikjkik   (3) This means that it is actually not necessary to compute both new values *ik and *jk simultaneously. It is enough to compute anew only one of them and the second one can be obtained with the use of the first. However, it is worth reminding that we have considered only those terms jkikij  ,, that connect our elements with each of the rest of the other 3N elements. So let’s consider now the terms ** ,,,, jkikjkikij RRRRR which denote the cost change caused by assignments of the flows to distances between the ele- ments .,, kji The variables jkikij RRR ,, indicate the cost change before ele- ments i and j are swapped, and the variables ** , jkik RR indicate the cost change after the exchange of i and j respectively ,ijijij R  ,ikijik R  ,jkjkjk R  ,*** ikikik R  .*** jkjkjk R  The values ** ,,,, jkikjkikij RRRRR are expressed by the following formulae:  ))(())(( )()()()()()()()( jkikkjkikjkijkikij ffddffddR  ;))(())(( )()()()()()()()( ijjijiijjjiijjii ffddffdd   (4)  ))(())(( )()()()()()()()( jikijkjiijikkjijik ffddffddR  |;))(())(( )()()()()()()()( jkkjikkijjkkkkii ffddffdd   (5)  ))(())(( )()()()()()()()( kjijijikjkjijikijk ffddffddR  ;))(())(( )()()()()()()()( ikkijkkjiikkkkjj ffddffdd   (6)  ))(())(( )()()()()()()()(* ijkjjkjijijkkjijik ffddffddR  S. Podolsky, Yu. Zorin ISSN 1681–6048 System Research & Information Technologies, 2015, № 2 118 ;))(())(( )()()()()()()()( kiikkiikiikkkkii ffddffdd   (7)  ))(())(( )()()()()()()()(* jikiikijijikkiijjk ffddffddR  .))(())(( )()()()()()()()( kjjkkjjkjjkkkkjj ffddffdd   (8) Hence, we can represent the equality (3) as following .)()()()()( **** ijijjkjkikikjkjkikik RRRRR  Now we can compute a new value *jk with the complexity )1(O using the computed anew value *ik and the old values jkikij  ,, .**** jkikijjkikikijjkikjk RRRRR  (9) After substitution of the expressions (4-8) for ** ,,,, jkikjkikij RRRRR into (10) which are the parts of (9) without deltas, we will obtain a cumbersome ex- pression (Appendix A) which probably has no benefits on small QAP instances because of a bunch of arithmetical operations .** jkikijjkik RRRRR  (10) To get rid of such complicated representation, an attempt of simplification has been made with the MATLAB script using simplify embedded function (Appendix B). As a result the final compact formula has been obtained to com- pute all the new values *jk with the time complexity .)1(O  ** ikijikjkjk ).)(( kjkijkjiikijkjkijkjiikij ffffffdddddd  EXPERIMENTAL RESULTS In order to verify the proposed approach the sample problems provided by E.Taillard and R.Kothari [8] have been chosen. The comparison of the following Taillard’s algorithms: the Ro-TS, the Fast ant system (FANT), Tabu Search with Local Optimization (TsLo) with the proposed algorithm (PzTabu) has been done. These algorithms have been chosen because of their high repute and availability of their source codes. All algorithms have been run 50 times on problems with dimensions varying from 20 to 81. The goal of the first part of the simulation was to define best and average solution costs obtained by all algorithms. The second part of simulation has been aimed on comparison of the average time and itera- tions number required by the algorithms to reach the result not worse than the best one obtained on the first stage. It may be considered as some measure of the algo- rithm convergence speed. All collected data is shown in Table. O(1) delta part computation technique for the quadratic assignment problem Системні дослідження та інформаційні технології, 2015, № 2 119 T a b l e . Statistical results of comparative analysis of the algorithms for the QAP solution Problem Algo-rithm Best cost Mean cost Mean time to reach the best cost (ms) Mean number of iterations to reach the best cost Ro-TS 703482 704449 1353 60175 FANT 708654 709059 3394 116938 TsLo 703482 710428 9481 17345 Tai20a PzTabu 703482 704296 870 39938 Ro-TS 1818442 1823408 6035 83488 FANT 1841298 1894374 68786 107317 TsLo 1825384 1825626 28016 16987 Tai30a PzTabu 1820934 1821489 3055 48741 Ro-TS 7265144 7269162 269441 1353407 FANT 7350644 7351273 553179 1759429 TsLo 7388361 7396978 777978 65547 Tai60a PzTabu 7270382 7274858 159796 879714 Ro-TS 91062 91086  241003 560526 FANT 91106 91323  289203 720510 TsLo 91596 91666 323767 15944 Sko81 PzTabu 91030 91061 135472 360150 From Table we can make sure that the proposed approach outperformed all compared algorithms in terms of operation speed and convergence. It should be noted that, as expected, the larger the size of the problem the greater is the advan- tage of the proposed approach. A relatively small number of iteration of the TsLo algorithm may be explained by the fact that its each iteration performs local opti- mization procedure which takes longer time. CONCLUSIONS The formula obtained in this paper can also be successfully applied to the other heuristics that use the whole neighbor solutions scanning such as Reactive Tabu Search [6], heuristic for sparse matrices [7] and others. The Ro-TS algorithm is the most representative and exploitable because the solution construction is per- formed only once and the rest of computational time is dedicated for neighbor solutions scanning and delta values update. As computational experiments proved, we can get significant performance increase in comparison with well- known algorithms including the Ro-TS by replacing a half of )(NO computation operations with the operations of )1(O time complexity. APPENDIX A. The final formula’s part before simplification  ** jkikijjkik RRRRR  ))(())(( )()()()()()()()( jikijkjiijikkjij fdd fddff   ))(())(( )()()()()()()()( jkkjikkijjkkkkii ffddffdd  S. Podolsky, Yu. Zorin ISSN 1681–6048 System Research & Information Technologies, 2015, № 2 120  ))(())(( )()()()()()()()( kjjiijikjkjijiki ffddffdd   ))(())(( )()()()()()()()( ikkijkkjiikkkkjj ffddffdd   ))(())(( )()()()()()()()( jkikkjkikjkijkik ffddffdd   ))(())(( )()()()()()()()( ijjijiijjjiijjii ffddffdd   ))(())(( )()()()()()()()( ijkjjkjijijkkjij ffddffdd   ))(())(( )()()()()()()()( kiikkiikiikkkkii ffddffdd   ))(())(( )()()()()()()()( jikiikijjiikkiji ffddffdd  ).)(())(( )()()()()()()()( kjjkkjjkjjkkkkjj ffddffdd   APPENDIX B. Matlab formula’s simplification script clc clear all echo off syms fii fjj fkk fij fji fik fki fjk fkj syms dii djj dkk dij dji dik dki djk dkj Rij = ... (dik - djk) * (fik - fjk) + ... % Missed g = k in delta ij (dki - dkj) * (fki - fkj) + ... % Missed g = k in delta ij (dii - djj) * (fii - fjj) + ... % Loopback (dij - dji) * (fij - fji); % Reverse flows direction Rik = ... (dij - dkj) * (fki - fji) + ... % Missed g = j in delta ik (dji - djk) * (fik - fij) + ... % Missed g = j in delta ik (dii - dkk) * (fkk - fjj) + ... % Loopback (dki - dik) * (fjk - fkj); % Reverse flows direction Rjk = ... (dki - dji) * (fij - fkj) + ... % Missed g = i in delta jk (dik - dij) * (fji - fjk) + ... % Missed g = i in delta jk (djj - dkk) * (fkk - fii) + ... % Loopback (dkj - djk) * (fik - fki); % Reverse flows direction R_ik = ... (dij - dkj) * (fkj - fij) + ... % Missed g = j in delta* ik (dji - djk) * (fjk - fji) + ... % Missed g = j in delta* ik (dii - dkk) * (fkk - fii) + ... % Loopback (dik - dki) * (fki - fik); % Reverse flows direction O(1) delta part computation technique for the quadratic assignment problem Системні дослідження та інформаційні технології, 2015, № 2 121 R_jk = ... (dji - dki) * (fki - fji) + ... % Missed g = i in delta* jk (dij - dik) * (fik - fij) + ... % Missed g = i in delta* jk (djj - dkk) * (fkk - fjj) + ... % Loopback (djk - dkj) * (fkj - fjk); % Reverse flows direction x = - Rik - Rjk + Rij + R_ik + R_jk; simplify(x) REFERENCES 1. Koopmans T.C., Beckmann M.J. Assignment problems and the location of economic activities // Econometrica. — 1957. — № 25. — P. 53–76. 2. Sahni S., Gonzalez T. P-complete approximation problems // Journal of the Associa- tion of Computing Machinery. — 1976. — № 23. — P. 555–565. 3. Taillard E. Robust tabu search for the quadratic assignment problem // Parallel Computing. — 1991. — № 17. — P. 443–455. 4. Tate D.M., Smith A.E. A genetic approach to the quadratic assignment problem // Computers and Operations Research. — 1995. — № 22. — P. 73–83. 5. Colorni A., Maniezzo V. The ant system applied to the quadratic assignment problem // IEEE Transactions on Knowledge and Data Engineering. — 1999. — 11. — P. 769–778. 6. Battiti R., Tecchiolli G. The reactive tabu search // ORSA Journal on Computing. — 1994. — № 6. — P. 126–140. 7. Paul G. An efficient implementation of the robust tabu search heuristic for sparse quadratic assignment problems // European Journal of Operational Research. — 2011. — № 209. — P. 215–218. 8. Kothari R., Ghosh D. Insertion based Lin-Kernighan heuristic for single row facility layout // Computers & Operations Research. — 2013. — № 40(1). — P. 129–136. Received 17.07.2013 From the Editorial Board: the article corresponds completely to submitted manu- script.