Optimization with Quantum Annealing
In this tutorial, we will explore how to use Quantum Annealing to solve a simple optimization problem using QiliSDK.
Note
If you haven’t already, it might be useful to check out these tutorials first: Quantum Basics, Quantum Circuits and Quantum Annealing.
The Problem
Say we have four people and we need to form two teams of two people each. We have some information about how well each pair of people work together, and we want to form the teams such that the total “compatibility” of the teams is maximized. The information about who likes who is as follows:
Alice |
Bob |
Carol |
Dave |
|
|---|---|---|---|---|
Alice |
N/A |
1 |
3 |
4 |
Bob |
1 |
N/A |
5 |
2 |
Carol |
3 |
5 |
N/A |
6 |
Dave |
4 |
2 |
6 |
N/A |
So this means that Alice would prefer to work with Dave, while Bob would prefer to work with Carol, and so on.
Mathematically, we can represent our two teams as binary variables, where \(x_i\) is 0 if person \(i\) is in team 0, and 1 if they are in team 1. The compatibility of team 0 can then be expressed as:
Where \(p_{ij}\) is the compatibility score between person \(i\) and person \(j\). Note that here each term is only non-zero if both \(x_i\) and \(x_j\) are 0, which means that both people are in team 0. Meanwhile the compatibility of team 1 can be expressed as:
Thus we want to maximize the total compatibility \(C = C_0 + C_1\), subject to the constraint that team 1 (and thus team 0) has exactly two people:
One issue with this formulation is that it has a constraint, which makes it more difficult to solve using quantum optimization algorithms. To fix this, we can use some techniques to turn this into a QUBO (Quadratic Unconstrained Binary Optimization) problem. In this case, since we just have a single equality constraint, we can add a squared version of it to the objective function as a penalty term, such that the penalty term is zero when the constraint is satisfied, and positive otherwise. For more details on this process, check out the QUBO reference. Now we can assume that the problem has been reformulated into the following unconstrained optimization problem:
Where \(c_{ij}\) are coefficients that encode the objective function as well as the constraint.
The Solution
Now that we have our problem formulated as a QUBO, we can use quantum annealing to find an approximate solution to it. Quantum annealing is a metaheuristic optimization algorithm that is inspired by the process of annealing in metallurgy, where a material is heated and then slowly cooled to remove defects and find a low-energy state. In quantum annealing, we start with a simple Hamiltonian whose ground state is easy to prepare and then slowly evolve it into a more complex Hamiltonian that encodes the optimization problem we want to solve. The hope is that if we do this slowly enough, the system will remain in its ground state throughout the evolution, and thus end up in the ground state of the final Hamiltonian, which corresponds to the optimal solution of our problem.
Assuming we have reformulated our problem into a QUBO, we can construct the problem Hamiltonian for quantum annealing, which is in the form:
Where \(Z(i)\) is the Pauli-Z operator acting on qubit \(i\), and \(c_{ij}\) are the coefficients from our QUBO formulation. Note that these variables \(Z(i)\) are related to the binary variables in our original problem by the transformation \(x_i = (1 - Z(i))/2\),
Our mixing Hamiltonian is typically chosen to be the transverse field Hamiltonian, which is given by:
The overall time-dependent Hamiltonian that we evolve is then given by:
Where \(A(t)\) and \(B(t)\) are functions that determine how we interpolate between the mixing Hamiltonian and the problem Hamiltonian over time. For simplicity here will just assume that we do a linear interpolation, such that \(A(t) = 1 - t/T\) and \(B(t) = t/T\), where \(T\) is the total annealing time.
The Implementation
To simulate the evolution of this time-dependent Hamiltonian, first we need to form our problem:
from qilisdk.core.model import QUBO, ObjectiveSense
from qilisdk.core.variables import BinaryVariable, EQ
num_people = 4
vars = [BinaryVariable(f"x{i}") for i in range(num_people)]
preferences = [[0, 1, 3, 4],
[1, 0, 5, 2],
[3, 5, 0, 6],
[4, 2, 6, 0]]
model = QUBO("team_formation_example")
team_1 = sum(
preferences[i][j] * vars[i] * vars[j]
for i in range(num_people)
for j in range(i+1, num_people)
)
team_0 = sum(
preferences[i][j] * (1 - vars[i]) * (1 - vars[j])
for i in range(num_people)
for j in range(i+1, num_people)
)
model.set_objective(
team_0 + team_1,
label="obj",
sense=ObjectiveSense.MAXIMIZE
)
model.add_constraint(
"team_size_constraint",
EQ(sum(vars[i] for i in range(num_people)), 2),
lagrange_multiplier=10
)
Then we use the model to form our problem Hamiltonian, and we can also define our mixing Hamiltonian and the schedule for our evolution:
from qilisdk.analog import X
from qilisdk.analog import Schedule
from qilisdk.core import QTensor
problem_hamiltonian = model.to_hamiltonian()
mixer_hamiltonian = sum(-1.0 * X(i) for i in range(num_people))
We also need to define the schedule - how the coefficients of the problem and mixing Hamiltonians change over time:
T = 100.0
schedule = Schedule(
hamiltonians={"problem": problem_hamiltonian, "mixer": mixer_hamiltonian},
coefficients={
"mixer": {(0, T): lambda t: 1 - t / T},
"problem": {(0, T): lambda t: t / T},
},
dt=0.1,
)
We will start in the ground state of our mixing Hamiltonian, which is the equal superposition state over all possible team assignments:
initial_state = QTensor.tensor_product([QTensor.ket(0) + QTensor.ket(1) for _ in range(num_people)]).normalized()
Finally, we initialize our quantum simulator, execute the evolution, and read out the results:
from qilisdk.functionals import AnalogEvolution
from qilisdk.readout import Readout
from qilisdk.backends import QiliSim
backend = QiliSim()
evolution = AnalogEvolution(schedule=schedule, initial_state=initial_state)
readout = Readout().with_sampling(1000).with_expectation([problem_hamiltonian])
results = backend.execute(evolution, readout)
print("Results:", results)
This will print something like the following:
Functional Results: [
Sampling Results: (
nshots=1000,
samples={'0110': 513, '1001': 487}
)
Expectation Value Results: (
expectation_values=[8.999696720325046],
)
]
As you can see from these results, the samples satisfy the constraint (i.e. have exactly two 1’s). More specifically, the two samples “0110” and “1001” (corresponding to Alice/Dave on a team and Bob/Carol on the other team) are the optimal solutions to our problem, with 9 being the highest obtainable compatibility score.