Solve the Travelling salesman problem

title

Travelling salesman problem

The traveling salesman problem for a graph $G = (V, E)$, where each edge $uv$ in the graph has a weight $W_{uv}$ associated to it, is to find the Hamiltonian cycle such that the sum of the weights of each edge in the cycle is minimized. The decision problem (does a path of total weight ≤ W exist?) is NP-complete.

How to make the TSP Quantum

To solve this problem we define binary variable $x_{ij}= \left\{ \begin{array}{ll} 1 & \text{if depot $i$ is at location $j$}\\ 0 & \text{otherwise} \end{array} \right.$

A path can then be expressed as a matrix: $X=\begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}$ that represents a path $0 \rightarrow 2 \rightarrow 1$.

The total length of the path is: $$D(x) = B\cdot \sum_{j=0}^{N - 1} \sum_{i,k=0}^{N - 1} w_{ik} x_{i,j} x_{k, j+1}$$

Attempting to minimize the function above will most likely give us a wrong solution, where we could be at two depots at the same location, or not even visit all depots. To solve this issue, we will add the constraints above as additional penalties to the distance function D(x). This can be achieved constructing a cost function C(x) as below: $$C(x) = A\cdot\sum_{j=0}^{N-1}\left(\sum_{i=0}^{N-1} x_{i j}-1\right)^{2}+A\cdot\sum_{i=0}^{N-1}\left(\sum_{j=0}^{N-1} x_{i j}-1\right)^{2}+ D(x)$$

$B$ shall be small enough that it is never favorable to violate the constraints; one such constraint is $0 < B \cdot \text{max}(W_{uv}) < A$ (we assume in complete generality $W_{uv} ≥ 0$ for each $(uv) ∈ E$, and $B=1$). If the traveling salesman does not have to return to his starting position, we can restrict the sum over $j$ from 0 to N−2 in $D(x)$ ( we dont care if last position and first position are connected). We can optimize the Hamiltonian as we may fix node 1 to appear first in the cycle.

Formulating TSP as Ising model

To solve this problem, we use a Hamiltonian $H = H_A + H_B$ , with $H_A$ the Hamiltonian given for the directed (or undirected) Hamiltonian cycles problem. $H_A = A\cdot\sum_{j=0}^{N-1}\left(\sum_{i=0}^{N-1} x_{i j}-1\right)^{2}+A\cdot\sum_{i=0}^{N-1}\left(\sum_{j=0}^{N-1} x_{i j}-1\right)^{2}$
The actual objective to minimize $H_B = B\cdot\sum_{j=0}^{N - 1} \sum_{i,k=0}^{N - 1} w_{ik} x_{i,j} x_{k, j+1}$ adds the weights for each edge.

We can use the Ising Model as a formulation.

title

Given graph $𝐺 = (𝑉, 𝐸)$ with edge weights $𝑊_{uv}$, find the hamiltonian cycle with minimum sum of edge weights:
$H=\alpha \sum_{j=0}^{N-1}\left(\sum_{i=0}^{N-1} x_{i j}-1\right)^{2}+\alpha \sum_{i=0}^{N-1}\left(\sum_{j=0}^{N-1} x_{i j}-1\right)^{2}+\beta \sum_{j=0}^{N - 1} \sum_{i,k=0}^{N - 1} w_{ik} x_{i,j} x_{k, j+1}$

$i$: number of depot
$j$: position in the cycle
$x_{ij} = 1$ if depot $i$ is located at position $j$ in the cycle, else 0
$w_{ik}$: weight on edge $(ik)$
$\alpha$: the scaling for the constraints (A1), (A2)
$\beta$: the scaling for the objective (B)

The qubo function consists of three terms, 2 terms in $H_A$ and one quadratic term in $H_B$:
$H_A = \alpha \sum_{j=0}^{N-1}\left(\sum_{i=0}^{N-1} x_{i j}-1\right)^{2}+\alpha \sum_{i=0}^{N-1}\left(\sum_{j=0}^{N-1} x_{i j}-1\right)^{2}$

The first term $\sum_{j=0}^{N-1}\left(\sum_{i=0}^{N-1} x_{i j}-1\right)^{2}$ is to model that each location in the cycle must be assigned to one unique depot (A1). As we will minimize the qubo function, we would like to get the lowest possible value, which is 0. This is true if $\left(\sum_{i=0}^{N-1} x_{i j}-1\right)=0$. Therefore $\sum_{i=0}^{N-1} x_{i j}=1$, this implies $x_{i j}=1$ for exactly one pair $(i, j)$ and 0 otherwise.

The second term $\sum_{i=0}^{N-1}\left(\sum_{j=0}^{N-1} x_{i j}-1\right)^{2}$ ensures that each depot is exactly on one location in the cycle (A2) (analog to first term).

Note that $(\sum_i x_{ij} - 1)^2 = - \sum_i x_{ij} + 2 \sum_{i > k} x_{ij} x_{kj} + const$

The third term $H_B = \sum_{j=0}^{N - 1} \sum_{i,k=0}^{N - 1} w_{ik} x_{i,j} x_{k, j+1}, \beta=1$ adds the weights for each edge $(ik)$ (B).

Now, we can map the problem of minimizing $C(x)$, to a quantum analogue by promoting the boolean variable $x_{ij}$ to a quantum operator as: $x_{ij} \rightarrow \frac{1-\hat{Z}_{ij}}{2} \text{, with }\hat{Z}_{ij}=I\bigotimes I\bigotimes...\bigotimes\sigma^z_{ij}\bigotimes...\bigotimes I$

With this notation then if qubit $i,j$ is in state $|1⟩$ the solution is at location $j$ with depot $i$. Using this operator we can transform the cost function $C(x)$ into a Hamiltonian for the system of $n^2$ qubits.

Use qss.algorithms.qubo to implement the QUBO function with our SDK.

Let's implement $H_B=\beta \sum_{j=0}^{N - 1} \sum_{i,k=0}^{N - 1} w_{ik} x_{i,j} x_{k, j+1}, \beta=1$

In the next step we implement $H_A=\alpha_1 \sum_{j=0}^{N-1}\left(\sum_{i=0}^{N-1} x_{i j}-1\right)^{2}+\alpha_2 \sum_{i=0}^{N-1}\left(\sum_{j=0}^{N-1} x_{i j}-1\right)^{2}$.

Note that $(\sum_i x_{ij} - 1)^2 = - \sum_i x_{ij} + 2 \sum_{i > k} x_{ij} x_{kj} + const$

$\alpha$ and $\beta$ must be balanced. If $\alpha$ is too big, the optimizer will try to minimize the contraints and forget about the cost objective $H_B$. If $\alpha$ is too small, then the optimizer will find solutions that violate constraints (A1) or (A2).

Let's pose a constraint that the agent must start from the 0-th point and goes to the first one: $-10(x_{00} + x_{11})$

Compute the Hamiltonian

Now we create our Hamiltonian opertator, that we will use to solve the problem by finding its lowest eigenvalue (energy state)

In the next step we use the VQE Variational Quantum Eigensolver algorithm to adapt the angles in the ansatz circuit so that the minimal eigenvalue of the hamiltonain is returned. in this case we use the state vector for an exact calculation of the eigenvalues.
We choose one appropriate optimizer method, to solve the eigenvalue problem, e.g.

With the following test program you can compare all optimizer methods against each other:

This takes some time 4-6h for 'cuqs', 10-15h for 'qss' !

Sample run of the above cell

NONE    87s per run
0 [-95.99999710751007, -95.99999981436119, -79.99999939882011, -79.99999941892904, -69.49999983440878, -79.9999997333704, -69.99999878620403, -79.99999991996103, -95.99999682436939, -79.999999935115]
CG      100s per run
1 [-95.99999987653527, -95.99999999790253, -79.99999999315659, -75.99999967108741, -69.49999991856231, -95.99999957124055, -75.99999949840787, -79.99999999466465, -95.9999995011978, -79.99999999618097]
COBYLA  82s per run
2 [-79.99979220645096, -79.99995290919352, -95.9997689253631, -79.99978443382881, -95.99996476231684, -79.99983333578456, -95.99948804494534, -75.99926352608807, -95.99959878974275, -69.4924319524998]
BFGS    86s per run
3 [-95.99999998907293, -95.9999999992776, -79.99999999504831, -79.99999999940411, -69.49999999967466, -79.99999999620968, -69.99999999359255, -79.99999999951663, -95.99999999927931, -79.99999999882874]
SLSQP   37s per run
4 [-95.99962121273015, -79.99975977217105, -75.99982034414984, -85.99974083277789, -63.99901346094046, -81.99951715606446, -95.99966475567095, -69.99965451555488, -95.9992291716867, -75.99989611630568]
Powell  35s per run
5 [-95.95319254022507, -95.97353711991252, -95.94687354505373, -95.87459624709459, -79.90218124648075, -69.99652045764691, -85.99263923618861, -79.98331034194987, -95.8990149557406, -79.99414466987284]
TNC     250 - 850s / per run
6 [-74.33845122673526, -93.80198517706091, -49.05059776278679, -34.06108509207752, -0.7835182515148205, -56.33127363103466, -46.39010288890995, -40.03917644229523, -50.13753622769087, -45.1888775536836]
nelder-mead  680-1030s per run
7 [-79.99959834377123, -95.99999308588454, -79.99916418508184, -95.99712996679304, -85.99280093244582, -82.82259365967829, -95.99858601491891, -95.99758087458108, -58.99963594376695, -70.00067087534984]

DWave implementation

References:

Ising formulations of many NP problems, Andrew Lucas, https://arxiv.org/pdf/1302.5843.pdf