Solving combinatorial problems with QAOA
This notebook demonstrates the process of solving a combinatorial problem described in the QUBO formulation using the QAOA algorithm.
The Quantum Approximate Optimization Algorithm is an algorithm of the family of Variational Quantum Algorithms, and thus requires three main components to be executed:
A cost function encoding the problem to solve, in this modeled by a cost hamiltonian \(H_C\).
An Ansatz describing the architecture of the parametrized quantum circuit, with a cost hamiltonian \(H_C\) and a mixer hamiltonian \(H_M\)
A classical optimizer use to find the best parameters of the circuit
In this example we will be tackling the traveling salesman problem (TSP) which aims at finding the cycle (path where the last vertex is the first one) with the lowest cost in a graph while visiting each vertex once. We will do this by using an approach based on QAOA.
[1]:
import networkx # used here to generate the graph
import numpy as np
import numpy.typing as npt
import random
Solving the Travelling Salesman Problem consists in finding the shortest hamiltonian path, a path where we visit each vertex of the graph only once, and return to the starting point. We chose to model such path as a list of indices of the vertices. First we will be generating the graph. For this example we will be taking a 3-vertex graph with positive weights.
[2]:
nbr_vertexs = 3
G = networkx.fast_gnp_random_graph(
nbr_vertexs, 1
) # generates n vertexs with every vertex connected to each other
cost_list = []
for u, v, w in G.edges(data=True):
w['weight'] = random.randint(3, 10)
cost_list.append(([u, v], w['weight']))
cost_list.append([[v, u], w['weight'] + random.randint(-2, 15)])
print(G)
print()
print(cost_list)
Graph with 3 nodes and 3 edges
[([0, 1], 5), [[1, 0], 12], ([0, 2], 5), [[2, 0], 20], ([1, 2], 3), [[2, 1], 7]]
Now that we have our graph with random weights, we need to get the adjacency matrix of this graph.
This is a matrix that represent all of the weights between each vertices. We define the cost of a path to be the sum of the distances between each vertex of the path, plus the distance to come back to the starting point.
[3]:
def generate_adjacency_matrix(
list_of_cost: list[tuple[list[int], int]], nb_cities: int
):
"""Example of list of cost = [[[0,1],23], [[1,2], 52], [[0,2], 12]]"""
w_ij = np.zeros((nb_cities, nb_cities), dtype=int)
for i in list_of_cost:
cost = i[1]
w_ij[i[0][0], i[0][1]] = cost
return w_ij
adj_matrix = generate_adjacency_matrix(cost_list, 3)
print(adj_matrix)
[[ 0 5 5]
[12 0 3]
[20 7 0]]
Now that we have our matrix we will be generating our QUBO expression. A QUBO (quadratic unconstrained binary optimization) is a quadratic polynomial expression of binary variables, it is used to represent different combinatorial problems.
The goal is to find the set of binary values that minimize the expression.
For example if we have the expression : \(2*x_0 - 3*x_1\\\) Then the minimum value will be -3 with \(x_0 = 0\) and \(x_1 = 1\).
For each path between two vertexs we need to define a QUBO monomial so we have 2*edges.
[4]:
from mpqp.execution.vqa.qubo import Qubo, QuboAtom
def generate_qubo_expression(adj_matrix: npt.NDArray[np.int64]) -> Qubo:
length = len(adj_matrix)
result = 0
for i in range(length):
for j in range(length):
if adj_matrix[i][j] == 0:
continue
current = adj_matrix[i][j] * QuboAtom(f'x{i},{j}')
if isinstance(result, int):
result = current
else:
result += current
return (
result if not isinstance(result, int) else 0 * QuboAtom("x")
) # if the adjacency matrix is empty
expr = generate_qubo_expression(adj_matrix)
print(expr)
initial_coeffs = expr.get_terms_and_coefs()
print(initial_coeffs)
5*x0,1+5*x0,2+12*x1,0+3*x1,2+20*x2,0+7*x2,1
[(5, ['x0,1']), (5, ['x0,2']), (12, ['x1,0']), (3, ['x1,2']), (20, ['x2,0']), (7, ['x2,1'])]
Now that we have our Qubo expression we need to implement the constraints of this problem, constraints here
The first one is the fact that we have to go through each cities exactly once, in this case we express it as xor operators between every path that end with the same vertex.
For example if the path \(x_{0,1}\) is chosen then the path \(x_{2,1}\) cannot be chosen because they both lead to the same vertex.
The second constraint is the inverse because it is not possible to start from the same city twice so once again the path \(x_{0,1}\) excludes the path \(x_{0,2}\).
The way penalties works is when minimizing the loss function (our QUBO) we can enforce some behavior so that the minimum of the function respects the rules of the problem. Here because we need to go through each cities once so it translates to xor operations with a negative coefficient hence if the xor is true it will decrease the value of the loss function.
Note : The coefficient of penalties can decide if the minimization finds the correct answer, if the penalties are too high then the minimization won’t find the difference between solutions as long as they respect the constraints. However, if the penalties are not high enough the minimum of the loss function can ignore the penalties.
[5]:
def add_constraints(expr: Qubo, nbr_vertexs: int) -> Qubo:
for i in range(nbr_vertexs):
if i != 0:
current = QuboAtom(f'x{i},0')
else:
current = QuboAtom(f'x{i},1')
for j in range(1, nbr_vertexs):
if i == j or (i == 0 and j == 1):
continue
path = QuboAtom(f'x{i},{j}')
print(f'{current.value} ^ {path.value}')
expr += -20 * (current ^ path)
for i in range(nbr_vertexs):
if i != 0:
current = QuboAtom(f'x0,{i}')
else:
current = QuboAtom(f'x1,{i}')
for j in range(1, nbr_vertexs):
if i == j or (i == 0 and j == 1):
continue
path = QuboAtom(f'x{j},{i}')
print(f'{current.value} ^ {path.value}')
expr += -20 * (current ^ path)
return expr
expr = add_constraints(expr, nbr_vertexs)
x0,1 ^ x0,2
x1,0 ^ x1,2
x2,0 ^ x2,1
x1,0 ^ x2,0
x0,1 ^ x2,1
x0,2 ^ x1,2
[7]:
from mpqp.execution.vqa.qaoa import QaoaMixer, QaoaMixerType, qaoa_solver
from mpqp.execution.devices import IBMDevice
from mpqp.execution.vqa import Optimizer
def test_results(res: str, coeffs: list[tuple[float, list[str]]]):
current_answer = 0
other = 0
for i in range(len(res)):
if res[i] == '0':
other += coeffs[i][0]
else:
current_answer += coeffs[i][0]
return current_answer <= other
mixer = QaoaMixer(QaoaMixerType.MIXER_X)
result = qaoa_solver(expr, 3, mixer, IBMDevice.AER_SIMULATOR, Optimizer.POWELL, [1.] * 6)
print(result)
print()
print(
f"The path chosen : {result.final_state} is the best path : {test_results(result.final_state, initial_coeffs)}"
)
Minimum cost: -96
Associated state: 011001
Associated values: {'x0,1': 0, 'x0,2': 1, 'x1,0': 1, 'x1,2': 0, 'x2,0': 0, 'x2,1': 1}
The path chosen : 011001 is the best path : True