The traveling salesman problem is a well-known NP-hard problem in combinatorial optimization supported by an extensive body of research [1]. Given a fully connected graph of n nodes and a n × n adjacency matrix M whose entries give the cost to travel from node i to node j and M _ii = 0, the goal is to output a roundtrip (cycle) that minimizes the cost of traversing all the nodes in the graph (each node is traversed exactly once and so has exactly one ingress and one egress). We restrict ourselves to the case when the cost to travel from node i to node j is same as the cost to travel from node j to node i (negative values of cost are permissible) and is further given by M _ij + M _ji . If we let X be a n × n symmetrical indicator matrix consisting of 0 s and 1 s (x _ij = 1 if node i is adjacent to node j in the output cycle; we will also call this as the nodes being connected in the output cycle or graph; this is not to be confused with whether the nodes are connected in the input graph which they always are), we seek to find the X that minimizes ∑M. ∗ X where .∗ denotes element-wise multiplication of matrices and ∑ denotes sum over elements of the matrix.

This paper describes a quantum solution to this problem that leverages adiabatic quantum computation. Adiabatic quantum computation is one of the two paradigms of quantum computation—gate-based model being the other one. In adiabatic quantum computation, a system of qubits is prepared in an initial ground state associated with the Hamiltonian H ₀ and the system is gradually evolved to a target state whose Hamiltonian is given by H ₁. If the time evolution is done slowly enough, the theorem of adiabatic quantum computation guarantees that the qubits will track the ground state of the evolving Hamiltonian and at the end they will be found in the ground state corresponding to H ₁. Mathematically, if H t = 1 − t H 0 + t H 1 (1) where t is some normalized unit of time, then H(0) = H ₀, H(1) = H ₁. H ₀ is the starting Hamiltonian (initial conditions) and at time t = 1 the qubits will assume values that minimize the energy associated with H ₁.

Up until now the discussion is general. D-Wave’s ¹ quantum computers are special purpose devices that aim to realize adiabatic quantum computation through quantum annealing where the Hamiltonian is constrained to be an Ising Hamiltonian. The energy (eigenvalues) associated with the target state H ₁ is given by: E I s i n g = − ∑ i h i s i − ∑ i ∑ j , j ≠ i J i j s i s j (2)

The Ising constraint is a natural consequence of how the device is built and wired. It is a system of superconducting qubits obeying the Ising model. In Eq. 2, h _i are biases and J _ij are coupling constants. s _i are the spins of the qubits s _i ∈ {−1, + 1}. The H ₀ state is realized by subjecting the system to a high transverse magnetic field in the x direction (H ₀ = X ⊗ … ⊗ X where X is Pauli spin matrix along x) which aligns all the spins in the x direction and thus sets the ground state to an equal superposition of all states in the Z basis and the external magnetic field is gradually turned off to attain the state corresponding to H ₁. The spins will automatically align themselves so as to minimize the energy given by Eq. 2. This process is known as quantum annealing.

A simple change of variables (x = (1 + s)/2) together with the observation that x ² = x when x ∈ {0, 1} transforms Eq. 2 into (ignoring constant term as it is of no consequence) E = x T Q x (3) where Q is a Hermitian matrix (symmetric if dealing with real numbers only). The minimization of Eq. 3 is a Quadratic Unconstrained Binary Optimization (QUBO) problem. D-Wave’s devices are tailored to solve this class of problems.

The objective function of the traveling salesman problem can then be written as: f = ∑ M . ∗ X = ∑ i = 0 i < n ∑ j = i + 1 j < n M i j + M j i x i j = ∑ k = 0 k < m B k b k (4) where M _ij + M _ji is the cost of traveling from i to j (or vice-versa) and we have used indices starting from 0 which are convenient for turning equations into computer code. The x _ij variables are mapped to b _k and m is given by ⁿ C ₂ = n (n − 1)/2. Thus for a problem of size n, we have to solve for m = ⁿ C ₂ unknowns. f is a linear function in b _k and the solution is trivial if there are no constraints on b _k (set b _k = 1 if B _k < 0 and 0 otherwise). However, there are constraints on what b _k can be. Each node is traversed only once and connects to only two other nodes in the output cycle. These constraints can be written as n equations (one equation for each i): ∀ i ∈ 0 , … , n − 1 ∑ j x i j = 2 (5)

These constraints have to be mapped to b _k (as we solve for b _k ultimately) and the function mapping i, j indices to k (and thus x _ij to b _k ) is given by: k ( i , j , n ) = undefined f o r i = j k ( j , i , n ) f o r i > j i n − i ( i + 1 ) / 2 + j − ( i + 1 ) f o r i < j (6)

So far what we have is a linear integer programming problem subject to linear constraints. To convert it into a QUBO, we formulate each constraint as: ∀ i ∈ 0 , … , n − 1 min 2 − ∑ j x i j 2 (7) and use the method of Lagrange multipliers to absorb the constraints into the objective function as follows: q = ∑ k = 0 k < m B k b k + ∑ i = 0 i < n λ i 2 − ∑ j , j ≠ i b k i , j 2 (8) where λ _i are Lagrange multipliers. The minimization of Eq. 8 is a QUBO problem. The Ising (and equivalently QUBO) formulation of many NP problems can be found in [2]. The QUBO class of problems is NP-hard.

We have used D-Wave’s 5,436 qubit Advantage 1.1 QPU [3] to solve instances of Eq. 8. We created 8 sample problems all of size n = 8 (corresponding to m = 28 binary variables) in which the entries of M were randomly sampled from [1, n] or {1, 2, …, 8}. All the Lagrange multipliers (λ _i ) were set equal to each other and further equal to the maximum of the absolute values in the matrix M (thus 8 in this case). The Lagrange multiplier acts as a weight given to the constraint. It should be set high enough to ensure the constraint is satisfied, but setting it too high obscures the real function we are trying to minimize. In addition to the Lagrange multipliers, there are some other parameters that can affect the solutions returned by the QPU such as the chain strength and annealing time. We use the default annealing time of 20 μs of the QPU. The chain strength parameter is related to minor embedding of the problem graph into the actual physical connectivity graph of the QPU. The qubits represented by the variables b _k are logical qubits of the problem and we assume in our problem formulation that any of them can interact with another. However, in the physical processor the qubits are laid out in a certain fixed pattern (graph) and to emulate the problem graph it has to be minor embedded into the physical graph. This is a hard problem in its own right (see [4] and references cited therein) and in general given m QUBO variables, we require O(m ²) physical qubits for minor embedding. That’s why even though the processor has 5,436 qubits, its compute power translates to the order of 5436 ≈ 73 logical qubits. The task of minor embedding the problem graph into the QPU is akin to compilation of a classical computer program to the instruction set of the CPU. We use the default function provided by D-Wave (minorminer.find_embedding) ² to calculate the minor embedding of our problem graph. Once, the problem graph is embedded into the physical graph, a chain of physical qubits maps to a logical qubit. Long chains should be avoided and can create problems when the constituent qubits of a chain assume different values - this is known as a chain break and can result in garbage solutions. The chain strength parameter is a parameter that avoids the chain from breaking. Like the Lagrange multiplier, it should be set high enough but setting it too high will hide the actual problem we want to solve. In our implementation we use the scaled ³ function provided by D-Wave to automatically set the chain strength without any manual ad-hoc fine tuning. With this, we found the auto-computed chain strength to be between 30.5 and 32 for all the problem cases and there were no chain breaks in the returned solutions. If the annealing process were perfect, the wavefunction would be in the ground state of target Hamiltonian in the end, but in practice that is not always the case. Therefore we perform multiple runs (samples)—100 runs for each problem in our case. We then iterate through the samples (observations) in order of lowest energy and select the first sample that satisfies the constraints of Eq. 5. This marks the solution to the problem by the quantum method.

To compare the results and performance of the quantum solution, we also solved each of the problem cases by trying out all possible combinations (brute-force). For a problem of size n, the total number of possible combinations is given by (n − 1)!/2. For n = 8 this translates to 2,520. Figure 1 shows graph of the original objective function ∑M. ∗ X against all the combinations for a sample problem. It has many local extrema and no recognizable structure. All the 8 problems are characterized by closely spaced extrema - the difference in the objective function between the best and second-best input configuration ranges from 1–3. This equals the energy difference δE between the ground and first-excited state of the quantum system. Quantum annealing performs best when this energy gap is high. When the energy gap is small, the system can jump to a excited state over course of the time evolution giving suboptimal solution in the end.

A classical optimizer for the traveling salesman problem was also developed. It uses the following simple heuristic: we start with a random initial configuration (cycle). Then, we pick 2 nodes at random and check if swapping their positions (order) would lead to an improvement of the objective function. If so, we swap and repeat the process. If we have tried all ⁿ C ₂ swaps and got no improvement of the objective function, the algorithm terminates. This algorithm is guaranteed to converge (and converges fast) but it does not always converge to the global optimum. The chances of converging to the global optimum decline with increase in size of the problem n.

Table 1 shows our results and compares the performance of all the 3 solvers. Overall we find the quantum solver performed worse than the classical solver. In terms of time too, the classical solver ran faster for n = 8. The quantum solver takes 20 μs × 100 samples = 2 ms to run and this time excludes the time spent in minor embedding of the problem graph. In comparison the classical solver implemented in Python took less than a ms to run for each of the problems on a 2018 Macbook Pro equipped with a 2.2 GHz Intel Core i7 processor and 16 GB 2400 MHz DDR4 RAM. n = 8 is also the maximum size of the problem the quantum solver was able to handle. When we tried problems with n = 9 (m = 36), the D-Wave QPU did not return any valid solution for 5 out of 8 test cases with the methods described in previous section, and for n = 10 (m = 45), the QPU returned a solution for only 1 out of 8 test cases. The solutions were also suboptimal compared to classical solver. They are not listed in the main paper but can be found in the Supplementary material accompanying the paper. One strength of our implementation is that no tweaks are needed to the code and the same code that works when entries of M are sampled from [1, n] also works as-is when the entries are sampled from [−a, +a].

One advantage of the classical solver is that it always respects the constraints of the problem (Eq. 5) whereas the quantum solution turns the constraints into soft-constraints that can be violated and have to be checked for satisfaction. Another advantage is that the quantum solver can return incorrect solutions in which there is more than 1 loop (cycle) e.g., it might configure 4 nodes to be in one loop and 4 others to be in another loop—there is no constraint in Eq. 5 preventing that. Each node is still connected to two other nodes. The classical solver does not suffer from this limitation.

For reference, D-Wave also provides a code sample to solve the TSP using a QPU. ⁴ We tried that code as well but it resulted in errors. The formulation provided here is better than D-Wave’s solution because it takes into account that x _ij = x _ji (consequently there are ⁿ C ₂ variables to solve for) whereas in D-Wave’s solution this constraint is treated as a soft-constraint using Lagrange multipliers and the total number of variables is n(n − 1). The classical optimizer we have described is the 2-opt algorithm for solving the TSP [5]. Several computer codes and packages for classical solutions to TSP can be found online. ^{5 ,  6 ,  7} A quantum study using D-Wave to solve the related Capacitated Vehicle Routing Problem (CVRP) can be found in [6].

Lastly, in this paper we restricted ourselves to the case when the cost of traveling from i to j is same as the cost of traveling from j to i. The restriction can be removed and the methods developed here can be applied to the general case as well. See the Appendix for details.

In conclusion, this paper presented a QUBO formulation of the traveling salesman problem which makes it amenable to quantum annealers and although our experiments with D-Wave QPU did not demonstrate any quantum advantage, we remain hopeful that with advancements in technology the quantum solution might outperform classical method esp. with large n although quantum annealers like the D-Wave QPU will find a challenge since the number of qubits required grows as m ² due to the minor embedding problem. Hence, simply adding more qubits to the QPU while keeping the rest of the architecture and design the same is unlikely to succeed. Better innovations will be required.