A pre-sales person is covering a lot of cities and would like to travel to each city once and finally back to the originating city in such a way as to minimize the cost of fuel given that the more distance covered, the greater is the fuel consumed. This is the Traveling Salesman Problem (TSP). It is very popular since it is very hard to get an exact solution when the number of cities is large.

This is the second of a series of blogs that will talk about Variational Quantum Eigensolver (VQE). In this blog, we will show how to map the TSP to the Ising Model that will enable us to use the VQE to come up with an approximate solution to the TSP.

Let’s visualize this using the sample diagram below.

What you see is a graph representation of cities as nodes. Edges between nodes represent a route between them. Each edge has a weight which stands for the distance between the cities connected by this edge.

Suppose our pre-sales guy/girl wants to optimize the cost of fuel. He/She will need to list all possible routes and compute the total distance for each route. A route here is a closed route which begins with a certain city, visits the remaining cities once and returns to the starting city. Here is a listing of all possible routes (and total distance of that route) given the above graph:

From the result above, we can see that the shortest path has a total distance of 18 units.

## Cost to Minimize

The closed routes in the Traveling Salesman problem are called Hamiltonian routes. It is specified by binary variables , where is the ith city and is the order of the city for a given instance of a route.

For example, a route

is represented by

The rows represent the city and the columns represent the position of the node in the current route. To determine the position of a city in a given route, find the row of the city and find the column that contains the value 1. Each cell of the matrix is a variable and has a value of either 1 or 0.

**Note**: Although we did not state explicitly, the route ends in the starting city. So in this example, after visiting city 0, the traveller will go back to city 2.

Notice that in the above matrix, we have the following constraints:

and

We can “unpack” this matrix into a vector such that , where is the number of cities. So in our example, we have and the variable , maps to . An instance of this vector is a route. Using this representation, the route

which is represented by

becomes the vector

Later, we will see that this is an eigenstate when we map the TSP to the Ising Model.

If we let be the distance between city i and city j, then the total distance for each route is given by

If we add the constraints for the rows and columns, we get

where the terms in underbraces are squared so that the lowest possible minimum value is zero. The factor can be chosen so that the constraints will be respected.

## Mapping to the Ising Model

Our objective in this blog is to map this problem to the Ising Model so that it’s amenable to be solved using Quantum Computing.

First, we map the routes to the qubits:

Next, we map the cost function above to the Ising Model using the following change of variables:

where is the Pauli Z operator.

See Appendix to see a sample of how terms expand.

By mapping this to the Ising Model, this becomes an eigenvalue problem

where is the Ising Hamiltonian and is the eigenstate corresponding to the eigenvalue . The eigenstate corresponds to the state . By finding the eigenstate that corresponds to the minimum eigenvalue, we will be able to solve our problem.

We can easily do this mapping using the IBM QISKit tsp module.

## An Example Using QISKit

For a TSP problem with N cities, the number of qubits is equal to . The number of basis states corresponding to this is . For a 4-city problem, the number of qubits is and the number of states is . Because of this, we just give an example of 3 cities.

Suppose we have the following distance matrix for 3 cities:

Using IBM QisKit, we can map this to the Ising model using the `get_tsp_qubitops`

method of the tsp class:

import numpy as np import qiskit.aqua.translators.ising.tsp as tsp from qiskit.aqua.algorithms import VQE, ExactEigensolver from qiskit.aqua.translators.ising.tsp import TspData w=[[ 0., 78., 45.],[ 78., 0., 100.],[ 45., 100., 0.]] np.array(w) tspdata=TspData("tmp", 3, [], w=np.array(w)) qubitOp, offset = tsp.get_tsp_qubitops(tspdata)

The variable `qubitOp`

is an Operator object and represents the Hamiltonian of the system (not to be confused by the Hamiltonian circuit of the TSP).

Here is the code to solve this problem using the ExactEigensolver (lifted from the Jupyter notebook of the TSP tutorial in QISKit):

ee = ExactEigensolver(qubitOp, k=1) result = ee.run() """ algorithm_cfg = { 'name': 'ExactEigensolver', } params = { 'problem': {'name': 'ising'}, 'algorithm': algorithm_cfg } result = run_algorithm(params,algo_input) """ print('energy:', result['energy']) x = tsp.sample_most_likely(result['eigvecs'][0]) print('feasible:', tsp.tsp_feasible(x)) z = tsp.get_tsp_solution(x) print('solution:', z) print('solution objective:', tsp.tsp_value(z, tspdata.w))

Executing the code above using gives a solution to the TSP:

energy: -600111.5 feasible: True solution: [0, 1, 2] solution objective: 223.0

## Conclusion

We have not yet discussed the Variational Quantum Eigensolver (VQE). That will be the last topic in this series. What we have done is understand how the TSP is mapped into the Ising Problem. We then used the ExactEigensolver to find the minimum eigenvalue that will correspond to a solution to the problem. In the next blog, we will show how the VQE can be used to get an approximation of a solution to the TSP problem after mapping it to the Ising Model.

## Appendix

Mapping the TSP to the Ising model involves transforming the cost function

using

where is the Pauli Z operator.

As an example, for , the second term above will look like:

We can expand the square as

For a given p, say p=0, we have:

Here is the code to generate the circuit for the above expansion (taken from the tsp source code):

for p in range(num_nodes): for i in range(num_nodes): for j in range(i): print("p= ", p, " i=", i, "j=",j) shift += penalty / 2 zp = np.zeros(num_qubits, dtype=np.bool) zp[i * num_nodes + p] = True print(i*num_nodes + p) pauli_list.append([-penalty / 2, Pauli(zp, zero)]) zp = np.zeros(num_qubits, dtype=np.bool) zp[j * num_nodes + p] = True print(j*num_nodes + p) pauli_list.append([-penalty / 2, Pauli(zp, zero)]) zp = np.zeros(num_qubits, dtype=np.bool) zp[i * num_nodes + p] = True zp[j * num_nodes + p] = True pauli_list.append([penalty / 2, Pauli(zp, zero)])

If you print the pauli_list, you will get the following:

IIIIIZIII -> Z3 IIIIIIIIZ -> Z0 IIIIIZIIZ -> Z3Z0 IIZIIIIII -> Z6 IIIIIIIIZ -> Z0 IIZIIIIIZ -> Z6Z0 IIZIIIIII -> Z6 IIIIIZIII -> Z3 IIZIIZIII -> Z3Z6

This matches our expansion above.