Odhad energie základního stavu Heisenbergova řetězce pomocí VQE

Odhadovaná spotřeba: 37 minut na procesoru Heron (POZNÁMKA: Jedná se pouze o odhad. Tvůj skutečný čas běhu se může lišit.)

Výstupy z učení

Po dokončení tohoto tutoriálu bys měl(a) rozumět následujícím informacím:

• Jak modelovat Heisenbergův spinový řetězec jako kvantový hamiltonián pomocí Qiskitu
• Jak používat optimalizátor SPSA k odhadu energie základního stavu kvantového systému
• Jak spouštět variační pracovní postupy na kvantovém hardwaru IBM® pomocí primitiv Qiskit Runtime a relací

Předpoklady

Doporučujeme se předem seznámit s těmito tématy:

• Základy kvantových informací
• Úvod do Qiskit vzorů
• Návrh variačních algoritmů

Pozadí

Heisenbergův spinový řetězec je jedním z nejrozšířeněji studovaných modelů v fyzice kondenzované hmoty a kvantovém magnetismu. Popisuje jednorozměrnou mřížku vzájemně působících kvantových spinů, kde sousední spiny jsou vázány prostřednictvím výměnných interakcí. Hamiltonián pro izotropní Heisenbergův model s vnějším magnetickým polem je dán:

$H = \sum_{\langle i,j \rangle} \left( J_x X_i X_j + J_y Y_i Y_j + J_z Z_i Z_j \right) + \sum_{i} h_i Z_i,$

kde $X_i$, $Y_i$ a $Z_i$ jsou Pauliho operátory působící na místě $i$, suma $\langle i,j \rangle$ probíhá přes nejbližší sousední páry, $J_x = J_y = J_z = 0.5$ jsou konstanty výměnné vazby (v tomto tutoriálu izotropní) a $h_i$ představuje vnější magnetické pole závislé na místě. V tomto tutoriálu jsou hodnoty magnetického pole náhodně vzorkovány z intervalu $[-1, 1]$. Poznamenejme, že v níže uvedené implementaci je množina párů „nejbližších sousedů" určena nativní vazbou hardwarového backendu pro prvních $N$ qubitů, která nemusí tvořit přísný lineární řetězec v závislosti na topologii zařízení.

Porozumění energii základního stavu tohoto hamiltoniánu má v fyzice zásadní význam. Základní stav zakódovává informace o kvantových fázových přechodech, struktuře entanglementu a magnetickém uspořádání. Klasicky se výpočet přesné energie základního stavu stává nesouditelným s rostoucím počtem spinů, protože dimenze Hilbertova prostoru roste exponenciálně jako $2^N$ pro $N$ spinů. To z něj činí přirozený kandidát pro kvantovou simulaci.

Variační kvantový eigensolver (VQE) je hybridní kvantově-klasický algoritmus navržený k odhadu energie základního stavu hamiltoniánu. Funguje tak, že na kvantovém počítači připraví parametrizovaný kvantový stav $|\psi(\theta)\rangle$ (zvaný ansatz) a změří střední hodnotu $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Klasický optimalizátor pak iterativně upravuje parametry $\theta$ tak, aby tuto energii minimalizoval, přičemž využívá variačního principu, který zaručuje, že naměřená energie je vždy horním odhadem skutečné energie základního stavu.

V tomto tutoriálu používáme ansatz efficient_su2 z knihovny obvodů Qiskitu, který konstruuje vrstvy jednobitových rotací a entanglujících hradel. Optimalizace je prováděna pomocí algoritmu Simultaneous Perturbation Stochastic Approximation (SPSA), který je vhodný pro hlučný kvantový hardware, protože odhaduje gradienty pomocí pouhých dvou vyhodnocení funkce za iteraci bez ohledu na počet parametrů.

Požadavky

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalováno následující:

• Qiskit SDK v2.0 nebo novější, s podporou vizualizace
• Qiskit Runtime v0.44 nebo novější (pip install qiskit-ibm-runtime)

Nastavení

# 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()

Příklad v malém měřítku

V této části postupujeme krok za krokem vzorem Qiskit v malém měřítku a vysvětlujeme klíčové komponenty při sestavování pracovního postupu.

Krok 1: Mapování klasických vstupů na kvantový problém

• Vstup: Počet spinů
• Výstup: Ansatz a hamiltonián modelující Heisenbergův řetězec

Sestrojíme ansatz a hamiltonián, které modelují 10-spinový Heisenbergův řetězec. V tomto kroku sestavíme 10-spinový Heisenbergův hamiltonián na propojovací mapě nejméně vytíženého backendu a připravíme ansatz efficient_su2.

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

service = QiskitRuntimeService()
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")

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru

• Vstup: Abstraktní Circuit, observabilní veličina
• Výstup: Cílový Circuit a observabilní veličina, optimalizované pro vybrané QPU

Použijeme funkci generate_preset_pass_manager z Qiskitu k automatickému vygenerování optimalizační rutiny pro náš obvod vzhledem k vybranému QPU. Volíme optimization_level=3, což poskytuje nejvyšší úroveň optimalizace z přednastavených pass managerů. Také zahrneme plánovací průchody ALAPScheduleAnalysis a PadDynamicalDecoupling pro potlačení chyb způsobených dekoherencí.

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)

Krok 3: Spuštění pomocí Qiskit primitiv

• Vstup: Cílový Circuit a observabilní veličina
• Výstup: Výsledky optimalizace

Minimalizujeme odhadovanou energii základního stavu systému optimalizací parametrů Circuitu. Použijeme primitivum Estimator z Qiskit Runtime k vyhodnocení účelové funkce během optimalizace.

Protože jsme v kroku 2 optimalizovali Circuit pro backend, můžeme se vyhnout transpilaci na serveru Runtime nastavením skip_transpilation=True a předáním optimalizovaného Circuitu. Pro tuto ukázku spustíme výpočet na QPU pomocí primitiv qiskit-ibm-runtime. Pro spuštění s primitivy qiskit založenými na stavovém vektoru nahraď blok kódu používající primitiva Qiskit Runtime komentovaným blokem.

V tomto tutoriálu používáme Simultaneous Perturbation Stochastic Approximation (SPSA), což je gradientní optimalizátor. Níže uvádíme stručný úvod do SPSA a kód pro jeho implementaci pomocí Qiskit v2.0.

Představení SPSA

Simultaneous Perturbation Stochastic Approximation (SPSA) [1] je optimalizační algoritmus, který aproximuje celý gradientní vektor pomocí pouhých dvou volání funkce v každé iteraci. Nechť $f:\mathbb{R}^p\rightarrow \mathbb{R}$ je účelová funkce s $p$ parametry k optimalizaci a $x_i\in \mathbb{R}^p$ je vektor parametrů v $i$-tém kroku iterace. Pro výpočet gradientu se vytvoří náhodný vektor $\Delta_i$ velikosti $p$, kde každý prvek $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, je rovnoměrně vzorkován z $\{-1, 1\}$. Každý prvek náhodného vektoru $\Delta_i$ je poté vynásoben malou hodnotou $c_i$, čímž vznikne náhodná perturbace. Gradient je pak odhadnut jako

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

Intuitivně, protože se při odhadování gradientu aplikuje náhodná perturbace, lze očekávat toleranci vůči malým odchylkám v přesných hodnotách $f$ způsobeným šumem. SPSA je zejména známé svou robustností vůči šumu a vyžaduje pouze dvě volání hardwaru v každé iteraci. Je proto jedním z nejvíce preferovaných optimalizátorů pro implementaci variačních algoritmů.

V tomto tutoriálu jsou hyperparametry pro $i$-tou iteraci, $a_i$ a $c_i$, vypočteny jako

$a_i = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{c}{(i+1)^\gamma},$

kde konstantní hodnoty jsou $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ a $\gamma = 0.4$. Tyto hodnoty jsou převzaty z [2]. Pro dobré výsledky SPSA je nutné vhodné ladění hyperparametrů.

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
        estimator.options.environment.job_tags = ["TUT_HSVQE"]
        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)

Zde nastavíme maxiter = 50. Protože každá iterace vyžaduje dvě volání funkce pro výpočet gradientu, celkový počet volání funkce bude $2 \times \text{maxiter}$. Hodnotu maxiter lze zvýšit pro lepší odhad energie.

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

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

Krok 4: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu

• Vstup: Odhady energie základního stavu během optimalizace
• Výstup: Odhadovaná energie základního stavu

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

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Příklad na velkém hardwaru

Příklad na velkém hardwaru není v tomto tutoriálu zahrnut. S rostoucím počtem qubitů narážíme u VQE na zásadní problémy způsobené fenoménem barren plateau: gradient účelové funkce klesá exponenciálně s velikostí systému, takže optimalizace je pro velké obvody prakticky neproveditelná. V kombinaci s hardwarovým šumem to znamená, že škálování VQE na větší spinové řetězce nepřináší spolehlivě reprodukovatelné výsledky. Přístupy, které tyto omezení překonávají, najdeš v části Další kroky níže.

Výzva

Nyní, když máš funkční implementaci VQE pro Heisenbergův řetězec, zkus následující:

1. Experimentuj s hloubkou ansatzu: Změň parametr reps ve funkci efficient_su2 (vyzkoušej například reps=1 a reps=3). Jak hloubka ansatzu ovlivňuje odhadovanou energii základního stavu a rychlost konvergence? Kdy začínáš pozorovat klesající výnosy nebo nestabilitu?
2. Lad hyperparametry SPSA: Uprav parametry plánu učební rychlosti (a, c, alpha, gamma, A) a sleduj, jak ovlivňují konvergenci. Dokážeš najít konfiguraci, která konverguje rychleji než výchozí hodnoty použité zde?
3. Porovnej topologie propojení: Místo nativní propojovací mapy backendu zkus sestavit jednoduchý lineární řetězec nejbližších sousedů a porovnej výsledky. Jak konektivita fyzického hardwaru ovlivňuje hloubku transpilovaného obvodu a výsledný odhad energie?

Reference

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

Další kroky

Doporučení

Pokud tě tato práce zaujala, možná tě budou zajímat následující materiály:

• Vyzkoušej Sample-based Quantum Diagonalization (SQD): Jak bylo ukázáno v tomto tutoriálu, VQE čelí při větším měřítku problémům kvůli barren plateaus a vysoké zátěži měřením. IBM vyvinulo Sample-based Quantum Diagonalization (SQD) jako škálovatelnější alternativu. Na rozdíl od VQE se SQD zcela vyhýbá variační optimalizaci; místo toho kvantový počítač generuje vzorky a klasický počítač projektuje hamiltonián na podprostor rozepjatý těmito vzorky a diagonalizuje jej. To poskytuje horní odhad energie základního stavu s výrazně menším počtem měření a bez náchylnosti k barren plateaus. Postupuj podle tutoriálu SQD a uvidíš tento přístup v praxi.
• Prozkoumej kurz Quantum Diagonalization Algorithms: Prohloubíš porozumění jak VQE, tak SQD, včetně jejich kompromisů, v kurzu Quantum diagonalization algorithms na IBM Quantum Learning.