{
"cells": [
{
"cell_type": "markdown",
"id": "eb419bc5-908c-4f0d-a4ff-c4d13f04e332",
"metadata": {
"tags": []
},
"source": [
"---\n",
"title: Ground-state energy estimation of the Heisenberg chain with VQE\n",
"description: Build, deploy, and run a Qiskit pattern for simulating a Heisenberg chain and estimating its ground-state energy.\n",
"---\n",
"\n",
"{/* cspell:ignore hyperparameters, forall, nabla, nparams */}\n",
"\n",
"# Ground-state energy estimation of the Heisenberg chain with VQE\n",
"\n",
"*Usage estimate: 37 minutes on a Heron processor (NOTE: This is an estimate only. Your runtime might vary.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "49d868bf",
"metadata": {},
"source": [
"## Background\n",
"\n",
"This tutorial shows how to build, deploy, and run a development workflow called a [Qiskit pattern](/docs/guides/intro-to-patterns) for simulating a Heisenberg chain and estimating its ground-state energy using the SPSA optimizer.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "bc52f763",
"metadata": {},
"source": [
"## Requirements\n",
"\n",
"Before starting this tutorial, ensure that you have the following installed:\n",
"\n",
"* Qiskit SDK v2.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
"* Qiskit Runtime v0.44 or later (`pip install qiskit-ibm-runtime`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a46e9e3e",
"metadata": {},
"source": [
"## Setup\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e7754922",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from typing import Sequence\n",
"\n",
"from qiskit import QuantumCircuit\n",
"from qiskit.quantum_info import SparsePauliOp\n",
"from qiskit.primitives import BaseEstimatorV2\n",
"from qiskit.circuit.library import XGate\n",
"from qiskit.circuit.library import efficient_su2\n",
"from qiskit.transpiler import PassManager\n",
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"from qiskit.transpiler.passes.scheduling import (\n",
" ALAPScheduleAnalysis,\n",
" PadDynamicalDecoupling,\n",
")\n",
"from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2\n",
"\n",
"\n",
"def visualize_results(results):\n",
" plt.plot(results[\"cost_history\"], lw=2)\n",
" plt.xlabel(\"Number of function evaluations\")\n",
" plt.ylabel(\"Energy\")\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"id": "132fb15f-10b4-4d7e-83d8-f512a6f675d1",
"metadata": {},
"source": [
"## Step 1: Map classical inputs to a quantum problem\n",
"\n",
"* Input: Number of spins\n",
"* Output: Ansatz and Hamiltonian modeling the Heisenberg chain\n",
"\n",
"Construct an ansatz and Hamiltonian that model a 10-spin Heisenberg chain. First, we import some generic packages and create some helper functions.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7e8d2f10-f1d6-4ec2-bac9-9db23499c9e1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"num_spins = 10\n",
"ansatz = efficient_su2(num_qubits=num_spins, reps=2)\n",
"\n",
"service = QiskitRuntimeService(\n",
" channel=\"ibm_cloud\",\n",
" token=\"\", # Replace with your actual API token\n",
" instance=\"\", # Replace with your instance name if needed\n",
")\n",
"backend = service.least_busy(\n",
" operational=True, min_num_qubits=num_spins, simulator=False\n",
")\n",
"\n",
"coupling = backend.target.build_coupling_map()\n",
"reduced_coupling = coupling.reduce(list(range(num_spins)))\n",
"\n",
"edge_list = reduced_coupling.graph.edge_list()\n",
"ham_list = []\n",
"\n",
"for edge in edge_list:\n",
" ham_list.append((\"ZZ\", edge, 0.5))\n",
" ham_list.append((\"YY\", edge, 0.5))\n",
" ham_list.append((\"XX\", edge, 0.5))\n",
"\n",
"for qubit in reduced_coupling.physical_qubits:\n",
" ham_list.append((\"Z\", [qubit], np.random.random() * 2 - 1))\n",
"\n",
"hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)\n",
"\n",
"ansatz.draw(\"mpl\", style=\"iqp\")"
]
},
{
"cell_type": "markdown",
"id": "ab79119b-5e56-49d8-a20e-1c8e665baec0",
"metadata": {},
"source": [
"## Step 2: Optimize problem for quantum hardware execution\n",
"\n",
"* Input: Abstract circuit, observable\n",
"* Output: Target circuit and observable, optimized for the selected QPU\n",
"\n",
"Use the `generate_preset_pass_manager` function from Qiskit to automatically generate an optimization routine for our circuit with respect to the selected QPU. We choose `optimization_level=3`, which provides the highest level of optimization of the preset pass managers. We also include `ALAPScheduleAnalysis` and `PadDynamicalDecoupling` scheduling passes to suppress decoherence errors.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "a0a5f1c8-5c31-4d9f-ae81-37bd67271d44",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"target = backend.target\n",
"pm = generate_preset_pass_manager(optimization_level=3, target=target)\n",
"pm.scheduling = PassManager(\n",
" [\n",
" ALAPScheduleAnalysis(durations=target.durations()),\n",
" PadDynamicalDecoupling(\n",
" durations=target.durations(),\n",
" dd_sequence=[XGate(), XGate()],\n",
" pulse_alignment=target.pulse_alignment,\n",
" ),\n",
" ]\n",
")\n",
"isa_ansatz = pm.run(ansatz)\n",
"isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)\n",
"isa_ansatz.draw(\"mpl\", scale=0.6, style=\"iqp\", fold=-1, idle_wires=False)"
]
},
{
"cell_type": "markdown",
"id": "9e889d0b-30b5-4e6b-84c9-d1f096abf132",
"metadata": {},
"source": [
"## Step 3: Execute using Qiskit primitives\n",
"\n",
"* Input: Target circuit and observable\n",
"* Output: Results of optimization\n",
"\n",
"Minimize the estimated ground-state energy of the system by optimizing the circuit parameters. Use the `Estimator` primitive from Qiskit Runtime to evaluate the cost function during optimization.\n",
"\n",
"Since we optimized the circuit for the backend in Step 2, we can avoid doing transpilation on the Runtime server by setting `skip_transpilation=True` and passing the optimized circuit. For this demo, we will run on a QPU using `qiskit-ibm-runtime` primitives. To run with `qiskit` statevector-based primitives, replace the block of code using Qiskit Runtime primitives with the commented block.\n",
"\n",
"In this tutorial we use Simultaneous Perturbation Stochastic Approximation (SPSA), which is a gradient-based optimizer. Next we give a brief introduction to it, and provide the code to implement SPSA using Qiskit v2.0.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "64f97e6f",
"metadata": {},
"source": [
"### Introducing SPSA\n",
"\n",
"Simultaneous Perturbation Stochastic Approximation (SPSA) [\\[1\\]](#references) is an optimization algorithm that approximates the entire gradient vector using only two function calls at each iteration. Let $f:\\mathbb{R}^p\\rightarrow \\mathbb{R}$ be the cost function with $p$ parameters to be optimized, and $x_i\\in \\mathbb{R}^p$ be the parameter vector at the $i^{th}$ step of the iteration. To compute the gradient, a random vector $\\Delta_i$ of size $p$ is created, where each element $\\Delta_{ij}$, $\\forall$ $j\\in \\{1,2,...,p\\}$, is uniformly sampled from $\\{-1, 1\\}$. Next, each element of the random vector $\\Delta_i$ is multiplied with a small value $c_i$ to create a random perturbation. The gradient is then estimated as\n",
"\n",
"$[\\nabla f(x_i)]_j \\approx \\frac{f(x_i + c_i \\Delta_i) - f(x_i - c_i \\Delta_i)}{2c_i\\Delta_{ij}}.$\n",
"\n",
"Intuitively, since a random perturbation is applied during the gradient estimation, it is expected that small deviations in the exact values of $f$ coming from noise can be tolerated and accounted for. In fact, SPSA is particularly known to be robust against noise, and requires only two hardware calls for each iteration. It is, therefore, one of the highly preferred optimizers for implementing variational algorithms.\n",
"\n",
"In this tutorial, the hyperparameters for the $i^{th}$ iteration, $a_i$ and $c_i$, are computed as $a_i = a/(A + i + 1)^\\alpha$ and $c_i = c / (i+1)^\\gamma$, where the constant values are taken as $A = 30$, $\\alpha = 0.9$, $a = 0.3$, $c = 0.1$, and $\\gamma = 0.4$. These values are selected from [\\[2\\]](#references). Appropriate tuning of hyperparameters is necessary for extracting a good performance out of SPSA.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "73a9352c",
"metadata": {},
"outputs": [],
"source": [
"def spsa(\n",
" fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100\n",
"):\n",
" nparams = len(x0)\n",
" x = np.copy(x0)\n",
"\n",
" for i in range(maxiter):\n",
" a_i = a / (A + i + 1) ** alpha\n",
" c_i = c / (i + 1) ** gamma\n",
" delta_i = np.random.choice([-1, 1], nparams)\n",
"\n",
" # two hardware calls\n",
" eval_1 = fun(x + c_i * delta_i, *args)\n",
" eval_2 = fun(x - c_i * delta_i, *args)\n",
"\n",
" # compute the gradient and update the parameters\n",
" grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)\n",
" x = x - a_i * grad\n",
"\n",
" return x"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "32ca3b6a",
"metadata": {},
"outputs": [],
"source": [
"def cost_func(\n",
" params: Sequence,\n",
" ansatz: QuantumCircuit,\n",
" hamiltonian: SparsePauliOp,\n",
" estimator: BaseEstimatorV2,\n",
" cost_history_dict: dict,\n",
") -> float:\n",
" \"\"\"Ground state energy evaluation.\"\"\"\n",
" energy = (\n",
" estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs\n",
" )\n",
"\n",
" cost_history_dict[\"iters\"] += 1\n",
" cost_history_dict[\"prev_vector\"] = list(params)\n",
" cost_history_dict[\"cost_history\"].append(float(energy[0]))\n",
"\n",
" print(\n",
" f\"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]\",\n",
" end=\"\\r\",\n",
" )\n",
"\n",
" return energy\n",
"\n",
"\n",
"def solve(x0, isa_ansatz, isa_observable, maxiter=150):\n",
" cost_history_dict = {\n",
" \"prev_vector\": None,\n",
" \"iters\": 0,\n",
" \"cost_history\": [],\n",
" \"y_min\": None,\n",
" }\n",
"\n",
" # Evaluate the problem using a QPU via Qiskit IBM Runtime\n",
" with Session(backend=backend) as session:\n",
" estimator = EstimatorV2(mode=session)\n",
" estimator.skip_transpilation = True\n",
" x_opt = spsa(\n",
" cost_func,\n",
" x0=x0,\n",
" args=(isa_ansatz, isa_observable, estimator, cost_history_dict),\n",
" maxiter=maxiter,\n",
" )\n",
"\n",
" y_min = cost_func(\n",
" x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict\n",
" )\n",
"\n",
" return y_min, cost_history_dict"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "f418b372",
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(42)\n",
"num_params = ansatz.num_parameters\n",
"params = 2 * np.pi * np.random.random(num_params)"
]
},
{
"cell_type": "markdown",
"id": "3bf42923",
"metadata": {},
"source": [
"Here we set the `maxiter = 50`. Note that since each iteration requires two calls to the function to compute the gradient, the total number of function calls will be $2 \\times \\text{maxiter}$. The `maxiter` can be increased to any higher value for better energy estimation.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1732ce37",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fx Iters. done: 101 [Current cost: -2.19621]\r"
]
}
],
"source": [
"maxiter = 50\n",
"spsa_min, spsa_history = solve(\n",
" params, isa_ansatz, isa_observable, maxiter=maxiter\n",
")"
]
},
{
"cell_type": "markdown",
"id": "33abbb3f-6245-4610-a05d-e2bc4cc551f0",
"metadata": {},
"source": [
"## Step 4: Post-process and return result in desired classical format\n",
"\n",
"* Input: Ground-state energy estimates during optimization\n",
"* Output: Estimated ground-state energy\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e5b58771-d543-4e75-9746-fbc7b28e4360",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimated ground state energy: [-2.19621239]\n"
]
}
],
"source": [
"print(f\"Estimated ground state energy: {spsa_min}\")"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "ecd7762a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"results = {\n",
" \"spsa\": spsa_history,\n",
"}\n",
"\n",
"visualize_results(spsa_history)"
]
},
{
"cell_type": "markdown",
"id": "217a9379",
"metadata": {},
"source": [
"## References\n",
"\n",
"\\[1] Spall, J. C. (2002). Implementation of the simultaneous perturbation algorithm for stochastic optimization.\n",
"IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.\n",
"\n",
"\\[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning: an open-source library for quantum machine learning tasks at scale on quantum hardware and classical simulators. arXiv:2505.17756.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8f14492c",
"metadata": {},
"source": [
"## Next steps\n",
"\n",
"\n",
" If you found this work interesting, you might be interested in the following material:\n",
"\n",
" * Find the ground-state energy of a sparse Hamiltonian by following the [SQD tutorial](/docs/tutorials/sample-based-quantum-diagonalization)\n",
" * Take the [Variational algorithm design](/learning/courses/variational-algorithm-design) course in IBM Learning\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"
}
},
"nbformat": 4,
"nbformat_minor": 5
}