Grundzustandsenergie-Schätzung der Heisenberg-Kette mit VQE

Geschätzte Laufzeit: 37 Minuten auf einem Heron-Prozessor (HINWEIS: Dies ist lediglich eine Schätzung. (Die Laufzeit kann variieren.)

Lernziele

• Wie man eine Heisenberg-Spinkette mithilfe von Qiskit als Quantenhamiltonian modelliert
• Wie man den SPSA-Optimierer zur Abschätzung der Grundzustandsenergie eines Quantensystems einsetzt
• So führen Sie Variations-Workflows auf der „ IBM® “-Quantenhardware mithilfe von „ IBM Quantum “-Primitiven und -Sitzungen aus

Voraussetzungen

• Grundlagen der Quanteninformation
• Einführung in Qiskit-Muster
• Entwurf eines Variationsalgorithmus

Hintergrund

Die Heisenberg-Spinkette ist eines der am intensivsten untersuchten Modelle in der Festkörperphysik und im Bereich des Quantenmagnetismus. Es beschreibt ein eindimensionales Gitter aus wechselwirkenden Quantenspins, in dem benachbarte Spins durch Austauschwechselwirkungen gekoppelt sind. Der Hamilton-Operator für das isotrope Heisenberg-Modell mit einem externen Magnetfeld lautet:

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

wobei $X_i$, $Y_i$ und $Z_i$ die auf die Position $i$ wirkenden Pauli-Operatoren sind, die Summe $\langle i,j \rangle$ über die Paare der nächsten Nachbarn läuft, $J_x = J_y = J_z = 0.5$ die Austauschkopplungskonstanten sind (in diesem Tutorial isotrop) und $h_i$ ein ortsabhängiges externes Magnetfeld darstellt. In diesem Tutorial werden die Magnetfeldwerte zufällig aus dem Intervall $[-1, 1]$ entnommen. Beachten Sie, dass in der folgenden Implementierung die Menge der „Nearest-Neighbor“-Paare durch die native Kopplung des Hardware-Backends zwischen den ersten $N$ Qubits bestimmt wird, die je nach Gerätetopologie möglicherweise keine streng lineare Kette bilden.

Das Verständnis der Grundzustandsenergie dieses Hamilton-Operators ist in der Physik von grundlegender Bedeutung. Der Grundzustand enthält Informationen über Quantenphasenübergänge, die Verschränkungsstruktur und die magnetische Ordnung. Im klassischen Sinne wird die Berechnung der exakten Grundzustandsenergie mit steigender Anzahl von Spins unlösbar, da die Dimension des Hilbertraums exponentiell mit $2^N$ für $N$ Spins wächst. Damit ist es ein naheliegender Kandidat für die Quantensimulation.

Der Variational Quantum Eigensolver (VQE) ist ein hybrid-quantentheoretisch-klassischer Algorithmus, der zur Schätzung der Grundzustandsenergie eines Hamilton-Operators entwickelt wurde. Dabei wird auf einem Quantencomputer ein parametrisierter Quantenzustand $|\psi(\theta)\rangle$ (ein sogenannter Ansatz) erzeugt und der Erwartungswert $\langle \psi(\theta) | H | \psi(\theta) \rangle$ gemessen. Ein klassischer Optimierer passt dann die Parameter $\theta$ iterativ an, um diese Energie zu minimieren, wobei er sich das Variationsprinzip zunutze macht, das garantiert, dass die gemessene Energie stets eine Obergrenze für die tatsächliche Grundzustandsenergie darstellt.

In diesem Tutorial verwenden wir den efficient_su2 Ansatz aus der Schaltkreisbibliothek von Qiskit, der Schichten aus Rotationen einzelner Quantenbits und Verschränkungsgates aufbaut. Die Optimierung erfolgt mithilfe des SPSA-Algorithmus (Simultaneous Perturbation Stochastic Approximation), der sich gut für verrauschte Quantenhardware eignet, da er Gradienten mit nur zwei Funktionsauswertungen pro Iteration schätzt, unabhängig von der Anzahl der Parameter.

Anforderungen

Vergewissern Sie sich vor Beginn dieses Tutorials, dass Sie die folgenden Komponenten installiert haben:

• Qiskit SDK v2.0 oder höher, mit Visual isierungsunterstützung
• Qiskit Runtime v0.44 oder später (pip install qiskit-ibm-runtime)

Einrichtung

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

Ein kleines Beispiel

In diesem Abschnitt gehen wir jeden Schritt des Qiskit-Musters im Kleinen durch und erläutern dabei die wichtigsten Komponenten, während wir den Arbeitsablauf aufbauen.

Schritt 1: Klassische Eingaben einem Quantenproblem zuordnen

• Eingabe: Anzahl der Spins
• Ausgabe: Ansatz und Hamiltonian zur Modellierung der Heisenberg-Kette

Entwickeln Sie einen Ansatz und einen Hamilton-Operator, die eine Heisenberg-Kette mit 10 Spins modellieren. In diesem Schritt werden wir einen Heisenberg-Hamilton-Operator mit 10 Spins über der Kopplungskarte des am wenigsten ausgelasteten Backends aufbauen und den efficient_su2 Ansatz vorbereiten.

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")

Output:

Schritt 2: Optimierung des Problems für die Ausführung auf Quantenhardware

• Eingabe: Abstrakter Schaltkreis, beobachtbar
• Ausgabe: Zielschaltung und Observable, optimiert für die ausgewählte QPU

Verwenden Sie die Funktion generate_preset_pass_manager von Qiskit, um automatisch eine Optimierungsroutine für unsere Schaltung in Bezug auf die ausgewählte QPU zu erstellen. Wir haben uns für optimization_level=3 entschieden, das den höchsten Optimierungsgrad der voreingestellten Passmanager bietet. Zur Unterdrückung von Dekohärenzfehlern werden auch ALAPScheduleAnalysis und PadDynamicalDecoupling in die Planung einbezogen.

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)

Output:

Schritt 3: Mit „ Qiskit primitives “ ausführen

• Eingabe: Zielschaltung und Beobachtungsobjekt
• Ausgabe: Ergebnisse der Optimierung

Minimieren Sie die geschätzte Grundzustandsenergie des Systems durch Optimierung der Schaltungsparameter. Verwenden Sie die Primitive Estimator „ IBM Quantum “, um die Kostenfunktion während der Optimierung zu berechnen.

Da wir den Schaltkreis in Schritt 2 für das Backend optimiert haben, können wir die Transpilation auf dem Quantum Compute-Server vermeiden, indem wir festlegen skip_transpilation=True und den optimierten Schaltkreis übergeben. Für diese Demo werden wir die Berechnung auf einer QPU unter Verwendung von „ IBM Quantum “-Primitiven durchführen. Um mit auf Zustandsvektoren basierenden qiskit Primitiven zu arbeiten, ersetzen Sie den Codeabschnitt, der „ IBM Quantum “-Primitiven verwendet, durch den auskommentierten Abschnitt.

In diesem Tutorial verwenden wir die Simultaneous Perturbation Stochastic Approximation (SPSA), einen gradientenbasierten Optimierer. Als Nächstes geben wir eine kurze Einführung in das Thema und stellen den Code zur Implementierung von SPSA mit Qiskit v2.0 zur Verfügung.

Wir stellen SPSA vor

Die simultane stochastische Approximation durch Störung (SPSA) [1] ist ein Optimierungsalgorithmus, der den gesamten Gradientenvektor mit nur zwei Funktionsaufrufen pro Iteration approximiert. Sei $f:\mathbb{R}^p\rightarrow \mathbb{R}$ die Kostenfunktion mit den zu optimierenden Parametern $p$ und sei $x_i\in \mathbb{R}^p$ der Parametervektor im Schritt $i^{th}$ der Iteration. Zur Berechnung des Gradienten wird ein Zufallsvektor $\Delta_i$ der Größe $p$ erzeugt, wobei jedes Element $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$ gleichverteilten Werten aus $\{-1, 1\}$ entnommen wird. Anschließend wird jedes Element des Zufallsvektors $\Delta_i$ mit einem kleinen Wert $c_i$ multipliziert, um eine zufällige Störung zu erzeugen. Der Gradient wird dann wie folgt geschätzt:

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

Da bei der Gradientenschätzung eine zufällige Störung einwirkt, ist intuitiv zu erwarten, dass kleine Abweichungen der exakten Werte von $f$, die auf Rauschen zurückzuführen sind, toleriert und berücksichtigt werden können. Tatsächlich ist SPSA dafür bekannt, besonders unempfindlich gegenüber Störsignalen zu sein, und benötigt für jede Iteration nur zwei Hardware-Aufrufe. Daher ist er einer der am meisten bevorzugten Optimierer für die Implementierung von Variationsalgorithmen.

In diesem Tutorial werden die Hyperparameter für die Iteration „ $i^{th}$ “, „ $a_i$ “ und „ $c_i$ “ wie folgt berechnet:

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

wobei die konstanten Werte wie folgt angenommen werden: $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ und $\gamma = 0.4$. Diese Werte wurden aus [2] entnommen. Eine angemessene Abstimmung der Hyperparameter ist erforderlich, um eine gute Leistung aus SPSA herauszuholen.

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)

Hier legen wir die maxiter = 50. Beachten Sie, dass jede Iteration zwei Funktionsaufrufe zur Berechnung des Gradienten erfordert, sodass die Gesamtzahl der Funktionsaufrufe $2 \times \text{maxiter}$ beträgt. Der Wert maxiter kann für eine genauere Energieschätzung auf einen beliebig höheren Wert erhöht werden.

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

Output:

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

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

• Eingabe: Schätzungen der Grundzustandsenergie während der Optimierung
• Ergebnis: Geschätzte Grundzustandsenergie

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

Output:

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Output:

Beispiel für groß angelegte Hardware

Ein Beispiel für eine groß angelegte Hardware-Lösung ist in diesem Tutorial nicht enthalten. Mit steigender Anzahl von Qubits steht die VQE aufgrund des karge Hochebene folgenden Phänomens vor erheblichen Herausforderungen: Der Gradient der Kostenfunktion nimmt mit zunehmender Systemgröße exponentiell ab, was die Optimierung bei großen Schaltkreisen praktisch unmöglich macht. In Verbindung mit Hardware-Rauschen bedeutet dies, dass die Skalierung von VQE auf längere Spinketten keine zuverlässig reproduzierbaren Ergebnisse liefert. Ansätze zur Überwindung dieser Einschränkungen finden Sie im Abschnitt „Nächste Schritte“ weiter unten.

Herausforderung

Nachdem Sie nun eine funktionierende VQE-Implementierung für die Heisenberg-Kette haben, versuchen Sie Folgendes:

1. reps=3Experimentieren Sie mit der Ansatz-Tiefe: Ändern Sie den reps Parameter in efficient_su2 (probieren reps=1 Sie beispielsweise und aus). Wie wirkt sich die Ansatztiefe auf die geschätzte Grundzustandsenergie und die Konvergenzgeschwindigkeit aus? Ab wann stellen Sie abnehmende Erträge oder Instabilität fest?
2. SPSA-Hyperparameter optimieren: Passen Sie die Parameter des Lernratenverlaufs (a, c, alpha, gamma, A) an und beobachten Sie, wie sich dies auf die Konvergenz auswirkt. Findest du eine Konfiguration, die schneller konvergiert als die hier verwendeten Standardeinstellungen?
3. Vergleichen Sie Kopplungstopologien: Anstatt die native Kopplungskarte des Backends zu verwenden, versuchen Sie, eine einfache lineare Kette nach dem Prinzip der nächsten Nachbarn zu erstellen, und vergleichen Sie die Ergebnisse. Inwiefern wirkt sich die Konnektivität der physikalischen Hardware auf die Tiefe der transpilierten Schaltung und die endgültige Energieschätzung aus?

Referenzen

[1] Spall, J. C. (2002). Implementierung des Algorithmus der simultanen Störung für die stochastische Optimierung. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817–823.

[2] Sahin, M. Emre u. a. (2025). Qiskit- Machine Learning : eine Open-Source-Bibliothek für groß angelegte Aufgaben im Bereich des quantenbasierten maschinellen Lernens auf Quantenhardware und klassischen Simulatoren. arXiv:2505.17756.

Nächste Schritte