Estymacja energii stanu podstawowego łańcucha Heisenberga metodą VQE

Szacowane zużycie zasobów: 37 minut na procesorze Heron (UWAGA: To tylko szacunek. Twój czas wykonania może się różnić.)

Cele dydaktyczne

Po ukończeniu tego samouczka powinieneś rozumieć następujące zagadnienia:

• Jak modelować łańcuch spinów Heisenberga jako Hamiltonian kwantowy przy użyciu Qiskit
• Jak używać optymalizatora SPSA do estymacji energii stanu podstawowego układu kwantowego
• Jak wykonywać wariacyjne przepływy pracy na sprzęcie kwantowym IBM® przy użyciu prymitywów Qiskit Runtime i sesji

Wymagania wstępne

Zaleca się zapoznanie z następującymi tematami:

• Podstawy informacji kwantowej
• Wprowadzenie do wzorców Qiskit
• Projektowanie algorytmów wariacyjnych

Tło

Łańcuch spinów Heisenberga jest jednym z najszerzej badanych modeli w fizyce materii skondensowanej i magnetyzmie kwantowym. Opisuje jednowymiarową sieć oddziałujących spinów kwantowych, gdzie sąsiednie spiny są sprzężone poprzez oddziaływania wymienne. Hamiltonian izotropowego modelu Heisenberga z zewnętrznym polem magnetycznym ma postać:

$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,$

gdzie $X_i$, $Y_i$ i $Z_i$ to operatory Pauliego działające na węzeł $i$, suma $\langle i,j \rangle$ przebiega po parach najbliższych sąsiadów, $J_x = J_y = J_z = 0.5$ to stałe sprzężenia wymiennego (izotropowe w tym samouczku), a $h_i$ reprezentuje zależne od węzła zewnętrzne pole magnetyczne. W tym samouczku wartości pola magnetycznego są losowo próbkowane z przedziału $[-1, 1]$. Należy zauważyć, że w poniższej implementacji zbiór par „najbliższych sąsiadów" jest wyznaczany przez natywne sprzężenie backendu sprzętowego wśród pierwszych $N$ kubitów, co w zależności od topologii urządzenia może nie tworzyć ścisłego łańcucha liniowego.

Zrozumienie energii stanu podstawowego tego Hamiltonianu ma fundamentalne znaczenie w fizyce. Stan podstawowy koduje informacje o kwantowych przejściach fazowych, strukturze splątania i porządkowaniu magnetycznym. Klasyczne obliczenie dokładnej energii stanu podstawowego staje się niewykonalne wraz ze wzrostem liczby spinów, ponieważ wymiar przestrzeni Hilberta rośnie wykładniczo jako $2^N$ dla $N$ spinów. To czyni go naturalnym kandydatem do symulacji kwantowej.

Variational Quantum Eigensolver (VQE) to hybrydowy algorytm kwantowo-klasyczny zaprojektowany do estymacji energii stanu podstawowego Hamiltonianu. Działa poprzez przygotowanie sparametryzowanego stanu kwantowego $|\psi(\theta)\rangle$ (zwanego ansatzem) na komputerze kwantowym i pomiar wartości oczekiwanej $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Klasyczny optymalizator następnie iteracyjnie dostosowuje parametry $\theta$ w celu minimalizacji tej energii, korzystając z zasady wariacyjnej, która gwarantuje, że mierzona energia jest zawsze górnym ograniczeniem prawdziwej energii stanu podstawowego.

W tym samouczku używamy ansatza efficient_su2 z biblioteki układów Qiskit, który konstruuje warstwy jednobitowych rotacji i bramek splątujących. Optymalizacja jest przeprowadzana przy użyciu algorytmu Simultaneous Perturbation Stochastic Approximation (SPSA), który jest dobrze przystosowany do zaszumionego sprzętu kwantowego, ponieważ szacuje gradienty przy użyciu tylko dwóch wywołań funkcji na iterację, niezależnie od liczby parametrów.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane następujące pakiety:

• Qiskit SDK w wersji 2.0 lub nowszej, z obsługą wizualizacji
• Qiskit Runtime w wersji 0.44 lub nowszej (pip install qiskit-ibm-runtime)

Konfiguracja

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

Przykład małej skali

W tej sekcji przechodzimy przez każdy krok wzorca Qiskit w małej skali, wyjaśniając kluczowe komponenty w miarę budowania przepływu pracy.

Krok 1: Mapowanie klasycznych danych wejściowych na problem kwantowy

• Wejście: Liczba spinów
• Wyjście: Ansatz i Hamiltonian modelujące łańcuch Heisenberga

Skonstruuj ansatz i Hamiltonian modelujące 10-spinowy łańcuch Heisenberga. W tym kroku zbudujemy 10-spinowy Hamiltonian Heisenberga na mapie sprzężeń najmniej zajętego backendu i przygotujemy 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: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

• Wejście: Abstrakcyjny Circuit, obserwowalny
• Wyjście: Docelowy Circuit i obserwowalny, zoptymalizowane dla wybranego QPU

Użyj funkcji generate_preset_pass_manager z Qiskit, aby automatycznie wygenerować procedurę optymalizacji naszego Circuit względem wybranego QPU. Wybieramy optimization_level=3, co zapewnia najwyższy poziom optymalizacji spośród wstępnie zdefiniowanych menedżerów przejść. Dodajemy również przejścia harmonogramowania ALAPScheduleAnalysis i PadDynamicalDecoupling, aby tłumić błędy dekoherencji.

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: Wykonanie przy użyciu prymitywów Qiskit

• Wejście: Docelowy Circuit i obserwowalny
• Wyjście: Wyniki optymalizacji

Minimalizuj szacowaną energię stanu podstawowego systemu poprzez optymalizację parametrów Circuit. Użyj prymitywu Estimator z Qiskit Runtime do obliczania funkcji kosztu podczas optymalizacji.

Ponieważ zoptymalizowaliśmy Circuit dla backendu w Kroku 2, możemy uniknąć transpilacji na serwerze Runtime, ustawiając skip_transpilation=True i przekazując zoptymalizowany Circuit. W tej demonstracji uruchomimy kod na QPU, używając prymitywów qiskit-ibm-runtime. Aby uruchomić z prymitywami opartymi na wektorach stanów qiskit, zastąp blok kodu używający prymitywów Qiskit Runtime skomentowanym blokiem.

W tym samouczku używamy Simultaneous Perturbation Stochastic Approximation (SPSA), czyli optymalizatora opartego na gradiencie. Poniżej przedstawiamy jego krótkie wprowadzenie oraz kod implementujący SPSA przy użyciu Qiskit v2.0.

Wprowadzenie do SPSA

Simultaneous Perturbation Stochastic Approximation (SPSA) [1] to algorytm optymalizacji, który przybliża cały wektor gradientu przy użyciu tylko dwóch wywołań funkcji na każdej iteracji. Niech $f:\mathbb{R}^p\rightarrow \mathbb{R}$ będzie funkcją kosztu z $p$ parametrami do optymalizacji, a $x_i\in \mathbb{R}^p$ wektorem parametrów na $i$-tym kroku iteracji. Aby obliczyć gradient, tworzony jest losowy wektor $\Delta_i$ o rozmiarze $p$, gdzie każdy element $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, jest próbkowany jednostajnie z $\{-1, 1\}$. Następnie każdy element losowego wektora $\Delta_i$ jest mnożony przez małą wartość $c_i$, tworząc losowe zaburzenie. Gradient jest wówczas szacowany 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}}.$

Intuicyjnie, ponieważ podczas estymacji gradientu stosowane jest losowe zaburzenie, oczekuje się, że małe odchylenia w dokładnych wartościach $f$ wynikające z szumu mogą być tolerowane i uwzględniane. W rzeczywistości SPSA jest szczególnie znane ze swojej odporności na szum i wymaga tylko dwóch wywołań sprzętowych na każdą iterację. Jest zatem jednym z najbardziej preferowanych optymalizatorów do implementacji algorytmów wariacyjnych.

W tym samouczku hiperparametry dla $i$-tej iteracji, $a_i$ i $c_i$, są obliczane jako

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

gdzie stałe wartości przyjęto jako $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ i $\gamma = 0.4$. Wartości te zaczerpnięto z [2]. Odpowiednie dostrojenie hiperparametrów jest niezbędne do uzyskania dobrej wydajności 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
        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)

Tutaj ustawiamy maxiter = 50. Zauważ, że ponieważ każda iteracja wymaga dwóch wywołań funkcji do obliczenia gradientu, łączna liczba wywołań funkcji wyniesie $2 \times \text{maxiter}$. Wartość maxiter można zwiększyć do dowolnej wyższej wartości, aby uzyskać lepszą estymację energii.

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

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

Krok 4: Przetwarzanie końcowe i zwrócenie wyniku w żądanym formacie klasycznym

• Wejście: Estymacje energii stanu podstawowego podczas optymalizacji
• Wyjście: Szacowana energia stanu podstawowego

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

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Przykład sprzętowy dużej skali

Przykład sprzętowy dużej skali nie jest uwzględniony w tym samouczku. Wraz ze wzrostem liczby kubitów VQE napotyka istotne trudności z powodu zjawiska jałowego plateau: gradient funkcji kosztu zanika wykładniczo wraz z rozmiarem układu, co sprawia, że optymalizacja staje się praktycznie niewykonalna dla dużych układów. W połączeniu z szumem sprzętowym oznacza to, że skalowanie VQE do większych łańcuchów spinów nie daje wiarygodnych i powtarzalnych wyników. Podejścia pozwalające przezwyciężyć te ograniczenia znajdziesz w sekcji Kolejne kroki poniżej.

Wyzwanie

Teraz, gdy masz działającą implementację VQE dla łańcucha Heisenberga, spróbuj wykonać następujące zadania:

1. Eksperymentuj z głębokością ansatza: Zmodyfikuj parametr reps w efficient_su2 (na przykład wypróbuj reps=1 i reps=3). Jak głębokość ansatza wpływa na szacowaną energię stanu podstawowego i szybkość zbieżności? W którym miejscu obserwujesz malejące zyski lub niestabilność?
2. Dostrajaj hiperparametry SPSA: Dostosuj parametry harmonogramu współczynnika uczenia (a, c, alpha, gamma, A) i obserwuj, jak wpływają na zbieżność. Czy możesz znaleźć konfigurację, która zbiega szybciej niż domyślne wartości użyte tutaj?
3. Porównaj topologie sprzężeń: Zamiast używać natywnej mapy sprzężeń backendu, spróbuj skonstruować prosty łańcuch liniowy z połączeniami najbliższych sąsiadów i porównaj wyniki. Jak łączność fizycznego sprzętu wpływa na głębokość przetranspilowanego układu i końcową estymację energii?

Odniesienia

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

Kolejne kroki

Rekomendacje

Jeśli ta praca wydała ci się interesująca, możesz zainteresować się następującymi materiałami:

• Wypróbuj Sample-based Quantum Diagonalization (SQD): Jak pokazano w tym samouczku, VQE napotyka trudności przy skalowaniu ze względu na jałowe plateau i duże zapotrzebowanie na pomiary. IBM opracował Sample-based Quantum Diagonalization (SQD) jako bardziej skalowalną alternatywę. W przeciwieństwie do VQE, SQD całkowicie unika optymalizacji wariacyjnej; zamiast tego komputer kwantowy generuje próbki, a komputer klasyczny rzutuje Hamiltonian na podprzestrzeń rozpiętą przez te próbki i go diagonalizuje. Zapewnia to górne ograniczenie energii stanu podstawowego przy znacznie mniejszej liczbie pomiarów i bez podatności na jałowe plateau. Skorzystaj z samouczka SQD, aby zobaczyć to podejście w działaniu.
• Odkryj kurs Quantum Diagonalization Algorithms: Pogłęb swoje rozumienie zarówno VQE, jak i SQD, w tym ich kompromisów, w kursie Algorytmy kwantowej diagonalizacji na IBM Quantum Learning.