import networkx as nx
# Create empty graph
G = nx.Graph()
# Add edges to the graph (also adds nodes)
G.add_edges_from([(1,2),(1,3),(2,4),(3,4),(3,5),(4,5)])
nx.draw(G)
QAOA in partice
As we have seen before, QAOA take the best of both worlds (ideally) considering the ansatz out of the AQC but variationaly training it to obtain the best scheduling function out of the free parameters.
Continuing with previous example….
import numpy as np
from collections import defaultdict
# Initialize our Q matrix
Q = defaultdict(int)
# Update Q matrix for every edge in the graph
for i, j in G.edges:
Q[(i,i)]+= -1
Q[(j,j)]+= -1
Q[(i,j)]+= 2
# Give it also a matrix shape to ease readability
size = G.number_of_nodes()
H_problem = np.zeros((size, size))
for k, v in Q.items():
if not isinstance(k, int):
H_problem[k[0]-1][k[1]-1] = v
H_problemarray([[-2., 2., 2., 0., 0.],
[ 0., -2., 0., 2., 0.],
[ 0., 0., -3., 2., 2.],
[ 0., 0., 0., -3., 2.],
[ 0., 0., 0., 0., -2.]])For this simple example, brute-forcing the solution can be done to have a reference for it.
import numpy as np
def objective_value(x: np.ndarray, w: np.ndarray) -> float:
"""Compute the value of a cut.
Args:
x: Binary string as numpy array.
w: Adjacency matrix.
Returns:
Value of the cut.
"""
return np.dot(np.dot(x.T, w), x)
def bitfield(n: int, L: int) -> list[int]:
""" Get the binary representation"""
result = np.binary_repr(n, L)
return [int(digit) for digit in result] # [2:] to chop off the "0b" part
# use the brute-force way to generate the oracle
L = G.number_of_nodes()
max = 2**L
sol = np.inf
for i in range(max):
cur = bitfield(i, L)
cur_v = objective_value(np.array(cur), H_problem)
if cur_v < sol:
sol = cur_v
print(f'Objective value computed by the brute-force method is {sol}')Objective value computed by the brute-force method is -5.0The main challenge would be understanding if this one is the only solution in our example
x = [0, 1, 1, 0, 0]
print(np.dot(np.dot(np.transpose(x), H_problem), x))-5.0x = [0, 1, 1, 0, 1]
print(np.dot(np.dot(np.transpose(x), H_problem), x))-5.0In fact, for this example we know there are more than one optimal solutions…
from qiskit.visualization import plot_histogram
solution = {
'01101': 0.25,
'10010': 0.25,
'10011': 0.25,
'01100': 0.25,
}
plot_histogram(solution)
We would like to fully understand our solution landscape. Quantum Computing may help us by mapping the solution to a more complex superposed state that reveals all potential solutions to our problem.
We will need to build a Hamiltonian that maps our problem.
from qiskit.quantum_info import SparsePauliOp
H_op = SparsePauliOp.from_list(
[
# h
("ZIIII", -2),
("IZIII", -2),
("IIZII", -3),
("IIIZI", -3),
("IIIIZ", -2),
# j
("ZZIII", 2),
("ZIZII", 2),
("IZIZI", 2),
("IIZIZ", 2),
("IIIZZ", 2),
]
)
print(f"Number of qubits: {H_op.num_qubits}")Number of qubits: 5In this case, we will start with the existing Qiskit functionality and do a drill-down to highlight what is done behind the scenes. Qiskit provides a ready to be used QAOA implementation that will only require the and optimizer and a way to extract information out of our quantum device (either Estimator or Sampler primitive).
from qiskit_algorithms.minimum_eigensolvers import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
optimizer = COBYLA()
sampler = Sampler()
qaoa = QAOA(sampler, optimizer, reps=2)
result = qaoa.compute_minimum_eigenvalue(H_op)/tmp/ipykernel_3017/2458342227.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()print(result){ 'aux_operators_evaluated': None,
'best_measurement': { 'bitstring': '01001',
'probability': np.float64(0.2919347898040691),
'state': 9,
'value': np.complex128(-10+0j)},
'cost_function_evals': np.int64(235),
'eigenstate': {0: np.float64(0.009586873659048), 1: np.float64(0.032574317072889), 2: np.float64(0.002906229790438), 3: np.float64(0.001539056301848), 4: np.float64(0.002906229790438), 5: np.float64(0.001539056301848), 6: np.float64(0.00056519557148), 7: np.float64(9.9947738642e-05), 8: np.float64(0.110558191572831), 9: np.float64(0.291760584421124), 10: np.float64(0.004837359712891), 11: np.float64(0.018755676256582), 12: np.float64(0.014337705300952), 13: np.float64(0.016088090251325), 14: np.float64(0.000691899425862), 15: np.float64(0.001011322269963), 16: np.float64(0.110558191572831), 17: np.float64(0.291760584421124), 18: np.float64(0.014337705300952), 19: np.float64(0.016088090251325), 20: np.float64(0.004837359712891), 21: np.float64(0.018755676256582), 22: np.float64(0.000691899425862), 23: np.float64(0.001011322269963), 24: np.float64(0.009533226937286), 25: np.float64(0.019508719961229), 26: np.float64(1.3202901272e-05), 27: np.float64(0.001494079511016), 28: np.float64(1.3202901272e-05), 29: np.float64(0.001494079511016), 30: np.float64(9.718422913e-06), 31: np.float64(0.000135205204302)},
'eigenvalue': np.float64(-7.702210890981521),
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7fa0b1f05610>,
'optimal_parameters': { ParameterVectorElement(β[0]): np.float64(5.464949090715174),
ParameterVectorElement(β[1]): np.float64(4.960460062899238),
ParameterVectorElement(γ[0]): np.float64(3.9391791330863413),
ParameterVectorElement(γ[1]): np.float64(1.178803173541103)},
'optimal_point': array([5.46494909, 4.96046006, 3.93917913, 1.17880317]),
'optimal_value': np.float64(-7.702210890981521),
'optimizer_evals': None,
'optimizer_result': <qiskit_algorithms.optimizers.optimizer.OptimizerResult object at 0x7fa0c1367080>,
'optimizer_time': 1.5053293704986572}from qiskit.visualization import plot_histogram
counts = {}
for key in result.eigenstate:
key_str = format(key, 'b').zfill(5)
counts[key_str] = result.eigenstate[key]
plot_histogram(counts)
from qiskit.quantum_info import Statevector
from qiskit.result import QuasiDistribution
def sample_most_likely(state_vector: QuasiDistribution | Statevector) -> np.ndarray:
"""Compute the most likely binary string from state vector.
Args:
state_vector: State vector or quasi-distribution.
Returns:
Binary string as an array of ints.
"""
if isinstance(state_vector, QuasiDistribution):
values = list(state_vector.values())
else:
values = state_vector
n = int(np.log2(len(values)))
k = np.argmax(np.abs(values))
x = bitfield(k, n)
x.reverse()
return np.asarray(x)
x = sample_most_likely(result.eigenstate)
print(x)
print(f'Objective value computed by QAOA is {objective_value(x, H_problem)}')[1 0 0 1 0]
Objective value computed by QAOA is -5.0from qiskit_algorithms.optimizers import SLSQP
optimizer = SLSQP()
qaoa = QAOA(sampler, optimizer, reps=4)
result = qaoa.compute_minimum_eigenvalue(H_op)
print(result){ 'aux_operators_evaluated': None,
'best_measurement': { 'bitstring': '01001',
'probability': np.float64(0.1932706955303222),
'state': 9,
'value': np.complex128(-10+0j)},
'cost_function_evals': 285,
'eigenstate': {0: np.float64(0.011305799273325), 1: np.float64(0.008387295761362), 2: np.float64(0.024275579531199), 3: np.float64(0.006249560281358), 4: np.float64(0.024275579531199), 5: np.float64(0.006249560281358), 6: np.float64(0.033057833032238), 7: np.float64(0.000266368470969), 8: np.float64(0.063503350506968), 9: np.float64(0.192309032378154), 10: np.float64(0.006308056514893), 11: np.float64(0.000851223959722), 12: np.float64(0.101498567702871), 13: np.float64(0.002600321525932), 14: np.float64(0.002606855969397), 15: np.float64(0.001962437215773), 16: np.float64(0.063503350506968), 17: np.float64(0.192309032378154), 18: np.float64(0.101498567702871), 19: np.float64(0.002600321525932), 20: np.float64(0.006308056514893), 21: np.float64(0.000851223959722), 22: np.float64(0.002606855969397), 23: np.float64(0.001962437215773), 24: np.float64(0.029283385169574), 25: np.float64(0.107498197914821), 26: np.float64(0.002675540405354), 27: np.float64(3.852936896e-05), 28: np.float64(0.002675540405354), 29: np.float64(3.852936896e-05), 30: np.float64(4.550530625e-06), 31: np.float64(0.000438459125925)},
'eigenvalue': np.float64(-7.353101752321757),
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7fa0ac7bf110>,
'optimal_parameters': { ParameterVectorElement(γ[1]): np.float64(4.548828250065513),
ParameterVectorElement(γ[2]): np.float64(-2.6864468837121773),
ParameterVectorElement(γ[3]): np.float64(-3.6679971855296025),
ParameterVectorElement(γ[0]): np.float64(-6.140006111232352),
ParameterVectorElement(β[3]): np.float64(2.1446195228385796),
ParameterVectorElement(β[2]): np.float64(3.222554048781957),
ParameterVectorElement(β[0]): np.float64(2.7464433448437005),
ParameterVectorElement(β[1]): np.float64(4.292780823213835)},
'optimal_point': array([ 2.74644334, 4.29278082, 3.22255405, 2.14461952, -6.14000611,
4.54882825, -2.68644688, -3.66799719]),
'optimal_value': np.float64(-7.353101752321757),
'optimizer_evals': None,
'optimizer_result': <qiskit_algorithms.optimizers.optimizer.OptimizerResult object at 0x7fa0cff72900>,
'optimizer_time': 2.273702621459961}counts = {}
for key in result.eigenstate:
key_str = format(key, 'b').zfill(5)
counts[key_str] = result.eigenstate[key]
plot_histogram(counts)
x = sample_most_likely(result.eigenstate)
print(x)
print(f'Objective value computed by QAOA is {objective_value(x, H_problem)}')[1 0 0 1 0]
Objective value computed by QAOA is -5.0optimizer = SLSQP(
maxiter=100000
)
qaoa = QAOA(sampler, optimizer, reps=1)
result = qaoa.compute_minimum_eigenvalue(H_op)
print(result){ 'aux_operators_evaluated': None,
'best_measurement': { 'bitstring': '01001',
'probability': np.float64(0.1481940026773182),
'state': 9,
'value': np.complex128(-10+0j)},
'cost_function_evals': 38,
'eigenstate': {0: np.float64(0.049683862130195), 1: np.float64(0.074154922931544), 2: np.float64(0.008731076933889), 3: np.float64(0.003870276043694), 4: np.float64(0.008731076933889), 5: np.float64(0.003870276043694), 6: np.float64(0.00128492491713), 7: np.float64(0.000170123829585), 8: np.float64(0.091089635773707), 9: np.float64(0.145557663369752), 10: np.float64(0.004844518650464), 11: np.float64(0.000784455359524), 12: np.float64(0.016048795577276), 13: np.float64(0.007854097184562), 14: np.float64(0.000653593486794), 15: np.float64(0.000103499896109), 16: np.float64(0.091089635773707), 17: np.float64(0.145557663369752), 18: np.float64(0.016048795577276), 19: np.float64(0.007854097184562), 20: np.float64(0.004844518650464), 21: np.float64(0.000784455359524), 22: np.float64(0.000653593486794), 23: np.float64(0.000103499896109), 24: np.float64(0.111017103835445), 25: np.float64(0.188436458771374), 26: np.float64(0.006211782518087), 27: np.float64(0.001688564507361), 28: np.float64(0.006211782518087), 29: np.float64(0.001688564507361), 30: np.float64(0.000288102298604), 31: np.float64(8.8582683684e-05)},
'eigenvalue': np.float64(-6.210057646396847),
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7fa0b0ca2f90>,
'optimal_parameters': { ParameterVectorElement(γ[0]): np.float64(7.0237858493936915),
ParameterVectorElement(β[0]): np.float64(0.9327190506300823)},
'optimal_point': array([0.93271905, 7.02378585]),
'optimal_value': np.float64(-6.210057646396847),
'optimizer_evals': None,
'optimizer_result': <qiskit_algorithms.optimizers.optimizer.OptimizerResult object at 0x7fa0b0d4c800>,
'optimizer_time': 0.1328577995300293}counts = {}
for key in result.eigenstate:
key_str = format(key, 'b').zfill(5)
counts[key_str] = result.eigenstate[key]
plot_histogram(counts)
x = sample_most_likely(result.eigenstate)
print(x)
print(f'Objective value computed by QAOA is {objective_value(x, H_problem)}')[1 0 0 1 1]
Objective value computed by QAOA is -5.0Let’s check what is going on behind it. We can check the ansatz that QAOA creates.
qaoa.ansatz.draw('mpl', fold=150)
qaoa.ansatz.decompose().draw('mpl', fold=150)
qaoa.ansatz.decompose(reps=2).draw('mpl', fold=150)
As we can see the ansatz is build following the Adiabatic Quantum Computing scheme where a block of \(X\) rotations is applied, and alternated with the problem hamiltonian (\(Z + ZZ\) interactions) as many times as repetitions are performed.
qaoa = QAOA(sampler, optimizer, reps=4)
result = qaoa.compute_minimum_eigenvalue(H_op)
print(result){ 'aux_operators_evaluated': None,
'best_measurement': { 'bitstring': '01001',
'probability': np.float64(0.2016985082946548),
'state': 9,
'value': np.complex128(-10+0j)},
'cost_function_evals': 384,
'eigenstate': {0: np.float64(0.007706890589662), 1: np.float64(0.068172655547573), 2: np.float64(0.029654423907845), 3: np.float64(0.002361049740202), 4: np.float64(0.029654423907845), 5: np.float64(0.002361049740202), 6: np.float64(0.073738195216247), 7: np.float64(0.002590380910088), 8: np.float64(0.009557277637467), 9: np.float64(0.183315435377866), 10: np.float64(0.010115526796152), 11: np.float64(0.003120300666106), 12: np.float64(0.063221065422187), 13: np.float64(0.006017369449093), 14: np.float64(0.002540148318749), 15: np.float64(0.007051902511363), 16: np.float64(0.009557277637467), 17: np.float64(0.183315435377866), 18: np.float64(0.063221065422187), 19: np.float64(0.006017369449093), 20: np.float64(0.010115526796152), 21: np.float64(0.003120300666106), 22: np.float64(0.002540148318749), 23: np.float64(0.007051902511363), 24: np.float64(0.083105955068743), 25: np.float64(0.103520231032441), 26: np.float64(0.004112007698028), 27: np.float64(0.000456117876314), 28: np.float64(0.004112007698028), 29: np.float64(0.000456117876314), 30: np.float64(0.0142960431147), 31: np.float64(0.003824397717799)},
'eigenvalue': np.float64(-6.293446338626337),
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7fa0a4429af0>,
'optimal_parameters': { ParameterVectorElement(β[1]): np.float64(-0.8481463411342136),
ParameterVectorElement(β[0]): np.float64(-9.961772667943281),
ParameterVectorElement(β[3]): np.float64(5.537169496503184),
ParameterVectorElement(γ[1]): np.float64(2.4047469775752717),
ParameterVectorElement(γ[2]): np.float64(6.719784554662614),
ParameterVectorElement(γ[3]): np.float64(-3.6116397606401276),
ParameterVectorElement(γ[0]): np.float64(-0.3664281133719904),
ParameterVectorElement(β[2]): np.float64(-8.962238627520238)},
'optimal_point': array([-9.96177267, -0.84814634, -8.96223863, 5.5371695 , -0.36642811,
2.40474698, 6.71978455, -3.61163976]),
'optimal_value': np.float64(-6.293446338626337),
'optimizer_evals': None,
'optimizer_result': <qiskit_algorithms.optimizers.optimizer.OptimizerResult object at 0x7fa0a446e810>,
'optimizer_time': 3.110597848892212}qaoa.ansatz.draw('mpl', fold=150)
qaoa.ansatz.decompose().decompose().draw('mpl', fold=150)
Obtained solution distribution really points out to the solutions that provide solutions to our problem, superposed in a state that is the one minimizing the outcome of the circuit with respect to the expectation value of the hamitlonian.
counts = {}
for key in result.eigenstate:
key_str = format(key, 'b').zfill(5)
counts[key_str] = result.eigenstate[key]
plot_histogram(counts)
Therefore, we are able to encode our problem into an Ising Hamiltonian, so that we can easily design a circuit to obtain the results from the potential solution landscape. We may need to set up the repetitions needed in order to map the optimal parameters to our variationaly trained scheduling function.
qaoa.ansatz.depth()10Small problems will render shallow circuits, but we need to be aware that introducing many repetitions of those blocks or layers, it will require more operations to be squeezed in, the lifespan of the cubits will need to be longer and the more operations we add, more noise it will enter into our final state.