QAOA
This notebook demonstrates the Classiq performance re. the Quantum Approximate Optimization Algorithm (QAOA), focusing on the Max Clique problem.
1. Calling the Built-in QAOA
This section calls the built-in QAOA of Classiq, constructing the corresponding quantum model from a combinatorial optimization Pyomo model.
1.1. Generating a Pyomo Model
import time
import networkx as nx
import numpy as np
import pyomo.environ as pyo
np.random.seed(2)
def define_max_clique_model(graph):
model = pyo.ConcreteModel()
# each x_i states if node i belongs to the cliques
model.x = pyo.Var(graph.nodes, domain=pyo.Binary)
x_variables = np.array(list(model.x.values()))
# define the complement adjacency matrix as the matrix where 1 exists for each non-existing edge
adjacency_matrix = nx.convert_matrix.to_numpy_array(graph, nonedge=0)
complement_adjacency_matrix = (
1
- nx.convert_matrix.to_numpy_array(graph, nonedge=0)
- np.identity(len(model.x))
)
# constraint that 2 nodes without an edge in the graph cannot be chosen together
model.clique_constraint = pyo.Constraint(
expr=x_variables @ complement_adjacency_matrix @ x_variables == 0
)
# maximize the number of nodes in the chosen clique
model.value = pyo.Objective(expr=sum(x_variables), sense=pyo.maximize)
return model
Setting a specific problem and some hyperparameters.
QAOA_NUM_LAYERS = 10
NUM_SHOTS = 1e4
NUM_QUBITS = 7
graph = nx.erdos_renyi_graph(NUM_QUBITS, 0.6, seed=79)
max_clique_model = define_max_clique_model(graph)
# transpilation_options = {"classiq": "custom", "qiskit": 3}
transpilation_options = {"classiq": "auto optimize", "qiskit": 1}
1.2 Constructing, Synthesizing, and Running a QAOA Model
from classiq import (
CustomHardwareSettings,
Preferences,
QuantumProgram,
construct_combinatorial_optimization_model,
execute,
set_execution_preferences,
set_preferences,
show,
synthesize,
)
from classiq.applications.combinatorial_optimization import OptimizerConfig, QAOAConfig
from classiq.execution import ExecutionPreferences
qaoa_config = QAOAConfig(num_layers=QAOA_NUM_LAYERS)
optimizer_config = OptimizerConfig(max_iteration=400, alpha_cvar=1)
qmod = construct_combinatorial_optimization_model(
pyo_model=max_clique_model,
qaoa_config=qaoa_config,
optimizer_config=optimizer_config,
)
execution_preferences = ExecutionPreferences(num_shots=NUM_SHOTS)
preferences = Preferences(
custom_hardware_settings=CustomHardwareSettings(basis_gates=["cx", "u"]),
transpilation_option=transpilation_options["classiq"],
)
qmod = set_execution_preferences(qmod, execution_preferences)
qmod = set_preferences(qmod, preferences=preferences)
qprog = synthesize(qmod)
res = execute(qprog).result()
circuit = QuantumProgram.from_qprog(qprog)
depth_classiq = circuit.transpiled_circuit.depth
cx_counts_classiq = circuit.transpiled_circuit.count_ops["cx"]
classiq_solving_time = res[0].value.time
interations_classiq = [
intermediate_result.iteration_number
for intermediate_result in res[0].value.intermediate_results
]
results_classiq = [
-intermediate_result.mean_all_solutions
for intermediate_result in res[0].value.intermediate_results
]
2. Comparing to Qiskit
We use qiskit version 0.39.4.
Qiskit has no module in which to specify a generic optimization problem; therefore, you have to do the preprocessing and post-processing yourself. Retrieve the Hamiltonian that enters into the VQE.
from typing import List
from classiq import Pauli as ClassiqPauli
def get_classiq_hamiltonian(execution_result) -> List[List[str]]:
hamiltonian_result = execution_result[1].value
parsed_pauli_list = list()
for pauli_term in hamiltonian_result:
pauli_str = "".join(ClassiqPauli(pauli).name for pauli in pauli_term["pauli"])
coefficient = pauli_term["coefficient"]
parsed_pauli_list.append([pauli_str, coefficient])
return parsed_pauli_list
Define a function for running QAOA on Qiskit and returning the results.
from qiskit import QuantumCircuit, transpile
from qiskit.algorithms.minimum_eigensolvers import QAOA
from qiskit.algorithms.optimizers import COBYLA
from qiskit.primitives import Estimator, Sampler
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.result import QuasiDistribution
def qiskit_qaoa(hamiltonian, qaoa_num_layers, num_qubits):
"""
Gets a Hamiltonian for QAOA and num of quantum layers, returning the most probable solution and its corresponding cost
as well as intermediate results
"""
counts = []
values = []
def store_intermediate_result(eval_count, parameters, mean, std):
# callable to store results
counts.append(eval_count)
values.append(mean)
def objective_value(x, hamiltonian):
# get objective value for a given computational basis state
qc = QuantumCircuit(num_qubits)
for k in range(num_qubits):
if x[k] == 1:
qc.x(k)
estimated_cost = estimator.run(qc, hamiltonian).result().values[0]
return -estimated_cost
def bitfield(n: int, L: int) -> list[int]:
# binary representation as list
result = np.binary_repr(n, L)
return [int(digit) for digit in result] # [2:] to chop off the "0b" part
def sample_most_likely(state_vector: QuasiDistribution | Statevector) -> np.ndarray:
"""Compute the most likely binary string from the 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)
estimator = Estimator(options={"shots": int(NUM_SHOTS)})
sampler = Sampler()
optimizer = COBYLA()
qaoa = QAOA(
sampler, optimizer, reps=qaoa_num_layers, callback=store_intermediate_result
)
start = time.time()
result = qaoa.compute_minimum_eigenvalue(hamiltonian)
solving_time = time.time() - start
x = sample_most_likely(result.eigenstate)
transpiled_circuit = transpile(
result.optimal_circuit,
basis_gates=["u", "cx"],
optimization_level=transpilation_options["qiskit"],
)
return (
transpiled_circuit.depth(),
transpiled_circuit.count_ops()["cx"],
solving_time,
x,
objective_value(x, hamiltonian),
counts,
values,
)
The same QAOA with Qiskit
pauli_list = get_classiq_hamiltonian(res)
hamiltonian = SparsePauliOp.from_list(pauli_list)
(
depth_qiskit,
cx_counts_qiskit,
time_qiskit,
most_probable_state,
cost,
iterations_qiskit,
results_qiskit,
) = qiskit_qaoa(hamiltonian, QAOA_NUM_LAYERS, NUM_QUBITS)
3. Plotting the Data
import matplotlib.pyplot as plt
classiq_color = "#F43764"
qiskit_color = "#6FA4FF"
plt.rcParams["font.family"] = "serif"
plt.rc("savefig", dpi=300)
plt.rcParams["axes.linewidth"] = 1
plt.rcParams["xtick.major.size"] = 5
plt.rcParams["xtick.minor.size"] = 5
plt.rcParams["ytick.major.size"] = 5
plt.rcParams["ytick.minor.size"] = 5
qiskit_label = f"qiskit: depth={depth_qiskit} \n cx-counts={cx_counts_qiskit}"
classiq_label = f"classiq: depth={depth_classiq},\n cx-counts={cx_counts_classiq}"
max_classiq_iter = len(interations_classiq)
plt.plot(
iterations_qiskit[0:max_classiq_iter],
np.real(results_qiskit[0:max_classiq_iter]),
"-",
label=qiskit_label,
linewidth=1.5,
color=qiskit_color,
)
plt.plot(
interations_classiq[0:max_classiq_iter],
results_classiq[0:max_classiq_iter],
"-",
label=classiq_label,
linewidth=1.5,
color=classiq_color,
)
# plt.ylim(0,90)
plt.ylabel("Energy", fontsize=16)
plt.xlabel("Iteration", fontsize=16)
plt.yticks(fontsize=16)
plt.xticks(fontsize=16)
plt.legend(loc="upper right", fontsize=16)
<matplotlib.legend.Legend at 0x7feedc946c50>
