{
"cells": [
{
"cell_type": "markdown",
"id": "bc51e7bf-e582-49ba-93f8-035624d56ccf",
"metadata": {},
"source": [
"---\n",
"title: Quantum approximate optimization algorithm\n",
"description: Solve max-cut using QAOA with a Qiskit pattern at utility scale.\n",
"---\n",
"\n",
"{/* cspell:ignore frameon popcount fval */}\n",
"\n",
"# Quantum approximate optimization algorithm\n",
"\n",
"*Usage estimate: 22 minutes on a Heron r3 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "learning-outcomes-prereqs",
"metadata": {},
"source": [
"## Learning outcomes\n",
"\n",
"After completing this tutorial, you can expect to understand the following information:\n",
"\n",
"* How to map a classical combinatorial optimization problem (max-cut) to a quantum Hamiltonian\n",
"* How to implement and run the Quantum Approximate Optimization Algorithm (QAOA) using IBM Quantum Compute Service sessions\n",
"* How to scale a QAOA workflow from a small simulator example to utility-scale hardware execution\n",
"\n",
"## Prerequisites\n",
"\n",
"It is recommended that you familiarize yourself with these topics:\n",
"\n",
"* [Basics of quantum circuits](/learning/courses/basics-of-quantum-information)\n",
"* [Variational algorithms](/learning/courses/variational-algorithm-design)\n",
"* [QAOA in depth](/learning/courses/quantum-computing-in-practice/utility-scale-qaoa) — for a comprehensive treatment of the QAOA algorithm and applying it at utility scale\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "de201dbb",
"metadata": {},
"source": [
"## Background\n",
"\n",
"The **Quantum Approximate Optimization Algorithm (QAOA)** is a hybrid quantum-classical iterative method for solving combinatorial optimization problems. In this tutorial, you will use QAOA to solve the **maximum-cut (max-cut)** problem — an NP-hard optimization problem with applications in clustering, network science, and statistical physics. Given a graph of nodes connected by edges, the goal is to partition the nodes into two sets such that the number of edges crossing the partition is maximized.\n",
"\n",
"\n",
"\n",
"### From classical optimization to quantum circuits\n",
"\n",
"Max-cut can be expressed as a classical binary optimization problem. Each node is assigned a binary variable $x_i \\in \\{0, 1\\}$ indicating which set it belongs to. The objective is to maximize the number of edges where the endpoints are in different sets:\n",
"\n",
"$$\n",
"\\max_{x \\in \\{0,1\\}^n} \\sum_{(i,j)} x_i + x_j - 2x_ix_j.\n",
"$$\n",
"\n",
"This is equivalently a **Quadratic Unconstrained Binary Optimization (QUBO)** problem of the form $\\min_x\\, x^T Q x$. Through a standard variable substitution ($x_i \\to (1 - Z_i)/2$), the QUBO can be rewritten as a **cost Hamiltonian** whose ground state encodes the optimal solution. In general, this Hamiltonian has both quadratic and linear terms:\n",
"\n",
"$$\n",
"H_C = \\sum_{ij} Q_{ij} \\, Z_i Z_j + \\sum_i b_i \\, Z_i.\n",
"$$\n",
"\n",
"For the unweighted max-cut problem considered here, the linear coefficients vanish ($b_i = 0$) and $Q_{ij} = 1$ for each edge, leaving the simpler form $H_C = \\sum_{(i,j) \\in E} Z_i Z_j$ that you will build in code below. The more general form above is what you would need to adapt this workflow to weighted graphs or other QUBO-expressible problems.\n",
"\n",
"### How QAOA works\n",
"\n",
"QAOA prepares candidate solutions by applying alternating layers of two operators to an initial superposition state $H^{\\otimes n}|0\\rangle$: the **cost operator** $e^{-i\\gamma_k H_C}$ and a **mixer operator** $e^{-i\\beta_k H_m}$. The angles $\\gamma_k$ and $\\beta_k$ are optimized in a classical feedback loop; the quantum computer evaluates the cost function, and a classical optimizer updates the parameters until convergence. This iterative loop runs within a Quantum Compute **session**, which keeps the quantum device reserved across iterations for lower latency.\n",
"\n",
"\n",
"\n",
"For a deeper treatment of QAOA theory, including the full QUBO-to-Hamiltonian derivation, see the [QAOA course module](/learning/courses/utility-scale-quantum-computing/variational-quantum-algorithms#1-introduction).\n",
"\n",
"In this tutorial you will first solve max-cut on a small five-node graph, then scale the same workflow to a 100-node utility-scale problem on real hardware. *Note on plan access:* This tutorial uses Quantum Compute [sessions](/docs/guides/execution-modes#session-mode), which are only available on the Premium Plan. If you are on the Open Plan, you cannot run this tutorial as written; instead, you will need to swap `Session` for [job mode](/docs/guides/execution-modes#job-mode) (as in, submit each iteration as an independent job rather than wrapping the optimization loop in with `Session(...)`). The workflow still runs, but each iteration incurs the full queue latency rather than reusing a reserved device. See [Overview of available plans](/docs/guides/plans-overview) for more information.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "381800e5",
"metadata": {},
"source": [
"## Requirements\n",
"\n",
"Before starting this tutorial, be sure you have the following installed:\n",
"\n",
"* Qiskit SDK v2.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
"* Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n",
"\n",
"In addition, you will need access to an instance on [IBM Quantum® Platform](/docs/guides/cloud-setup).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f5307376",
"metadata": {},
"source": [
"## Setup\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "37b3acfc",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import rustworkx as rx\n",
"from rustworkx.visualization import mpl_draw as draw_graph\n",
"import numpy as np\n",
"from scipy.optimize import minimize\n",
"from collections import defaultdict\n",
"from typing import Sequence\n",
"\n",
"\n",
"from qiskit.quantum_info import SparsePauliOp\n",
"from qiskit.circuit.library import QAOAAnsatz\n",
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"\n",
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator\n",
"from qiskit_ibm_runtime import SamplerV2 as Sampler"
]
},
{
"cell_type": "markdown",
"id": "68fd0b4f-baa4-45dc-9f4c-d9cdff01a651",
"metadata": {},
"source": [
"## Small-scale example\n",
"\n",
"This section walks through each step of the QAOA workflow on a small five-node max-cut instance. Despite the \"small-scale\" label, this example still runs on real IBM Quantum hardware — the code selects a backend with 127 or more qubits and executes the circuit there.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "0cbca6cb",
"metadata": {},
"source": [
"Initialize your problem by creating a graph with $n=5$ nodes.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "6ced6bea",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n_small = 5\n",
"\n",
"graph = rx.PyGraph()\n",
"graph.add_nodes_from(np.arange(0, n_small, 1))\n",
"edge_list = [\n",
" (0, 1, 1.0),\n",
" (0, 2, 1.0),\n",
" (0, 4, 1.0),\n",
" (1, 2, 1.0),\n",
" (2, 3, 1.0),\n",
" (3, 4, 1.0),\n",
"]\n",
"graph.add_edges_from(edge_list)\n",
"draw_graph(graph, node_size=600, with_labels=True)"
]
},
{
"cell_type": "markdown",
"id": "a06e4386-d7bd-4914-9baa-36a5cc60e3ab",
"metadata": {},
"source": [
"### Step 1: Map classical inputs to a quantum problem\n",
"\n",
"Map the classical graph into quantum **circuits** and **operators**. As described in the [Background](#background), for unweighted max-cut the cost Hamiltonian reduces to $H_C = \\sum_{(i,j) \\in E} Z_i Z_j$, and QAOA uses a parametrized ansatz circuit to prepare candidate ground states of $H_C$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a5b9e551-38a1-4543-b9f1-caaefb0ef3a9",
"metadata": {},
"source": [
"#### Build the cost Hamiltonian\n",
"\n",
"Convert the graph edges into Pauli $Z_iZ_j$ terms to construct $H_C$ (see [Background](#background) for the derivation).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "52d1ba92",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'ZIIIZ', 'IIZZI', 'IZZII', 'ZZIII'],\n",
" coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n"
]
}
],
"source": [
"def build_max_cut_paulis(\n",
" graph: rx.PyGraph,\n",
") -> list[tuple[str, list[int], float]]:\n",
" \"\"\"Convert graph edges to a list of ZZ Pauli terms.\n",
"\n",
" The returned list is in the sparse format expected by\n",
" ``SparsePauliOp.from_sparse_list``: each element is\n",
" ``(pauli_string, qubit_indices, coefficient)``.\n",
" \"\"\"\n",
" pauli_list = []\n",
" for edge in list(graph.edge_list()):\n",
" weight = graph.get_edge_data(edge[0], edge[1])\n",
" pauli_list.append((\"ZZ\", [edge[0], edge[1]], weight))\n",
" return pauli_list\n",
"\n",
"\n",
"max_cut_paulis = build_max_cut_paulis(graph)\n",
"cost_hamiltonian = SparsePauliOp.from_sparse_list(max_cut_paulis, n_small)\n",
"print(\"Cost Function Hamiltonian:\", cost_hamiltonian)"
]
},
{
"cell_type": "markdown",
"id": "33f71b0d-4a2a-4082-8c1a-ce9d2b769048",
"metadata": {},
"source": [
"#### Build the QAOA ansatz circuit\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "00431c46-30c2-40f9-99df-40baf8da98f6",
"metadata": {},
"source": [
"Use `QAOAAnsatz` to construct the parametrized QAOA circuit from the cost Hamiltonian. Here we use `reps=2` (two QAOA layers, four parameters: $\\beta_0, \\beta_1, \\gamma_0, \\gamma_1$).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "7bd8c6d4-f40f-4a11-a440-0b26d9021b53",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=2)\n",
"circuit.measure_all()\n",
"\n",
"circuit.draw(\"mpl\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "315c495a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(β[1]), ParameterVectorElement(γ[0]), ParameterVectorElement(γ[1])])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"circuit.parameters"
]
},
{
"cell_type": "markdown",
"id": "82f70daa-ff68-447a-8064-8b7df7a646cf",
"metadata": {},
"source": [
"### Step 2: Optimize problem for quantum hardware execution\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c08be444-e3ed-4178-a10b-414069b1b411",
"metadata": {},
"source": [
"Transpile the abstract circuit into hardware-native instructions. This step handles qubit mapping, gate decomposition, routing, and error suppression. See the transpilation [documentation](/docs/guides/transpile) for more information.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "3f28a422-805c-4d3d-b5f6-62539e9133bd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"service = QiskitRuntimeService()\n",
"backend = service.least_busy(\n",
" operational=True, simulator=False, min_num_qubits=127\n",
")\n",
"print(backend)\n",
"\n",
"# Create pass manager for transpilation. Level 3 is the most aggressive\n",
"# preset: slower to transpile, but produces shorter circuits that are\n",
"# more robust to hardware noise.\n",
"pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n",
"\n",
"candidate_circuit = pm.run(circuit)\n",
"candidate_circuit.draw(\"mpl\", fold=False, idle_wires=False)"
]
},
{
"cell_type": "markdown",
"id": "4e75cad7-f599-4937-b5fe-f4d01f53423c",
"metadata": {},
"source": [
"### Step 3: Execute using Qiskit primitives\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9b99ce67-f121-4244-b62a-536be38fea86",
"metadata": {},
"source": [
"The QAOA optimization loop runs inside a Quantum Compute [session](/docs/guides/execution-modes) to keep the device reserved across iterations. An Estimator evaluates $\\langle H_C \\rangle$ at each step, and a classical optimizer (COBYLA) updates the parameters until convergence.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "00b2b0f1-9bad-4ad3-b93e-5cbf40395dbf",
"metadata": {},
"source": [
"Define initial parameters and run the optimization loop:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "afa5747f-44dc-4e41-a875-7b6f896f13e2",
"metadata": {},
"outputs": [],
"source": [
"# QAOA doesn't prescribe principled default angles — any bounded choice\n",
"# works as a warm start for problems this small. beta and gamma are\n",
"# periodic (beta in [0, pi] and gamma in [0, 2*pi] modulo the underlying\n",
"# Pauli-rotation periods), and pi/2 and pi are just midpoints of those\n",
"# ranges. For harder problems you would typically warm start from known\n",
"# good angles or transfer parameters from smaller instances.\n",
"initial_gamma = np.pi\n",
"initial_beta = np.pi / 2\n",
"init_params = [initial_beta, initial_beta, initial_gamma, initial_gamma]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "3e64a862",
"metadata": {},
"outputs": [],
"source": [
"def cost_func_estimator(params, ansatz, hamiltonian, estimator):\n",
" # transform the observable defined on virtual qubits to\n",
" # an observable defined on all physical qubits\n",
" isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)\n",
"\n",
" pub = (ansatz, isa_hamiltonian, params)\n",
" job = estimator.run([pub])\n",
"\n",
" results = job.result()[0]\n",
" cost = results.data.evs\n",
"\n",
" objective_func_vals.append(cost)\n",
"\n",
" return cost"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "2df241a9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" message: Return from COBYLA because the trust region radius reaches its lower bound.\n",
" success: True\n",
" status: 0\n",
" fun: -2.0402211719947774\n",
" x: [ 3.041e+00 1.212e+00 2.081e+00 4.471e+00]\n",
" nfev: 36\n",
" maxcv: 0.0\n"
]
}
],
"source": [
"objective_func_vals = [] # Global variable\n",
"with Session(backend=backend) as session:\n",
" # If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`\n",
" estimator = Estimator(mode=session)\n",
" estimator.options.default_shots = 1000\n",
"\n",
" # Set simple error suppression/mitigation options\n",
" estimator.options.dynamical_decoupling.enable = True\n",
" estimator.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
" estimator.options.twirling.enable_gates = True\n",
" estimator.options.twirling.num_randomizations = \"auto\"\n",
" estimator.options.environment.job_tags = [\"TUT_QAOA\"]\n",
"\n",
" result = minimize(\n",
" cost_func_estimator,\n",
" init_params,\n",
" args=(candidate_circuit, cost_hamiltonian, estimator),\n",
" method=\"COBYLA\",\n",
" tol=1e-2,\n",
" )\n",
" print(result)"
]
},
{
"cell_type": "markdown",
"id": "01d6b81c",
"metadata": {},
"source": [
"The optimizer was able to reduce the cost and find better parameters for the circuit.\n",
"\n",
"A smoothly decreasing curve that plateaus is the signature of convergence. A noisy, non-monotonic curve usually indicates that something upstream needs attention; common causes are too few shots per evaluation (high estimator variance), poor initial parameters, or a circuit whose depth is dominated by hardware noise. COBYLA is derivative-free and fairly robust to moderate noise, but when the noise swamps the actual cost improvements per step, its linear-approximation model can no longer tell real descent from random jitter and the optimizer wanders.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e14ecc92",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(12, 6))\n",
"plt.plot(objective_func_vals)\n",
"plt.xlabel(\"Iteration\")\n",
"plt.ylabel(\"Cost\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "1f9c8a9c",
"metadata": {},
"source": [
"Assign the optimized parameters and sample the final distribution using the Sampler primitive.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "2989e76e-4296-4dd8-b065-2b8fced064cf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"optimized_circuit = candidate_circuit.assign_parameters(result.x)\n",
"optimized_circuit.draw(\"mpl\", fold=False, idle_wires=False)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "d8f0e302",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{18: 0.039, 5: 0.0665, 20: 0.0973, 29: 0.0063, 9: 0.0899, 13: 0.0379, 2: 0.0047, 1: 0.0153, 11: 0.0932, 14: 0.0327, 12: 0.0314, 25: 0.0193, 21: 0.0398, 6: 0.0224, 4: 0.0197, 10: 0.0387, 3: 0.0181, 26: 0.07, 17: 0.0327, 19: 0.0332, 22: 0.0914, 24: 0.007, 0: 0.0033, 8: 0.0066, 30: 0.0158, 28: 0.0169, 27: 0.0222, 16: 0.0073, 7: 0.0057, 23: 0.0062, 15: 0.0054, 31: 0.0041}\n"
]
}
],
"source": [
"# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`\n",
"sampler = Sampler(mode=backend)\n",
"sampler.options.default_shots = 10000\n",
"\n",
"# Set simple error suppression/mitigation options\n",
"sampler.options.dynamical_decoupling.enable = True\n",
"sampler.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
"sampler.options.twirling.enable_gates = True\n",
"sampler.options.twirling.num_randomizations = \"auto\"\n",
"\n",
"sampler.options.environment.job_tags = [\"TUT_QAOA\"]\n",
"\n",
"pub = (optimized_circuit,)\n",
"job = sampler.run([pub], shots=int(1e4))\n",
"counts_int = job.result()[0].data.meas.get_int_counts()\n",
"counts_bin = job.result()[0].data.meas.get_counts()\n",
"shots = sum(counts_int.values())\n",
"final_distribution_int = {key: val / shots for key, val in counts_int.items()}\n",
"final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}\n",
"print(final_distribution_int)"
]
},
{
"cell_type": "markdown",
"id": "dace5fed-5555-4f1c-9109-7f5a31832d04",
"metadata": {},
"source": [
"### Step 4: Post-process and return result in desired classical format\n",
"\n",
"Extract the most likely bitstring from the sampled distribution. This represents the best cut found by QAOA.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "d4f7fc70-883f-4b6b-8e92-2fc4afbbea46",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result bitstring: [0, 0, 1, 0, 1]\n"
]
}
],
"source": [
"# auxiliary functions to sample most likely bitstring\n",
"def to_bitstring(integer, num_bits):\n",
" result = np.binary_repr(integer, width=num_bits)\n",
" return [int(digit) for digit in result]\n",
"\n",
"\n",
"keys = list(final_distribution_int.keys())\n",
"values = list(final_distribution_int.values())\n",
"most_likely = keys[np.argmax(np.abs(values))]\n",
"most_likely_bitstring = to_bitstring(most_likely, len(graph))\n",
"most_likely_bitstring.reverse()\n",
"\n",
"print(\"Result bitstring:\", most_likely_bitstring)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "650875e9-adbc-43bd-9505-556be2566278",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.rcParams.update({\"font.size\": 10})\n",
"final_bits = final_distribution_bin\n",
"values = np.abs(list(final_bits.values()))\n",
"top_4_values = sorted(values, reverse=True)[:4]\n",
"positions = []\n",
"for value in top_4_values:\n",
" positions.append(np.where(values == value)[0])\n",
"fig = plt.figure(figsize=(11, 6))\n",
"ax = fig.add_subplot(1, 1, 1)\n",
"plt.xticks(rotation=45)\n",
"plt.title(\"Result Distribution\")\n",
"plt.xlabel(\"Bitstrings (reversed)\")\n",
"plt.ylabel(\"Probability\")\n",
"ax.bar(list(final_bits.keys()), list(final_bits.values()), color=\"tab:grey\")\n",
"for p in positions:\n",
" ax.get_children()[int(p[0])].set_color(\"tab:purple\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "207443f2-34d9-424a-a6d7-44707ef1488b",
"metadata": {},
"source": [
"#### Visualize best cut\n",
"\n",
"From the optimal bitstring, you can then visualize this cut on the original graph.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "33135970-8bc4-4fb2-ab87-08726a432ce4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# auxiliary function to plot graphs\n",
"def plot_result(G, x):\n",
" colors = [\"tab:grey\" if i == 0 else \"tab:purple\" for i in x]\n",
" pos, _default_axes = rx.spring_layout(G), plt.axes(frameon=True)\n",
" rx.visualization.mpl_draw(\n",
" G, node_color=colors, node_size=100, alpha=0.8, pos=pos\n",
" )\n",
"\n",
"\n",
"plot_result(graph, most_likely_bitstring)"
]
},
{
"cell_type": "markdown",
"id": "2119575f-f3cf-45bc-ae2b-93c046391eb6",
"metadata": {},
"source": [
"Now, calculate the value of the cut:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "2f6a73c4-f5ae-4647-a0dd-d77a13f66388",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The value of the cut is: 5\n"
]
}
],
"source": [
"def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:\n",
" assert len(x) == len(\n",
" list(graph.nodes())\n",
" ), \"The length of x must coincide with the number of nodes in the graph.\"\n",
" return sum(\n",
" x[u] * (1 - x[v]) + x[v] * (1 - x[u])\n",
" for u, v in list(graph.edge_list())\n",
" )\n",
"\n",
"\n",
"cut_value = evaluate_sample(most_likely_bitstring, graph)\n",
"print(\"The value of the cut is:\", cut_value)"
]
},
{
"cell_type": "markdown",
"id": "76a7241e",
"metadata": {},
"source": [
"For a graph this small, the true optimum is easy to brute-force, so you can double-check the results by comparing the QAOA result against the exact answer.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "0a3b5267",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Classical optimum (brute force): 5\n",
"QAOA cut value: 5\n"
]
}
],
"source": [
"# Classical baseline: enumerate all 2**n_small bitstrings and take the best cut.\n",
"def brute_force_max_cut(graph: rx.PyGraph) -> tuple[int, list[int]]:\n",
" n = len(list(graph.nodes()))\n",
" best_cut = -1\n",
" best_x: list[int] = []\n",
" for i in range(2**n):\n",
" x = [(i >> k) & 1 for k in range(n)]\n",
" cut = evaluate_sample(x, graph)\n",
" if cut > best_cut:\n",
" best_cut = int(cut)\n",
" best_x = x\n",
" return best_cut, best_x\n",
"\n",
"\n",
"classical_best, classical_x = brute_force_max_cut(graph)\n",
"print(f\"Classical optimum (brute force): {classical_best}\")\n",
"print(f\"QAOA cut value: {cut_value}\")"
]
},
{
"cell_type": "markdown",
"id": "large-scale-header",
"metadata": {},
"source": [
"## Large-scale hardware example\n",
"\n",
"You have access to many devices with over 100 qubits on IBM Quantum Platform. Select one on which to solve max-cut on a 100-node weighted graph. This is a \"utility-scale\" problem. The workflow follows the same steps as above, applied to a much larger graph.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "large-scale-steps-intro",
"metadata": {},
"source": [
"### End-to-end workflow at utility scale\n",
"\n",
"All four steps are shown below, applied to the 100-node graph. The structure is the same as the small-scale walkthrough: map, transpile, execute, post-process — but with a larger problem and split across the four cells below for clarity.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "large-scale-helpers",
"metadata": {},
"outputs": [],
"source": [
"# Precomputed parity lookup table: _PARITY[b] = +1 if popcount(b) is even, else -1.\n",
"# We use this to vectorize expectation-value evaluation across all Pauli terms.\n",
"_PARITY = np.array(\n",
" [-1 if bin(i).count(\"1\") % 2 else 1 for i in range(256)],\n",
" dtype=np.complex128,\n",
")\n",
"\n",
"\n",
"def evaluate_sparse_pauli(state: int, observable: SparsePauliOp) -> complex:\n",
" \"\"\"Expectation value of a SparsePauliOp on a single computational-basis state.\n",
"\n",
" For a Z-only observable (which QAOA cost Hamiltonians are, after the\n",
" QUBO-to-Hamiltonian mapping), the eigenvalue of each Pauli term on a\n",
" computational-basis state is simply (-1)**popcount(z_mask AND state),\n",
" i.e., the parity of the bitwise-AND of the term's Z-support and the\n",
" measured bitstring.\n",
"\n",
" This routine packs the Z-support of every Pauli term into bytes, ANDs\n",
" them against the measured state in a single vectorized op, and looks up\n",
" the parity in _PARITY. For a 100-qubit / ~hundreds-of-terms Hamiltonian\n",
" over 10_000 samples, this is dramatically faster than calling\n",
" SparsePauliOp.expectation_value per sample.\n",
" \"\"\"\n",
" packed_uint8 = np.packbits(observable.paulis.z, axis=1, bitorder=\"little\")\n",
" state_bytes = np.frombuffer(\n",
" state.to_bytes(packed_uint8.shape[1], \"little\"), dtype=np.uint8\n",
" )\n",
" reduced = np.bitwise_xor.reduce(packed_uint8 & state_bytes, axis=1)\n",
" return np.sum(observable.coeffs * _PARITY[reduced])\n",
"\n",
"\n",
"def best_solution(samples, hamiltonian):\n",
" \"\"\"Return the sampled bitstring (as int) with the lowest Hamiltonian cost.\"\"\"\n",
" min_cost = float(\"inf\")\n",
" min_sol = None\n",
" for bit_str in samples.keys():\n",
" candidate_sol = int(bit_str)\n",
" fval = evaluate_sparse_pauli(candidate_sol, hamiltonian).real\n",
" if fval <= min_cost:\n",
" min_cost = fval\n",
" min_sol = candidate_sol\n",
" return min_sol\n",
"\n",
"\n",
"def _plot_cdf(objective_values: dict, ax, color):\n",
" x_vals = sorted(objective_values.keys(), reverse=True)\n",
" y_vals = np.cumsum([objective_values[x] for x in x_vals])\n",
" ax.plot(x_vals, y_vals, color=color)\n",
"\n",
"\n",
"def plot_cdf(dist, ax, title):\n",
" _plot_cdf(dist, ax, \"C1\")\n",
" ax.vlines(min(list(dist.keys())), 0, 1, \"C1\", linestyle=\"--\")\n",
" ax.set_title(title)\n",
" ax.set_xlabel(\"Objective function value\")\n",
" ax.set_ylabel(\"Cumulative distribution function\")\n",
" ax.grid(alpha=0.3)\n",
"\n",
"\n",
"def samples_to_objective_values(samples, hamiltonian):\n",
" \"\"\"Convert the samples to values of the objective function.\"\"\"\n",
" objective_values = defaultdict(float)\n",
" for bit_str, prob in samples.items():\n",
" candidate_sol = int(bit_str)\n",
" fval = evaluate_sparse_pauli(candidate_sol, hamiltonian).real\n",
" objective_values[fval] += prob\n",
" return objective_values"
]
},
{
"cell_type": "markdown",
"id": "4dc0c8ff",
"metadata": {},
"source": [
"**Step 1**: Build the graph, cost Hamiltonian, and ansatz.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "94190344",
"metadata": {},
"outputs": [],
"source": [
"# Step 1: build the 100-node graph, cost Hamiltonian, and QAOA ansatz.\n",
"n_large = 100\n",
"graph_100 = rx.PyGraph()\n",
"graph_100.add_nodes_from(np.arange(0, n_large, 1))\n",
"elist = []\n",
"for edge in backend.coupling_map:\n",
" if edge[0] < n_large and edge[1] < n_large:\n",
" elist.append((edge[0], edge[1], 1.0))\n",
"graph_100.add_edges_from(elist)\n",
"\n",
"max_cut_paulis_100 = build_max_cut_paulis(graph_100)\n",
"cost_hamiltonian_100 = SparsePauliOp.from_sparse_list(\n",
" max_cut_paulis_100, n_large\n",
")\n",
"\n",
"circuit_100 = QAOAAnsatz(cost_operator=cost_hamiltonian_100, reps=1)\n",
"circuit_100.measure_all()"
]
},
{
"cell_type": "markdown",
"id": "c91f0c16",
"metadata": {},
"source": [
"**Step 2**: Transpile for the selected hardware backend.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "2b59da0e",
"metadata": {},
"outputs": [],
"source": [
"# Step 2: transpile for hardware.\n",
"pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n",
"candidate_circuit_100 = pm.run(circuit_100)"
]
},
{
"cell_type": "markdown",
"id": "eb2c4d66",
"metadata": {},
"source": [
"**Step 3**: Run the QAOA optimization loop inside a session, then sample.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "e5aceab3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" message: Return from COBYLA because the trust region radius reaches its lower bound.\n",
" success: True\n",
" status: 0\n",
" fun: -17.172689238986344\n",
" x: [ 2.574e+00 4.166e+00]\n",
" nfev: 28\n",
" maxcv: 0.0\n"
]
}
],
"source": [
"# Step 3: run the QAOA optimization loop on the device, then sample the\n",
"# final distribution with the optimized parameters.\n",
"initial_gamma = np.pi\n",
"initial_beta = np.pi / 2\n",
"init_params = [initial_beta, initial_gamma]\n",
"\n",
"objective_func_vals = [] # Global variable\n",
"with Session(backend=backend) as session:\n",
" estimator = Estimator(mode=session)\n",
" estimator.options.default_shots = 1000\n",
"\n",
" # Set simple error suppression/mitigation options\n",
" estimator.options.dynamical_decoupling.enable = True\n",
" estimator.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
" estimator.options.twirling.enable_gates = True\n",
" estimator.options.twirling.num_randomizations = \"auto\"\n",
" estimator.options.environment.job_tags = [\"TUT_QAOA\"]\n",
"\n",
" result = minimize(\n",
" cost_func_estimator,\n",
" init_params,\n",
" args=(candidate_circuit_100, cost_hamiltonian_100, estimator),\n",
" method=\"COBYLA\",\n",
" )\n",
" print(result)\n",
"\n",
"# Assign optimal parameters and sample the final distribution.\n",
"optimized_circuit_100 = candidate_circuit_100.assign_parameters(result.x)\n",
"\n",
"sampler = Sampler(mode=backend)\n",
"sampler.options.default_shots = 10000\n",
"\n",
"# Set simple error suppression/mitigation options\n",
"sampler.options.dynamical_decoupling.enable = True\n",
"sampler.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
"sampler.options.twirling.enable_gates = True\n",
"sampler.options.twirling.num_randomizations = \"auto\"\n",
"\n",
"# Add a unique tag to the job execution\n",
"sampler.options.environment.job_tags = [\"TUT_QAOA\"]\n",
"\n",
"pub = (optimized_circuit_100,)\n",
"job = sampler.run([pub], shots=int(1e4))\n",
"\n",
"counts_int = job.result()[0].data.meas.get_int_counts()\n",
"shots = sum(counts_int.values())\n",
"final_distribution_100_int = {\n",
" key: val / shots for key, val in counts_int.items()\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "7f0c5980",
"metadata": {},
"source": [
"**Step 4**: Post-process the sampled distribution to extract the best cut.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "010571f7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result bitstring: [1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0]\n",
"The value of the cut is: 156\n"
]
}
],
"source": [
"# Step 4: find the best-cost sample and evaluate its cut value.\n",
"best_sol_100 = best_solution(final_distribution_100_int, cost_hamiltonian_100)\n",
"best_sol_bitstring_100 = to_bitstring(int(best_sol_100), len(graph_100))\n",
"best_sol_bitstring_100.reverse()\n",
"\n",
"print(\"Result bitstring:\", best_sol_bitstring_100)\n",
"\n",
"cut_value_100 = evaluate_sample(best_sol_bitstring_100, graph_100)\n",
"print(\"The value of the cut is:\", cut_value_100)"
]
},
{
"cell_type": "markdown",
"id": "large-scale-convergence-md",
"metadata": {},
"source": [
"Check that the cost minimized in the optimization loop has converged, and visualize results.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "large-scale-viz",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot convergence\n",
"plt.figure(figsize=(12, 6))\n",
"plt.plot(objective_func_vals)\n",
"plt.xlabel(\"Iteration\")\n",
"plt.ylabel(\"Cost\")\n",
"plt.show()\n",
"\n",
"# Visualize the cut\n",
"plot_result(graph_100, best_sol_bitstring_100)\n",
"\n",
"# Plot cumulative distribution function\n",
"result_dist = samples_to_objective_values(\n",
" final_distribution_100_int, cost_hamiltonian_100\n",
")\n",
"fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
"plot_cdf(result_dist, ax, backend.name)"
]
},
{
"cell_type": "markdown",
"id": "69ebc85b-6a29-4671-8d16-1ac97f089607",
"metadata": {},
"source": [
"## Next steps\n",
"\n",
"If you found this work interesting, you might be interested in the following material:\n",
"\n",
"\n",
" * [Advanced techniques for QAOA](/docs/tutorials/advanced-techniques-for-qaoa) — explores advanced strategies for improving QAOA performance\n",
" * [Multi-objective optimization challenge](https://github.com/qiskit-community/qdc-challenges-2025/blob/main/challenges/Track_B/qmoo/qmoo_qdc25.ipynb) — put your skills to the test with this community challenge on multi-objective quantum optimization\n",
" * [Transpilation documentation](/docs/guides/transpile) for fine-tuning circuit optimization\n",
" * [Error suppression and mitigation](/docs/guides/error-mitigation-and-suppression-techniques) for improving hardware results\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3"
},
"hours": 1,
"qpuSeconds": 1320
},
"nbformat": 4,
"nbformat_minor": 5
}