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.