Solving QUBO problems with QAOA in Qiskit¶
In this jupyter notebook, we are going use Qiskit to help us solve the QAOA problem.
Using QAOA with Hamiltonians¶
First, we import Pauli matrix from qiskit.quantum_info and QAOA from
from qiskit.primitives import Sampler
from qiskit.quantum_info import Pauli
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
sampler = Sampler()
H1 = Pauli('Z')^Pauli('Z')
circuit = QAOA(sampler = sampler, optimizer=COBYLA())
#circuit.draw(output = 'mpl')
/var/folders/kz/_mr3r3b55qd2r5hd025yvpfw0000gn/T/ipykernel_95126/2134802490.py:6: DeprecationWarning: The class ``qiskit.primitives.sampler.Sampler`` is deprecated as of qiskit 1.2. It will be removed no earlier than 3 months after the release date. All implementations of the `BaseSamplerV1` interface have been deprecated in favor of their V2 counterparts. The V2 alternative for the `Sampler` class is `StatevectorSampler`. sampler = Sampler()
Here is the guide through the latest qiskit optimization package: https://qiskit-community.github.io/qiskit-optimization/getting_started.html
First, let's install the required qiskit optimization package by running the following code.
#!pip install qiskit-optimization
Requirement already satisfied: qiskit-optimization in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (0.6.1) Requirement already satisfied: qiskit>=0.44 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (1.3.1) Requirement already satisfied: qiskit-algorithms>=0.2.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (0.3.1) Requirement already satisfied: scipy>=1.9.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (1.14.1) Requirement already satisfied: numpy>=1.17 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (2.0.2) Requirement already satisfied: docplex!=2.24.231,>=2.21.207 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (2.29.241) Requirement already satisfied: setuptools>=40.1.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (75.6.0) Requirement already satisfied: networkx>=2.6.3 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit-optimization) (3.4.2) Requirement already satisfied: six in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from docplex!=2.24.231,>=2.21.207->qiskit-optimization) (1.17.0) Requirement already satisfied: rustworkx>=0.15.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (0.15.1) Requirement already satisfied: sympy>=1.3 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (1.12.1) Requirement already satisfied: dill>=0.3 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (0.3.9) Requirement already satisfied: python-dateutil>=2.8.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (2.9.0.post0) Requirement already satisfied: stevedore>=3.0.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (5.4.0) Requirement already satisfied: typing-extensions in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (4.12.2) Requirement already satisfied: symengine<0.14,>=0.11 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from qiskit>=0.44->qiskit-optimization) (0.13.0) Requirement already satisfied: pbr>=2.0.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from stevedore>=3.0.0->qiskit>=0.44->qiskit-optimization) (6.1.0) Requirement already satisfied: mpmath<1.4.0,>=1.1.0 in /opt/anaconda3/envs/pcchouCR97/lib/python3.13/site-packages (from sympy>=1.3->qiskit>=0.44->qiskit-optimization) (1.3.0)
Let's define a quadratic problem with three binary variables. Our quadratic problem has a linear part linear= {'y':-1} and quadratic part quadratic = {('x','y'):2, ('z', 'y'): -4} and linear constraint qp.linear_constraint(linear = {'x':1, 'y':2, 'z':3}, sense= "<=", rhs = 5).
from qiskit_optimization import QuadraticProgram
# Define variables
qp = QuadraticProgram()
qp.binary_var('x')
qp.binary_var('y')
qp.binary_var('z')
qp.minimize(linear= {'y':-1}, quadratic = {('x','y'):2, ('z', 'y'): -4}) # Apply objective function.
qp.linear_constraint(linear = {'x':1, 'y':2, 'z':3}, sense= "<=", rhs = 5) # Apply linear constraint.
# Print problem
print(qp.export_as_lp_string())
\ This file has been generated by DOcplex \ ENCODING=ISO-8859-1 \Problem name: CPLEX Minimize obj: - y + [ 4 x*y - 8 y*z ]/2 Subject To c0: x + 2 y + 3 z <= 5 Bounds 0 <= x <= 1 0 <= y <= 1 0 <= z <= 1 Binaries x y z End
As you can observe from the result, a package from IBM called CPLEX was used to implemented to help with solving optimization problems with classic methods.
Once we have a QuadraticProgram object, we can use the solver provided by Qiskit. Let's import MinimumEigenOptimizer and NumPyMinimumEigensolver from qiskit_optimization.algorithms and qiskit_algorithms, respectively.
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import NumPyMinimumEigensolver
np_solver = NumPyMinimumEigensolver()
np_optimizer = MinimumEigenOptimizer(np_solver)
result = np_optimizer.solve(qp)
print(result)
fval=-5.0, x=0.0, y=1.0, z=1.0, status=SUCCESS
Here, you will see result fval=-5.0, x=0.0, y=1.0, z=1.0, status=SUCCESS. As you can see, we obtain the optimal value of the function as well as its x, y, and z values.
In the same fashion, we can use QAOA to solve the problem with following instructions.
You may face the following issue:
The qiskit.utils.QuantumInstance is a utility class that allows the joint configuration of the circuit transpilation and execution steps, and provides functions at a higher level of abstraction for a more convenient integration with algorithms. These include measurement error mitigation, splitting and combining execution to conform to job limits, and ensuring reliable circuit execution with additional job management tools.
Below are the alternatives:
class qiskit.primitives.Estimator(*, options=None)class qiskit.primitives.Sampler(*, options=None)(Deprecated since version 1.2)
QAOA¶
class QAOA(sampler, optimizer, *, reps=1, initial_state=None, mixer=None, initial_point=None, aggregation=None, callback=None)[source]
Tutorial: https://qiskit-community.github.io/qiskit-algorithms/tutorials/05_qaoa.html
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit_aer.primitives import Estimator, Sampler
from qiskit_algorithms.utils import algorithm_globals
sampler = Sampler()
algorithm_globals.random_seed = 10598
qaoa = QAOA(optimizer=COBYLA(), sampler = sampler, reps = 1)
qaoa_optimizer = MinimumEigenOptimizer(qaoa)
result1 = qaoa_optimizer.solve(qp)
print('result1:',result1)
print('Variable order', [var.name for var in result1.variables])
for s in result1.samples:
print(s)
result1: fval=-5.0, x=0.0, y=1.0, z=1.0, status=SUCCESS Variable order ['x', 'y', 'z'] SolutionSample(x=array([0., 1., 1.]), fval=np.float64(-5.0), probability=0.1025390625, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([0., 1., 0.]), fval=np.float64(-1.0), probability=0.1171875, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([0., 0., 0.]), fval=np.float64(0.0), probability=0.1376953125, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([0., 0., 1.]), fval=np.float64(0.0), probability=0.140625, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([1., 0., 0.]), fval=np.float64(0.0), probability=0.1611328125, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([1., 0., 1.]), fval=np.float64(0.0), probability=0.1982421875, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([1., 1., 0.]), fval=np.float64(1.0), probability=0.0703125, status=<OptimizationResultStatus.SUCCESS: 0>) SolutionSample(x=array([1., 1., 1.]), fval=np.float64(-3.0), probability=0.072265625, status=<OptimizationResultStatus.INFEASIBLE: 2>)
We print out a listing of different solutions that are part of the final, optimal state found by QAOA. Each item of the list includes the assignment the energy or function value, the probability of obtaining the corrsponding basis state when measuring the QAOA state, and whether the solution is feasible or not.
We can also access the full infomation about the QAOA execution by using the following:
print(result1.min_eigen_solver_result)
{ 'aux_operators_evaluated': None,
'best_measurement': { 'bitstring': '000110',
'probability': 0.0380859375,
'state': 6,
'value': np.complex128(-52+0j)},
'cost_function_evals': 30,
'eigenstate': {56: 0.0009765625, 26: 0.0029296875, 19: 0.0048828125, 4: 0.001953125, 54: 0.001953125, 6: 0.009765625, 24: 0.0048828125, 2: 0.0029296875, 27: 0.005859375, 9: 0.0146484375, 11: 0.005859375, 44: 0.0478515625, 39: 0.005859375, 30: 0.0205078125, 36: 0.001953125, 5: 0.03515625, 49: 0.001953125, 29: 0.005859375, 0: 0.005859375, 60: 0.0009765625, 8: 0.021484375, 34: 0.0205078125, 62: 0.013671875, 32: 0.0283203125, 48: 0.044921875, 46: 0.01171875, 53: 0.0009765625, 15: 0.017578125, 28: 0.0322265625, 12: 0.05078125, 40: 0.005859375, 10: 0.0146484375, 18: 0.0185546875, 59: 0.0068359375, 1: 0.0078125, 37: 0.0595703125, 43: 0.0107421875, 50: 0.0419921875, 22: 0.0087890625, 23: 0.01171875, 13: 0.0068359375, 57: 0.033203125, 45: 0.00390625, 51: 0.009765625, 42: 0.015625, 7: 0.033203125, 20: 0.0048828125, 38: 0.0146484375, 16: 0.025390625, 55: 0.00390625, 25: 0.048828125, 61: 0.0205078125, 21: 0.0654296875, 35: 0.0068359375, 3: 0.01953125, 41: 0.0546875, 14: 0.021484375},
'eigenvalue': np.float64(-18.9921875),
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x120e0f6b0>,
'optimal_parameters': { ParameterVectorElement(β[0]): np.float64(-0.3399070054274298),
ParameterVectorElement(γ[0]): np.float64(-4.1087983813741635)},
'optimal_point': array([-0.33990701, -4.10879838]),
'optimal_value': np.float64(-18.9921875),
'optimizer_evals': None,
'optimizer_result': <qiskit_algorithms.optimizers.optimizer.OptimizerResult object at 0x120ee5fd0>,
'optimizer_time': 1.0629889965057373}
We can see that these assignments include slack variables used in transform from constrained to unconstrained problem. Let's obtain the corrsponding QUBO problem with the following code:
from qiskit_optimization.converters import QuadraticProgramToQubo
qp_to_qubo = QuadraticProgramToQubo()
qubo = qp_to_qubo.convert(qp)
print(qubo.export_as_lp_string())
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: CPLEX
Minimize
obj: - 80 x - 161 y - 240 z - 80 c0@int_slack@0 - 160 c0@int_slack@1
- 160 c0@int_slack@2 + [ 16 x^2 + 68 x*y + 96 x*z + 32 x*c0@int_slack@0
+ 64 x*c0@int_slack@1 + 64 x*c0@int_slack@2 + 64 y^2 + 184 y*z
+ 64 y*c0@int_slack@0 + 128 y*c0@int_slack@1 + 128 y*c0@int_slack@2
+ 144 z^2 + 96 z*c0@int_slack@0 + 192 z*c0@int_slack@1
+ 192 z*c0@int_slack@2 + 16 c0@int_slack@0^2
+ 64 c0@int_slack@0*c0@int_slack@1 + 64 c0@int_slack@0*c0@int_slack@2
+ 64 c0@int_slack@1^2 + 128 c0@int_slack@1*c0@int_slack@2
+ 64 c0@int_slack@2^2 ]/2 + 200
Subject To
Bounds
0 <= x <= 1
0 <= y <= 1
0 <= z <= 1
0 <= c0@int_slack@0 <= 1
0 <= c0@int_slack@1 <= 1
0 <= c0@int_slack@2 <= 1
Binaries
x y z c0@int_slack@0 c0@int_slack@1 c0@int_slack@2
End
Congrats! We successfully convert our problem into a QUBO problem where slack variables and penalty terms have been introduced.
It is also worth to notice that qiskit_optimization.converters module in Qiskit also provide InequalityToEquality, IntegerToBinary, and LinearEqualityToPenalty functions. QuadraticProgramToQubo calls them to convert quadratic programs with constraints into QUBO instance, by first introducing slack variables to transform inequalities into equalities, then transforming the integer slack variables into binary once, and finally, replacing the equality constraints with penalty terms.
Run QAOA on the real quantum machine¶
Let's run our problem on the real quantum machine! But, before that, to avoid delay when access quantum device queue many time, we need to get the Hamiltonian first. After having the Hamiltonian, we can use it directly in a QAOA program to submit our request to the real quantum machine just once. Let's see the following code to obtain out Hamiltonian directly!
H1, offset = qubo.to_ising()
print("The Hamiltonian is", H1)
print("The constant term is",offset)
The Hamiltonian is SparsePauliOp(['IIIIZI', 'IIIIIZ', 'IIIZII', 'IIZIII', 'IZIIII', 'ZIIIII', 'IIIIZZ', 'IIIZIZ', 'IIZIIZ', 'IZIIIZ', 'ZIIIIZ', 'IIIZZI', 'IIZIZI', 'IZIIZI', 'ZIIIZI', 'IIZZII', 'IZIZII', 'ZIIZII', 'IZZIII', 'ZIZIII', 'ZZIIII'],
coeffs=[ -7. +0.j, -4.5+0.j, -11. +0.j, -4. +0.j, -8. +0.j, -8. +0.j,
8.5+0.j, 12. +0.j, 4. +0.j, 8. +0.j, 8. +0.j, 23. +0.j,
8. +0.j, 16. +0.j, 16. +0.j, 12. +0.j, 24. +0.j, 24. +0.j,
8. +0.j, 8. +0.j, 16. +0.j])
The constant term is 47.0
Ok, it seems like we can use H1 to solve the problem with QAOA runtime program and even recover the energy by adding the offset term. It seems like a lot of work! Luckly, Qislit provide a simpler way.
The classes VQEClient, QAOAClient, and VQERuntimeResult are removed. Instead, users should migrate their code to use the Qiskit Runtime Primitives with session.
from qiskit_ibm_runtime import QiskitRuntimeService, Session
from qiskit.primitives import Sampler
from qiskit_algorithms import QAOA
service = QiskitRuntimeService(channel="ibm_quantum", #ibm_cloud #https://quantum.ibm.com/
token = '***your own token***')
backend = service.least_busy(operational=True, simulator=False)
#algorithm_globals.random_seed = 10598
session = Session(backend = backend)
#print(session)
sampler = Sampler()
qaoa = QAOA(optimizer=COBYLA(), sampler = sampler, reps = 1)
qaoa_opt = MinimumEigenOptimizer(qaoa)
result2 = qaoa_opt.solve(qp)
print('result2:',result2)
/var/folders/kz/_mr3r3b55qd2r5hd025yvpfw0000gn/T/ipykernel_95126/3017666165.py:6: DeprecationWarning: The class ``qiskit.primitives.sampler.Sampler`` is deprecated as of qiskit 1.2. It will be removed no earlier than 3 months after the release date. All implementations of the `BaseSamplerV1` interface have been deprecated in favor of their V2 counterparts. The V2 alternative for the `Sampler` class is `StatevectorSampler`. sampler = Sampler()
result2: fval=-5.0, x=0.0, y=1.0, z=1.0, status=SUCCESS