Efficient quantum algorithm for solving travelling salesman problem: An IBM quantum experience
Karthik Srinivasan, Saipriya Satyajit, Bikash K. Behera, Prasanta K. Panigrahi
EEfficient quantum algorithm for solving travelling salesman problem: An IBMquantum experience
Karthik Srinivasan , ∗ Saipriya Satyajit , ∗ Bikash K. Behera , and Prasanta K. Panigrahi Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India Department of Physics, Indian Institute of Technology Bombay, Mumbai 400076, India and Department of Physical Sciences, Indian Institute of Science Educationand Research Kolkata, Mohanpur 741246, West Bengal, India
The famous Travelling Salesman Problem(TSP) is an important category of optimizationproblems [1] that is mostly encountered in var-ious areas of science and engineering. Studyingoptimization problems motivates to develop ad-vanced techniques more suited to contemporarypractical problems [2–7]. Among those, especiallythe NP hard problems provide an apt platform todemonstrate supremacy of quantum over classicaltechnologies in terms of resources and time. TSPis one such NP hard problem in combinatorialoptimization [8, 9] which takes exponential timeorder for solving by brute force method. Here wepropose a quantum algorithm to solve the trav-elling salesman problem using phase estimationtechnique. We approach the problem by encodingthe given distances between the cities as phases.We construct unitary operators whose eigenvec-tors are the computational basis states and eigen-values are various combinations of these phases.Then we apply phase estimation algorithm to cer-tain eigenstates which give us all the total dis-tances possible for all the routes. After obtainingthe distances we can search through this infor-mation using the quantum search algorithm forfinding the minimum [10] to find the least possi-ble distance as well the route taken. This pro-vides us a quadratic speedup over the classicalbrute force method for a large number of cities.In this paper, we illustrate an example of the trav-elling salesman problem by taking four cities andpresent the results by simulating the codes in theIBM’s quantum simulator.
Travelling Salesman Problem (TSP) is a classic opti-mization problem [1] in the field of computer science.It belongs to an intriguing class of ‘hard’ optimizationproblems called NP hard [11, 12]. The problem involvesa salesman who has to travel N cities, visiting each cityonce and reaching ultimately at the same city where hestarted. This cycle of visiting each city once, where eachcity represents a unique vertex in a graph, and returningto the starting city is known as a Hamiltonian cycle [13].Each city is connected to other cities with a specific costassociated to each connection. The cost gives an idea of ∗ These authors contributed equally to this work. how difficult it is to take the corresponding route. Theaim of the salesman is to minimize the cost of travel, sat-isfying the above described conditions. Even if we breakthe travelling salesman problem into smaller components,each component will be at least as complex as the initialproblem. This is why it belongs to the class of NP hardproblems. The most expensive and simplest classical so-lution to the problem is to find the solution by bruteforce method. However, the problem becomes impossi-ble to solve when a large number of cities are taken. ForN cities, (N-1)! possible iterations are needed to searchfor the solution, which shoots up very fast as N increases.Other classical approaches to solve the problem includebranch and bound algorithms [14–17], heuristics [18–21]and other methods [22–24]. Using branch and bound al-gorithms the problem has been solved for around 86,000cities, but the success in branch and bound algorithmsdepends on certain factors which if not satisfied give usthe same complexity as the brute force method. Heuris-tics approach is based on providing a set of rules on op-timal selection of next city to travel. But this does notgive optimal solution in every case as heuristics result inapproximations.With the advent of the era of quantum technologiespossibilities of solving this problem with quantum com-puters has come to the limelight with the aim to tacklea much bigger problem of proving P = NP class. Severalquantum algorithms [25] have also been proposed aim-ing at the same. Goswami et al. [26] have presented aframework for efficiently solving the approximate travel-ling salesman problem. A quantum heuristic algorithmhas been proposed by Bang et al. [27] to solve the travel-ing salesman problem by generalizing the Grover search.Moylett et al. [28] have given a proof of the quantumquadratic speed up for the Travelling Salesman Problemfor bounded-degree graphs. The above mentioned algo-rithms work only when certain conditions are satisfied,however our algorithm combined with the quantum algo-rithm for finding the minimum by Durr and Hoyer [10]gives a quadratic speedup over the classical brute forcemethod without further conditions on the problem.The classical algorithms to solve the problem take in-put in the form of a matrix say A such that [ A ] ij = φ ij ,where φ ij is the cost/distance/time or any other quan-tity taken to travel from city i to city j. This quantityfor the overall journey has to be minimized. Without lossof generality, we take the quantity as cost in the currentwork. In the problem, the main motivation to take input a r X i v : . [ qu a n t - ph ] M a y as phases stems from the following two facts; first thematrix made of the given distances using the above pro-cedure is not unitary in general which implies that theimplementation and manipulation of the operator is notpossible on a quantum computer. Second, the phases willget added when we multiply them or take tensor prod-ucts of states with these phases as coefficients, that is thedistances will get added as phases which is required forthe search. Hence we represent the input as a martix B,where B ij = e i ( φ ij ) . FIG. 1. The figure illustrates the Travelling Salesman Prob-lem for four cities. The cities are represented by circles named1, 2, 3 and 4 and the cost for each path connecting the citiesare labelled on them. Here, φ ij represents the cost of travelfrom city i to city j . This is the most general case of travellingsalesman where the cost for travelling from city i to city j isnot the same as the cost for travelling from city j to city i . Next major step in our algorithm is phase estimationwhich is the backbone of our algorithm. In fact the so-lution of Travelling Salesman Problem (TSP) we presenthere seems to be a perfect application of phase estimationalgorithm. The next step in the algorithm is to constructunitaries U j from the matrix B as[ U j ] kk = 1 √ N [ B ] jk , (1)where N is the number of cities and j, k ≥ ∈ [1 , N ]. The rest of the elements in the matrix U j areinitialized to zero. It is cleary observed that U j is a diago-nal unitary matrix. Using these unitaries, we now createanother unitary matrix which is the tensor products ofall the unitaries; U = U ⊗ U ⊗ ...U N . The eigenvalues ofthis unitary matrix U, are estimated using the phase esti-mation algorithm. The phases can be easily normalizedto be bound within 0 and 2 π once we know the rangeof distances between the cities which is given to us inthe problem. Here, U is a diagonal matrix since it is atensor product of N diagonal matrices. This means thatthe eigenstates of this matrix U are computational basisstates with eigen values as the corresponding diagonalelements. Hence there will be (N-1)! diagonal elements of interest to us out of the N N elements. These elementswill have the total cost of the (N-1)! possible Hamilto-nian cycles as phases. This implies that there are (N-1)!eigenstates of U with eigenvalues being the total cost ofthe corresponding Hamiltonian cycle as phase. It is easyto calculate the position of these elements in this matrixU which will have this “useful” information. Since fora given number of cities, we know the location of theelements with total cost, we can prepare computationalbasis eigenstates corresponding to the location of each ofthese elements and extract the phase using phase estima-tion. We get the phases in form of binary output fromphase estimation algorithm, then we can easily performthe quantum algorithm for finding the minimum [10] tofind the minimum cost and the corresponding route thatis to be taken for that particular cost.In the travelling salesman problem with four citiesshown in the Fig. 1, we can take a total of six routeswhich satisfy the conditions for being a Hamiltonian cy-cle. If the cost of travelling from city i to city j is thesame as that of city j to city i (the symmetric case) thenwe will get only three routes with unique total costs. Wehave taken the symmetric case for implementing on IBMquantum experience. We have implemented this problemusing the circuit given in the Fig. 2 on the IBM quantumexperience custom topology. However in the simulationwe took six qubits for phase estimation instead of fourapart from the eigenstate qubits. We repeated this cir-cuit six times, once for each eigenstate. The eigenstatescorresponding to the different routes, the theoretical ex-pectation for the values of φ ij and the experimental ob-servations are given in the table I. TABLE I. The result of simulation for (4-1)! eigenstates andthe theoretical expectations are presented in the T.bael Thevalues taken as costs for travelling are given in the Supple-mentary Information SectionSL No. Eigenstate Theoretical Experimental1. 11000110 100100 1000002. 01101100 100100 1010003. 10001101 100000 0100004. 01110010 100000 0100005. 11100001 011000 1010006. 10110100 011000 010000
Even though the experimental results do not exactlycoincide with the theoretical expectations, we can rectifythis by taking more qubits for phase estimation, whichis its inherent property. To successfully obtain the phaseaccurate up to n bits with probability of success at least1 − (cid:15) , the number of qubits we need for phase estimation(t) is given by the following expression, t = n + (cid:38) log (cid:32) (cid:15) (cid:33)(cid:39) (2) FIG. 2. The above figure shows the quantum circuit for phaseestimation of eigenstate | (cid:105) which corresponds to theroute (Hamiltonian cycle) 1 → → → →
1. We havetaken here six qubits for estimating the phase and rest eightqubits for preparing the corresponding eigenstate (all of whichare initialized in the | (cid:105) state) by the use of pauli X gates onspecific qubits. The unitary U is the one described earlier i.e.the tensor product unitaries U ⊗ U ⊗ ...U N . The part ofthe circuit in red represents initialization of the eigenstates.The part of the circuit in blue performs the phase estimationmethod. Using the proposed algorithm, we are able to create adatabase of all possible routes that can be taken alongwith the distance of each. If one devices a quantum al- gorithm to find the minimum element in an unsorted ar-ray, which is faster than the one we currently have, thenwe can use that algorithm to find the minimum. Thisgives our algorithm a flexibility, which can be exploitedin the future to solve the travelling salesman problemmuch more efficiently.Even though our algorithm deals with a very generalcase, there are certain cases which cannot be directlysolved using our algorithm. These are the cases wherethere are restrictions on routes connecting cities. Forinstance, city i does not have a route connecting it tocity j. This can be thought of as the distance betweenthose cities being infinite. Since our algorithm requiresdistances that can be normalized such that the total dis-tance for the longest route is less than 2 π , this does notbode well for us. Fortunately, there is a way to dealwith this exception. We can take the distance betweenthe concerned cities as a very large distance such thatthe routes containing the path connecting city i and jwill have a total distance which will certainly exceed theminimum distance. Supplementary Information is available in the on-line version of the paper.
Acknowledgments
B.K.B. is financially supportedby DST Inspire Fellowship. S.S. and K.S. acknowl-edge the support of HBCSE and TIFR for conductingNational Initiative on Undergraduate Sciences (NIUS)Physics camp. We are extremely grateful to IBM teamand IBM QE project. The discussions and opinions de-veloped in this paper are only those of the authors anddo not reflect the opinions of IBM or IBM QE team.
Author contributions
K.S. and S.S. contributedequally to this work. K.S. and S.S. came up with the effi-cient algorithm to solve the TSP problem. K.S. and S.S.designed the quantum circuit and simulated using IBMquantum experience. K.S., S.S. and B.K.B. contributedto the composition of the manuscript. K.S., S.S. andB.K.B. has done the work under the guidance of P.K.P.
Author information
The authors declare no com-peting financial interests. Correspondence and requestsfor materials should be addressed to P.K.P. ([email protected]). [1] Mott, B., Job, J., Vlimant, J. R., Lidar, D. & Spiropulu,M. Solving a Higgs optimization problem with quantumannealing for machine learning.
Nature , 375–379(2017).[2] Westerlund, T., Harjunkoski, I. & Isaksson, J. Solvinga production optimization problem in a paper-convertingmill with MILP.
Comput. Chem. Eng. , 563–570 (1998).[3] Deb K. & Sinha A. Solving Bilevel Multi-Objective Opti-mization Problems Using Evolutionary Algorithms. Evol.Multi-Criterion Optim. , 110–124 (2009).[4] Wang, Z. On Solving Convex Optimization Problems withLinear Ascending Constraints. arXiv preprint 1212.4701v7 (2014). [5] Fioretto, F., Pontelli, E. & William, Y. Distributed Con-straint Optimization Problems and Applications: A Sur-vey.
J. Artif. Intell. Res. (2018).[6] Zhao, S. & Zhou, J. Solutions to Constrained OptimalControl Problems with Linear Systems.
J. Optim. TheoryAppl.
Springer (2003).[8] Dai, H., Khalil, E. B., Zhang, Y., Dilkina, B. & Song,L. Learning Combinatorial Optimization Algorithms overGraphs. arXiv preprint 1704.01665v4 (2018)[9] Hoffman, K. L. Combinatorial optimization: Current suc- cesses and directions for the future.
J. Comput. Appl.Math. , 341–360 (2000).[10] Durr, C., & Hoyer, P. A Quantum Algorithm for Findingthe Minimum. arXiv preprint quant-ph/9607014 (1996).[11] Kumlander D. NP-Hard Graph Problems AlgorithmsTesting Guidelines: Artificial Intelligence Principles andTesting as a Service.
Innovative Tech. in InstructionTechnol., E-learning, E-Assess., and Educ. Springer, Dor-drecht npj Quantum Inf. , 35 (2017).[13] Berge, Claude & Minieka, Edward Graphs and hyper-graphs. North-Holland publishing company Amsterdam (1973).[14] Lawler, E. L. & Wood, D. E. Branch-and-Bound Meth-ods: A Survey.
Operations Res. , 699–719 (1966).[15] Padberg, M. & Rinaldi, G. Optimization of a 532-citysymmetric traveling salesman problem by branch and cut. Operations Res. Lett. , 1–7 (1987).[16] Little, J. D. C., Murty, K. G., Sweeney, D. W. & Karel,C. An Algorithm for the Traveling Salesman Problem. Operations Res. , 972–989 (1963).[17] Held, M. & Karp, R. M. The Traveling-Salesman Prob-lem and Minimum Spanning Trees. Operations Res. ,1138–1162 (1970).[18] Lin, S. & Kernighan, W. An Effective Heuristic Algo-rithm for the Traveling-Salesman Problem. OperationsRes. , 498–516 (1973).[19] Geem, Z. W., Kim, J. H. & Loganathan, G. V. A new heuristic optimization algorithm: Harmony search. Simu-lation , 60-68 (2001).[20] Laporte, G., Martello, S. The selective travelling sales-man problem. Discrete Appl. Math. , 193-207 (1990).[21] Karg, R. L. & Thompson, G. L. A Heuristic Approach toSolving Travelling Salesman Problems. Manage. Sci. ,225–248 (1964).[22] Bhide, S., John, N., & Kabuka, M. R. A real -time so-lution for the traveling salesman problem using a booleanneural network. IEEE Int. Conf. Neural Networks Conf.Proc.
Appl. Soft.Comput. , 867-878 (2014).[24] Rana, S. & Srivastava, R. S. Solving Travelling SalesmanProblem Using Improved Genetic Algorithm. Indian J.Sci. Technol. , (2017).[25] Kieu, T. D. The Travelling Salesman Problem and Adi-abatic Quantum Computation: An Algorithm. arXivpreprint quant-ph/1801.07859 (2018).[26] Goswami, D., Karnick, H., Jain, P. & Maji, H. K. To-wards Efficiently Solving Quantum Traveling SalesmanProblem. arXiv preprint quant-ph/0411013 (2004).[27] Bang, J., Yoo, S., Lim, J., Ryu, J., Lee, C. & Lee, J.Quantum heuristic algorithm for traveling salesman prob-lem. J. Korean Phys. Soc. , 1944–1949 (2012).[28] Moylett, D. J., Linden, N. & Montanaro, A. Quantumspeedup of the traveling-salesman problem for bounded-degree graphs. Phys. Rev. A , 032323 (2017). SUPPLEMENTARY INFORMATION: EFFICIENT QUANTUM ALGORITHM FOR SOLVINGTRAVELLING SALESMAN PROBLEMI. 4 CITY EXAMPLE EXPLANATION
The distance matrix for the Travelling salesman problem consisting of four cities is given by D = 12 e iφ e iφ e iφ e iφ e iφ e iφ e iφ e iφ e φ e iφ e iφ e iφ (3)where φ = φ = π/ , φ = φ = π/ , φ = φ = π/ , φ = φ = π/ , φ = φ = π/ , φ = φ = π/ , and the unitaries U j with j = 1, 2, 3 and 4 are as follows U = | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | (4) U = e iφ | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | (5) U = e iφ | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | (6) U = e iφ | (cid:105)(cid:104) | + e iφ | (cid:105)(cid:104) | + e φ | (cid:105)(cid:104) | + | (cid:105)(cid:104) | (7)We constructed the unitaries U j by decomposing it as follows; U j = e ia e ib e ic
00 0 0 e id = (cid:34) e i ( c − a ) (cid:35) ⊗ (cid:34) e ia e ib (cid:35) e i ( d + c − a − b ) (8)Note that this is the unitary U j , by putting specific values of a, b, c, d we can find decomposition for each U , U , U , U . For phase estimation we need controlled-( U ⊗ U ⊗ U ⊗ U ) which is same as ( C − U ⊗ C − U ⊗ C − U ⊗ C − U ), where C − U represents controlled- U . For realizing controlled U j we implemented controlled eachelement in the decomposition of U j in equation 8. FIG. S1. The figure illustrates our implementation of controlled U j gate with the use of single qubit unitaries U(a) as describedearlier. Here we implement a controlled unitary (C-A) which is equivalent to implementing controlled unitaries that belong tothe decomposition of A. In this figure, x = (d+c-a-b)/2 II. THE EIGENSTATES
Fig. 2 shows the part of the circuit where we estimate the phase corresponding to the route going through cities1, 2, 3 and 4 in the same order and returning to one. Starting from any other city but following the same orderas mentioned will also give us the same eigenstate. The eigenstate for any route can be calculated as follows. In aparticular route, if we are going from city i to j , then each i is uniquely mapped to j . Hence we can write i as afunction of j i.e, i ( j ). The eigenstate corresponding to that particular route is, | ψ (cid:105) = ⊗ j | i ( j ) − (cid:105) (9)where j goes from 1 to n.Once these eigenstates are calculated, circuits similar to Fig. 1, with eigenstate qubits initiated to the rest of theeigenstates, can be run in parallel. Then we can search through this database using the quantum search algorithm inthe order of O ( (cid:112) ( N − O ( (cid:112) (( N − /
2) steps.
III. SUBROUTINES IN IBM QUANTUM EXPERIENCE - CUSTOM TOPOLOGY
In the simulation, we took advantage of the “Add subroutine” under advanced option present in the custom topology[1] and built these unitaries. The qasm code written for the simulation of the circuit (Fig. S2) is given below. i n c l u d e ” q e l i b 1 . i n c ” ;qreg q [ 1 4 ] ; c r e g c [ 6 ] ; − U1 , C − U2 , C − U3 & C − U4gate uni ( a , b , c , d ) x , y , z { cu1 ( c − a ) x , y ;u1 ( a ) x ; cu1 ( b − a ) x , z ;ccx x , y , z ; cu1 ( ( d − c+a − b ) /2) x , z ;ccx x , y , z ; cu1 ( ( d − c+a − b ) /2) x , y ;cu1 ( ( d − c+a − b ) /2) x , z ; } gate bigU ( a , b , c , d , e , f , g , h , i , j , k , l ) m, n , o , p , q , r , s , t , u { uni ( 0 , a , b , c ) m, n , o ; uni ( d , 0 , e , f ) m, p , q ;uni ( g , h , 0 , i ) m, r , s ; uni ( j , k , l , 0 ) m, t , u ; } − U with the v a l u e s given gate finU m, n , o , p , q , r , s , t , u { bigU ( p i /2 , p i /8 , p i /4 , p i /2 , p i /4 , p i /4 , p i /8 , p i /4 , p i /8 , p i /4 , p i /4 , p i /8) m, n , o , p , q , r , s , t , u ; } − Uˆ2gate finU2 m, n , o , p , q , r , s , t , u { finU m, n , o , p , q , r , s , t , u ;finU m, n , o , p , q , r , s , t , u ; } − Uˆ4gate finU4 m, n , o , p , q , r , s , t , u { finU2 m, n , o , p , q , r , s , t , u ;finU2 m, n , o , p , q , r , s , t , u ; } − Uˆ8gate finU8 m, n , o , p , q , r , s , t , u { finU4 m, n , o , p , q , r , s , t , u ;finU4 m, n , o , p , q , r , s , t , u ; } − Uˆ16gate finU16 m, n , o , p , q , r , s , t , u { finU8 m, n , o , p , q , r , s , t , u ;finU8 m, n , o , p , q , r , s , t , u ; } − Uˆ32gate finU32 m, n , o , p , q , r , s , t , u { finU16 m, n , o , p , q , r , s , t , u ;finU16 m, n , o , p , q , r , s , t , u ; } { finU f , n , o , p , q , r , s , t , u ;finU2 e , n , o , p , q , r , s , t , u ; finU4 d , n , o , p , q , r , s , t , u ;finU8 c , n , o , p , q , r , s , t , u ; finU16 b , n , o , p , q , r , s , t , u ;finU32 a , n , o , p , q , r , s , t , u ; } { cu1 ( − p i /2) x , y ; } gate r3 x , y , z { cu1 ( − p i /4) x , z ; r2 y , z ; } gate r4 w, x , y , z { cu1 ( − p i /8) w, z ; r3 x , y , z ; } gate r5 v , w, x , y , z { cu1 ( − p i /16) v , z ; r4 w, x , y , z ; } gate r6 u , v , w, x , y , z { cu1 ( − p i /32) u , z ; r5 v , w, x , y , z ; } h q [ 1 ] ;h q [ 2 ] ; h q [ 3 ] ;h q [ 4 ] ; h q [ 5 ] ;x q [ 6 ] ; x q [ 8 ] ;x q [ 9 ] ; x q [ 1 1 ] ;fU q [ 0 ] , q [ 1 ] , q [ 2 ] , q [ 3 ] , q [ 4 ] , q [ 5 ] , q [ 6 ] , q [ 7 ] , q [ 8 ] , q [ 9 ] , q [ 1 0 ] , q [ 1 1 ] , q [ 1 2 ] , q [ 1 3 ] ; h q [ 0 ] ;r2 q [ 0 ] , q [ 1 ] ; h q [ 1 ] ;r3 q [ 0 ] , q [ 1 ] , q [ 2 ] ; h q [ 2 ] ;r4 q [ 0 ] , q [ 1 ] , q [ 2 ] , q [ 3 ] ; h q [ 3 ] ;r5 q [ 0 ] , q [ 1 ] , q [ 2 ] , q [ 3 ] , q [ 4 ] ; h q [ 4 ] ;r6 q [ 0 ] , q [ 1 ] , q [ 2 ] , q [ 3 ] , q [ 4 ] , q [ 5 ] ; h q [ 5 ] ;measure q [ 0 ] − > c [ 0 ] ; measure q [ 1 ] − > c [ 1 ] ;measure q [ 2 ] − > c [ 2 ] ; measure q [ 3 ] − > c [ 3 ] ;measure q [ 4 ] − > c [ 4 ] ; measure q [ 5 ] − > c [ 5 ] ; TSP qasm.py