Grundzustandsenergieschätzung der Heisenberg-Kette mit VQE

Geschätzter Ressourcenverbrauch: 37 Minuten auf einem Heron-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann abweichen.)

Lernziele

Nach Abschluss dieses Tutorials kannst du folgende Informationen verstehen:

• Wie man eine Heisenberg-Spin-Kette als Quanten-Hamiltonian mit Qiskit modelliert
• Wie man den SPSA-Optimierer verwendet, um die Grundzustandsenergie eines Quantensystems zu schätzen
• Wie man variationelle Workflows auf IBM®-Quantenhardware mit Qiskit Runtime-Primitiven und Sessions ausführt

Voraussetzungen

Es wird empfohlen, dass du dich mit diesen Themen vertraut machst:

• Grundlagen der Quanteninformation
• Einführung in Qiskit-Patterns
• Design variationeller Algorithmen

Hintergrund

Die Heisenberg-Spin-Kette ist eines der am weitesten erforschten Modelle in der Physik der kondensierten Materie und des Quantenmagnetismus. Sie beschreibt ein eindimensionales Gitter wechselwirkender Quantenspins, bei dem benachbarte Spins durch Austauschwechselwirkungen gekoppelt sind. Der Hamiltonian 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 Pauli-Operatoren sind, die auf Gitterplatz $i$ wirken, die Summe $\langle i,j \rangle$ über nächste-Nachbar-Paare 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]$ gesampelt. Beachte, dass in der folgenden Implementierung die Menge der „nächsten Nachbarn" durch die native Kopplung des Hardware-Backends unter den ersten $N$ Qubits bestimmt wird, die je nach Gerätetopologie nicht unbedingt eine strenge lineare Kette bildet.

Das Verständnis der Grundzustandsenergie dieses Hamiltonians ist von grundlegender Bedeutung in der Physik. Der Grundzustand codiert Informationen über Quantenphasenübergänge, Verschränkungsstruktur und magnetische Ordnung. Klassisch wird die Berechnung der exakten Grundzustandsenergie mit zunehmender Anzahl von Spins unlösbar, da die Hilbertraum-Dimension exponentiell als $2^N$ für $N$ Spins skaliert. Das macht sie zu einem natürlichen Kandidaten für die Quantensimulation.

Der Variational Quantum Eigensolver (VQE) ist ein hybrider Quanten-Klassik-Algorithmus zur Schätzung der Grundzustandsenergie eines Hamiltonians. Er bereitet einen parametrisierten Quantenzustand $|\psi(\theta)\rangle$ (Ansatz genannt) auf einem Quantencomputer vor und misst den Erwartungswert $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Ein klassischer Optimierer passt dann iterativ die Parameter $\theta$ an, um diese Energie zu minimieren, wobei das Variationsprinzip genutzt wird, das garantiert, dass die gemessene Energie stets eine obere Schranke zur wahren Grundzustandsenergie ist.

In diesem Tutorial verwenden wir den efficient_su2-Ansatz aus Qiskits Circuit-Bibliothek, der Schichten von Ein-Qubit-Rotationen und Verschränkungs-Gates konstruiert. Die Optimierung wird mit dem Simultaneous Perturbation Stochastic Approximation (SPSA)-Algorithmus durchgeführt, der sich gut für verrauschte Quantenhardware eignet, weil er Gradienten mit nur zwei Funktionsauswertungen pro Iteration unabhängig von der Parameteranzahl schätzt.

Voraussetzungen

Bevor du mit diesem Tutorial beginnst, stelle sicher, dass du Folgendes installiert hast:

• Qiskit SDK v2.0 oder höher, mit Unterstützung für Visualisierung
• Qiskit Runtime v0.44 oder höher (pip install qiskit-ibm-runtime)

Setup

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

Kleinskaliges Beispiel

In diesem Abschnitt gehen wir jeden Schritt des Qiskit-Patterns in kleinem Maßstab durch und erläutern die wichtigsten Komponenten, während wir den Workflow aufbauen.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

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

Konstruiere einen Ansatz und einen Hamiltonian, die eine 10-Spin-Heisenberg-Kette modellieren. In diesem Schritt erstellen wir einen 10-Spin-Heisenberg-Hamiltonian über die Kopplungskarte des am wenigsten ausgelasteten Backends und bereiten den efficient_su2-Ansatz vor.

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

Schritt 2: Problem für die Ausführung auf Quantenhardware optimieren

• Eingabe: Abstrakter Circuit, Observable
• Ausgabe: Ziel-Circuit und Observable, optimiert für den ausgewählten QPU

Verwende die Funktion generate_preset_pass_manager aus Qiskit, um automatisch eine Optimierungsroutine für unseren Circuit bezüglich des ausgewählten QPU zu erzeugen. Wir wählen optimization_level=3, das die höchste Optimierungsstufe der vordefinierten Pass-Manager bietet. Außerdem fügen wir die Scheduling-Passes ALAPScheduleAnalysis und PadDynamicalDecoupling ein, um Dekohärenzfehler zu unterdrücken.

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)

Schritt 3: Ausführung mit Qiskit-Primitiven

• Eingabe: Ziel-Circuit und Observable
• Ausgabe: Ergebnisse der Optimierung

Minimiere die geschätzte Grundzustandsenergie des Systems, indem du die Circuit-Parameter optimierst. Verwende das Estimator-Primitiv aus Qiskit Runtime, um die Kostenfunktion während der Optimierung auszuwerten.

Da wir den Circuit für das Backend in Schritt 2 optimiert haben, können wir die Transpilation auf dem Runtime-Server vermeiden, indem wir skip_transpilation=True setzen und den optimierten Circuit übergeben. Für diese Demo führen wir die Berechnungen auf einem QPU mit qiskit-ibm-runtime-Primitiven aus. Um stattdessen qiskit-Zustandsvektor-basierte Primitive zu verwenden, ersetze den Codeblock mit den Qiskit-IBM-Runtime-Primitiven durch den auskommentierten Block.

In diesem Tutorial verwenden wir Simultaneous Perturbation Stochastic Approximation (SPSA), einen gradientenbasierten Optimierer. Im Folgenden geben wir eine kurze Einführung dazu und stellen den Code zur Implementierung von SPSA mit Qiskit v2.0 bereit.

Einführung in SPSA

Simultaneous Perturbation Stochastic Approximation (SPSA) [1] ist ein Optimierungsalgorithmus, der den gesamten Gradientenvektor mit nur zwei Funktionsauswertungen pro Iteration approximiert. Sei $f:\mathbb{R}^p\rightarrow \mathbb{R}$ die Kostenfunktion mit $p$ zu optimierenden Parametern und $x_i\in \mathbb{R}^p$ der Parametervektor beim $i^{th}$ Schritt der Iteration. Zur Berechnung des Gradienten wird ein zufälliger Vektor $\Delta_i$ der Größe $p$ erstellt, wobei jedes Element $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, gleichmäßig aus $\{-1, 1\}$ gesampelt wird. Anschließend wird jedes Element des zufälligen Vektors $\Delta_i$ mit einem kleinen Wert $c_i$ multipliziert, um eine zufällige Störung zu erzeugen. Der Gradient wird dann geschätzt als

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

Intuitiv gilt: Da während der Gradientenschätzung eine zufällige Störung angewendet wird, können kleine Abweichungen in den exakten Werten von $f$ aufgrund von Rauschen toleriert und berücksichtigt werden. Tatsächlich ist SPSA besonders bekannt für seine Robustheit gegenüber Rauschen und benötigt nur zwei Hardware-Aufrufe pro Iteration. Er ist daher einer der bevorzugten Optimierer für die Implementierung variationeller Algorithmen.

In diesem Tutorial werden die Hyperparameter für den $i^{th}$ Schritt, $a_i$ und $c_i$, berechnet als

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

wobei die Konstantenwerte $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ und $\gamma = 0.4$ sind. Diese Werte wurden aus [2] entnommen. Eine angemessene Abstimmung der Hyperparameter ist notwendig, um eine gute Leistung von SPSA zu erzielen.

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 setzen wir maxiter = 50. Beachte, dass da jede Iteration zwei Funktionsaufrufe zur Berechnung des Gradienten benötigt, die Gesamtzahl der Funktionsaufrufe $2 \times \text{maxiter}$ beträgt. maxiter kann auf einen höheren Wert erhöht werden, um eine bessere Energieschätzung zu erzielen.

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

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

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

• Eingabe: Grundzustandsenergieschätzungen während der Optimierung
• Ausgabe: Geschätzte Grundzustandsenergie

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

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Großskaliges Hardware-Beispiel

Ein großskaliges Hardware-Beispiel ist in diesem Tutorial nicht enthalten. Mit zunehmender Qubit-Anzahl stößt VQE auf erhebliche Herausforderungen durch das Barren-Plateau-Phänomen: der Gradient der Kostenfunktion nimmt exponentiell mit der Systemgröße ab, wodurch die Optimierung für große Schaltkreise praktisch undurchführbar wird. In Kombination mit Hardware-Rauschen bedeutet dies, dass die Skalierung von VQE auf größere Spin-Ketten keine zuverlässig reproduzierbaren Ergebnisse liefert. Für Ansätze, die diese Einschränkungen überwinden, sieh den Abschnitt „Nächste Schritte" unten.

Challenge

Jetzt, da du eine funktionierende VQE-Implementierung für die Heisenberg-Kette hast, probiere Folgendes aus:

1. Experimentiere mit der Ansatz-Tiefe: Ändere den reps-Parameter in efficient_su2 (probiere zum Beispiel reps=1 und reps=3). Wie wirkt sich die Ansatz-Tiefe auf die geschätzte Grundzustandsenergie und Konvergenzgeschwindigkeit aus? Ab welchem Punkt beobachtest du abnehmende Erträge oder Instabilität?
2. Optimiere SPSA-Hyperparameter: Passe die Lernraten-Planungsparameter (a, c, alpha, gamma, A) an und beobachte, wie sie die Konvergenz beeinflussen. Kannst du eine Konfiguration finden, die schneller konvergiert als die hier verwendeten Standardwerte?
3. Vergleiche Kopplungstopologien: Statt die native Kopplungskarte des Backends zu verwenden, versuche eine einfache nächste-Nachbar-Linearkette zu konstruieren und vergleiche die Ergebnisse. Wie wirkt sich die Konnektivität der physischen Hardware auf die transpilierte Schaltkreistiefe und die endgültige Energieschätzung aus?

Referenzen

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

Nächste Schritte

Empfehlungen

Wenn du diese Arbeit interessant fandest, könnte dich folgendes Material interessieren:

• Sample-based Quantum Diagonalization (SQD) ausprobieren: Wie in diesem Tutorial gezeigt, hat VQE im großen Maßstab Herausforderungen durch Barren Plateaus und hohen Messaufwand. IBM hat Sample-based Quantum Diagonalization (SQD) als skalierbarere Alternative entwickelt. Im Gegensatz zu VQE vermeidet SQD die variationelle Optimierung vollständig; stattdessen generiert ein Quantencomputer Proben und ein klassischer Computer projiziert den Hamiltonian auf einen Unterraum, der durch diese Proben aufgespannt wird, und diagonalisiert ihn. Dies liefert eine obere Schranke für die Grundzustandsenergie mit deutlich weniger Messungen und ohne Anfälligkeit für Barren Plateaus. Folge dem SQD-Tutorial, um diesen Ansatz in der Praxis zu sehen.
• Den Kurs Quantum Diagonalization Algorithms erkunden: Vertiefe dein Verständnis von VQE und SQD, einschließlich ihrer Abwägungen, im Kurs Quantum diagonalization algorithms auf IBM Quantum Learning.