Оцінка енергії основного стану ланцюга Гайзенберга за допомогою VQE

Орієнтовне використання: 37 хвилин на процесорі Heron (ПРИМІТКА: це лише орієнтовна оцінка. Ваш час виконання може відрізнятися.)

Передумови

У цьому посібнику показано, як створити та запустити процес розробки, що називається шаблоном Qiskit, для моделювання ланцюга Гайзенберга та оцінки його енергії основного стану за допомогою оптимізатора SPSA.

Вимоги

Перед початком роботи з цим посібником переконайся, що у тебе встановлено наступне:

• Qiskit SDK v2.0 або новішої версії, з підтримкою візуалізації
• Qiskit Runtime v0.44 або новішої версії (pip install qiskit-ibm-runtime)

Налаштування

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-ibm-runtime

import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimatorV2
from qiskit.circuit.library import XGate
from qiskit.circuit.library import efficient_su2
from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.scheduling import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2

def visualize_results(results):
    plt.plot(results["cost_history"], lw=2)
    plt.xlabel("Number of function evaluations")
    plt.ylabel("Energy")
    plt.show()

Крок 1: Перетворення класичних вхідних даних у квантову задачу

• Вхідні дані: кількість спінів
• Вихідні дані: анзац та гамільтоніан, що моделюють ланцюг Гайзенберга

Побудуйте анзац та гамільтоніан, які моделюють ланцюг Гайзенберга з 10 спінами. Спочатку ми імпортуємо деякі загальні пакети та створюємо кілька допоміжних функцій.

num_spins = 10
ansatz = efficient_su2(num_qubits=num_spins, reps=2)

service = QiskitRuntimeService(
    channel="ibm_cloud",
    token="<YOUR_API_TOKEN>",  # Replace with your actual API token
    instance="<YOUR_INSTANCE_NAME>",  # Replace with your instance name if needed
)
backend = service.least_busy(
    operational=True, min_num_qubits=num_spins, simulator=False
)

coupling = backend.target.build_coupling_map()
reduced_coupling = coupling.reduce(list(range(num_spins)))

edge_list = reduced_coupling.graph.edge_list()
ham_list = []

for edge in edge_list:
    ham_list.append(("ZZ", edge, 0.5))
    ham_list.append(("YY", edge, 0.5))
    ham_list.append(("XX", edge, 0.5))

for qubit in reduced_coupling.physical_qubits:
    ham_list.append(("Z", [qubit], np.random.random() * 2 - 1))

hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)

ansatz.draw("mpl", style="iqp")

Крок 2: Оптимізація задачі для виконання на квантовому обладнанні

• Вхідні дані: абстрактна схема, спостережувана величина
• Вихідні дані: цільова схема та спостережувана величина, оптимізовані для обраного QPU

Використовуйте функцію generate_preset_pass_manager з Qiskit для автоматичної генерації процедури оптимізації нашої схеми відносно обраного QPU. Ми обираємо optimization_level=3, що забезпечує найвищий рівень оптимізації серед попередньо налаштованих менеджерів проходів. Ми також додаємо проходи планування ALAPScheduleAnalysis та PadDynamicalDecoupling для придушення помилок декогеренції.

target = backend.target
pm = generate_preset_pass_manager(optimization_level=3, target=target)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(durations=target.durations()),
        PadDynamicalDecoupling(
            durations=target.durations(),
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)
isa_ansatz = pm.run(ansatz)
isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)
isa_ansatz.draw("mpl", scale=0.6, style="iqp", fold=-1, idle_wires=False)

Крок 3: Виконання за допомогою примітивів Qiskit

• Вхідні дані: цільова схема та спостережувана величина
• Вихідні дані: результати оптимізації

Мінімізуйте оцінену енергію основного стану системи шляхом оптимізації параметрів схеми. Використовуйте примітив Estimator з Qiskit Runtime для обчислення функції вартості під час оптимізації.

Оскільки ми оптимізували схему для бекенду на кроці 2, ми можемо уникнути транспіляції на сервері Runtime, встановивши skip_transpilation=True та передавши оптимізовану схему. Для цієї демонстрації ми будемо виконувати обчислення на QPU за допомогою примітивів qiskit-ibm-runtime. Щоб виконати обчислення з примітивами на основі вектора стану qiskit, замініть блок коду, що використовує примітиви Qiskit Runtime, на закоментований блок.

У цьому посібнику ми використовуємо метод Simultaneous Perturbation Stochastic Approximation (SPSA) — оптимізатор на основі градієнта. Нижче ми коротко його представимо та надамо код для реалізації SPSA за допомогою Qiskit v2.0.

Знайомство з SPSA

Simultaneous Perturbation Stochastic Approximation (SPSA) [1] — це алгоритм оптимізації, який апроксимує весь вектор градієнта, використовуючи лише два виклики функції на кожній ітерації. Нехай $f:\mathbb{R}^p\rightarrow \mathbb{R}$ — функція вартості з $p$ параметрами для оптимізації, а $x_i\in \mathbb{R}^p$ — вектор параметрів на $i$-му кроці ітерації. Для обчислення градієнта створюється випадковий вектор $\Delta_i$ розміру $p$, де кожен елемент $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, рівномірно вибирається з $\{-1, 1\}$. Далі кожен елемент випадкового вектора $\Delta_i$ множиться на мале значення $c_i$ для створення випадкового збурення. Градієнт тоді оцінюється як

$[\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}}.$

Інтуїтивно, оскільки випадкове збурення застосовується під час оцінювання градієнта, очікується, що невеликі відхилення в точних значеннях $f$, що виникають від шуму, можуть бути допустимими та враховуватися. Насправді SPSA особливо відомий своєю стійкістю до шуму та вимагає лише двох апаратних викликів для кожної ітерації. Тому він є одним з найбільш переважних оптимізаторів для реалізації варіаційних алгоритмів.

У цьому посібнику гіперпараметри для $i$-ї ітерації, $a_i$ та $c_i$, обчислюються як $a_i = a/(A + i + 1)^\alpha$ та $c_i = c / (i+1)^\gamma$, де константні значення взяті як $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ та $\gamma = 0.4$. Ці значення вибрані з [2]. Відповідне налаштування гіперпараметрів необхідне для отримання хорошої продуктивності від SPSA.

def spsa(
    fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100
):
    nparams = len(x0)
    x = np.copy(x0)

    for i in range(maxiter):
        a_i = a / (A + i + 1) ** alpha
        c_i = c / (i + 1) ** gamma
        delta_i = np.random.choice([-1, 1], nparams)

        # two hardware calls
        eval_1 = fun(x + c_i * delta_i, *args)
        eval_2 = fun(x - c_i * delta_i, *args)

        # compute the gradient and update the parameters
        grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)
        x = x - a_i * grad

    return x

def cost_func(
    params: Sequence,
    ansatz: QuantumCircuit,
    hamiltonian: SparsePauliOp,
    estimator: BaseEstimatorV2,
    cost_history_dict: dict,
) -> float:
    """Ground state energy evaluation."""
    energy = (
        estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs
    )

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = list(params)
    cost_history_dict["cost_history"].append(float(energy[0]))

    print(
        f"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]",
        end="\r",
    )

    return energy

def solve(x0, isa_ansatz, isa_observable, maxiter=150):
    cost_history_dict = {
        "prev_vector": None,
        "iters": 0,
        "cost_history": [],
        "y_min": None,
    }

    # Evaluate the problem using a QPU via Qiskit IBM Runtime
    with Session(backend=backend) as session:
        estimator = EstimatorV2(mode=session)
        estimator.skip_transpilation = True
        x_opt = spsa(
            cost_func,
            x0=x0,
            args=(isa_ansatz, isa_observable, estimator, cost_history_dict),
            maxiter=maxiter,
        )

        y_min = cost_func(
            x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict
        )

    return y_min, cost_history_dict

np.random.seed(42)
num_params = ansatz.num_parameters
params = 2 * np.pi * np.random.random(num_params)

Тут ми встановлюємо maxiter = 50. Зверніть увагу, що оскільки кожна ітерація вимагає двох викликів функції для обчислення градієнта, загальна кількість викликів функції становитиме $2 \times \text{maxiter}$. Значення maxiter можна збільшити до будь-якого більшого значення для кращої оцінки енергії.

maxiter = 50
spsa_min, spsa_history = solve(
    params, isa_ansatz, isa_observable, maxiter=maxiter
)

Fx Iters. done: 101 [Current cost: -2.19621]

Крок 4: Постобробка та повернення результату у бажаному класичному форматі

• Вхідні дані: оцінки енергії основного стану під час оптимізації
• Вихідні дані: оцінена енергія основного стану

print(f"Estimated ground state energy: {spsa_min}")

Estimated ground state energy: [-2.19621239]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Посилання

[1] Spall, J. C. (2002). Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[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.

Наступні кроки

Рекомендації

Якщо ця робота вас зацікавила, вас можуть зацікавити такі матеріали:

• Знайдіть енергію основного стану розрідженого гамільтоніана, дотримуючись посібника SQD
• Пройдіть курс Variational algorithm design на IBM Learning