Mehrproduktformeln zur Reduzierung von Trotter-Fehlern

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

Lernziele

• Wie Multi-Product-Formeln (MPFs) den Trotter-Fehler bei der Hamilton-Simulation durch die Kombination von Erwartungswerten aus mehreren flachen Schaltkreisen reduzieren
• Wann MPFs gegenüber Standardproduktformeln Vorteile bieten und wann sie nicht das richtige Mittel der Wahl sind
• So berechnen Sie statische und dynamische MPF-Koeffizienten mithilfe des qiskit_addon_mpf Pakets
• So führen Sie einen MPF-Workflow von Anfang bis Ende auf einer „ IBM Quantum® “-Hardware durch, einschließlich Transpilation, Fehlerbehebung und Nachbearbeitung

Voraussetzungen

• Kompilierungsmethoden für Hamilton-Simulationsschaltungen – Einführung in Trotter-Schaltungen (Produktformel) in Qiskit.
• Produktformeln in Qiskit, insbesondere die SuzukiTrotter Synthese-Klassen und LieTrotter .
• Qiskit primitives und die Estimator-Schnittstelle.

Hintergrund

Was sind Mehrproduktformeln?

Bei der Simulation von Quantensystemen auf einem Quantencomputer besteht eine zentrale Aufgabe darin, den Zeitentwicklungsoperator $e^{-iHt}$ für einen Hamilton-Operator $H$ zu approximieren. Der Standardansatz nutzt Produktformeln (PFs), die auch als Trotter-Suzuki-Zerlegungen bekannt sind. Diese zerlegen $H = \sum_{a=1}^d F_a$ in Terme, deren einzelne Unitarien $e^{-iF_a t}$ effizient zu implementieren sind, und approximieren dann die vollständige Entwicklung als geordnetes Produkt dieser einfacheren Unitarien.

Die Produktformel erster Ordnung (Lie-Trotter) lautet:

was zu einem quadratischen Fehler führt: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. Symmetrische Formeln höherer Ordnung $S_{2\chi}(t)$, wobei $\chi$ die Ordnung der symmetrischen Produktformel bezeichnet (siehe Ref. [1] ), konvergieren schneller als $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, allerdings auf Kosten tieferer Schaltkreise pro Schritt.

Um den Fehler bei einem $\chi$ fester Ordnung zu verringern, wird die gesamte Evolutionszeit $t$ üblicherweise in $k$ kleinere Trotter-Schritte unterteilt. Jeder Schritt approximiert $e^{-iHt/k}$ durch eine Produktformel, und die Schritte werden aneinandergereiht:

Bei einer symmetrischen Formel d $2\chi$ -ter Ordnung skaliert der verbleibende Trotter-Fehler dann wie folgt: $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Eine Erhöhung von $k$ führt also zu einer raschen Unterdrückung des Trotter-Fehlers – allerdings wird dadurch auch die Schaltung linear tiefer, was auf rauschbehafteter Hardware zu mehr kumuliertem Gatterrauschen führt. Genau diese Spannung zwischen dem Trotter-Fehler (begünstigt größere $k$ ) und dem Hardware-Rauschen (begünstigt kleinere $k$ ) soll durch Multiproduktformeln aufgelöst werden. Beachten Sie, dass es bei MPFs darum geht, Ergebnisse aus verschiedenen Auswahlen von „ $k$ “ in einer festgelegten Reihenfolge $\chi$ zu kombinieren – sie ändern die Reihenfolge der zugrunde liegenden Produktformel nicht.

Multi-Produkt-Formeln (MPFs) [1] bilden eine gewichtete lineare Kombination von Erwartungswerten, die aus mehreren flacheren Trotter-Schaltkreisen gewonnen werden, wobei jeder eine unterschiedliche Anzahl von Trotter-Schritten verwendet $k_1, k_2, \ldots, k_r$ (eine Reihe von $r$ Schrittzahlen):

wobei $\langle A \rangle_{k_j}(t)$ der Erwartungswert einer Observablen $A$ zum Zeitpunkt $t$ ist, der anhand eines Trotter-Kreises mit $k_j$ Schritten geschätzt wurde, und die Koeffizienten $\{x_j\}_{j=1}^r$ so gewählt sind, dass sich die führenden Trotter-Fehlerterme in der Kombination aufheben. Wir werden in Schritt 4 auf diesen Ausdruck zurückkommen, wo wir ihn explizit auswerten, um unsere Trotter-Ergebnisse zusammenzufassen. Der entscheidende praktische Aspekt ist, dass die tiefste Schaltung im MPF nur $k_{\max}$ Schritte benötigt, was deutlich weniger ist als der einzelne Schritt $k$, der erforderlich wäre, um denselben effektiven Trotter-Fehler direkt zu erreichen. Durch die flacheren Schaltkreise eignet sich der MPF-Ansatz besser für rauschbehaftete Hardware.

Wie werden die Koeffizienten ermittelt?

Es gibt zwei Familien von MPF-Koeffizienten:

Statische Koeffizienten sind unabhängig vom Hamilton-Operator, vom Anfangszustand und von der Entwicklungszeit. Sie werden durch die Lösung eines linearen Gleichungssystems $Ax = b$ ermittelt, das die Aufhebung der führenden Trotter-Fehlerterme gewährleistet. Für eine Reihe von Trotter-Schritten $\{k_j\}_{j=1}^r$ in Verbindung mit einer symmetrischen Produktformel der Ordnung $2\chi$ führt die Entwicklung des Trotter-Fehlers in umgekehrten Potenzen von $k_j$ zu Nebenbedingungen der Form:

wobei die ganzzahligen Exponenten $\{\eta_n\}$ die Ordnungen der aufeinanderfolgenden Trotter-Fehlerterme für die gewählte Produktformel sind. Für eine symmetrische PF der Ordnung $2\chi$ gilt, dass der Hauptfehler in $\left[S_{2\chi}(t/k)\right]^k$ wie folgt skaliert: $1/k^{2\chi}$, mit nachfolgenden Korrekturen gemäß $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ – die Exponenten lauten also $\eta_n = 2\chi + 2n$. Bei nicht-symmetrischen PFs tragen sowohl ungerade als auch gerade Potenzen bei, und es gilt $\eta_n = 2\chi + n$. Die vollständige Herleitung findet sich in Ref. [1]. Die erste Gleichung im obigen System gewährleistet die Unverzerrtheit (die MPF gibt den exakten Erwartungswert im Grenzfall „ $k_j \to \infty$ “ wieder), und die übrigen $r-1$ Gleichungen heben nacheinander die ersten $r-1$ Trotter-Fehlerterme auf. Wenn die resultierende $L_1$ -Norm $\|x\|_1$ zu groß ist (was das Abtastrauschen verstärkt), können Sie stattdessen eine approximative Optimierung lösen, die $\|x\|_1$ begrenzt und gleichzeitig $\|Ax - b\|$ minimiert.

Die dynamischen Koeffizienten [2], [3] hängen zudem vom Hamilton-Operator, vom Anfangszustand und von der Entwicklungszeit $t$ ab. Sie minimieren den Abstand in der Frobenius-Norm zwischen dem tatsächlichen, zeitlich entwickelten Zustand und der MPF-Näherung:

wobei $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ die Gram-Matrix der Überlappungen zwischen den durch Trotter-Evolution erzeugten Zuständen für verschiedene Schrittzahlen $k_i, k_j$ ist und $L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ die Überlappung mit dem (annähernd) exakten Zustand misst. In diesem Tutorial werden diese Größen mithilfe von Tensor-Netzwerk-Methoden effizient berechnet, insbesondere mithilfe der Backends „ TeNPy-based “ in qiskit_addon_mpf.

Wann sollten MPFs eingesetzt werden?

MPFs sind am vorteilhaftesten, wenn:

• Die Schaltungstiefe ist der Engpass. Wenn die durch die Hardware bedingten Rauschgrenzen die mögliche Tiefe einschränken, verwenden Sie MPFs, um bei flacheren Schaltungen eine höhere effektive Trotter-Genauigkeit zu erzielen.
• Sie benötigen genaue Erwartungswerte, keine vollständige Zustandsvorbereitung. MPFs arbeiten auf der Ebene der Erwartungswerte – sie kombinieren klassische Zahlen, keine Quantenzustände. Sie eignen sich daher ideal für die Schätzung anhand von Beobachtungsdaten bei Verwendung der Estimator-Primitive.
• Du kombinierst eine überschaubare Anzahl von Trotter-Schritten. In der Regel reicht es aus, $r = 3$ – $5$ e unterschiedliche Schrittzahlen $k_j$ zu kombinieren, um mehrere führende Trotter-Fehlerterme zu eliminieren und gleichzeitig $\|x\|_1$ überschaubar zu halten.

Wann MPFs möglicherweise nicht helfen

• Sehr kurze Entwicklungszeiten. Wenn $t$ so klein ist, dass bereits eine einzige Trotter-Formel niedriger Ordnung ausreichend genau ist, entfällt der Mehraufwand für die Ausführung mehrerer Durchläufe.
• Aufgaben zur Vorbereitung auf die Prüfung. MPFs erzeugen einen korrigierten Erwartungswert, keinen korrigierten Quantenzustand. Wenn Sie den aktuellen, zeitlich entwickelten Zustand benötigen (beispielsweise als Eingabe für eine andere Quanten-Subroutine), kommen MPFs nicht in Frage.
• Trotter-Schrittzahlen, die gegen das Konvergenzregime verstoßen. Bei der Herleitung des statischen Koeffizienten wird jeder einzelne Term $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ als Reihe in $t/k_j$ entwickelt; diese Entwicklung konvergiert nur dann gut, wenn $t/k_{\min} \lesssim 1$ gilt. Wird $k_{\min}$ für die gegebenen $t$ zu klein gewählt, liegt der flachste Schaltkreis weit außerhalb des perturbativen Bereichs, die Fehlerterme höherer Ordnung, die durch die MPF nicht aufgehoben werden, werden groß, und die Aufhebung kann große Koeffizienten erfordern. Die „ $L_1$ “-Norm $\|x\|_1$ ist ein praktisches Diagnosekriterium: Wenn $\|x\|_1 \gg 1$ gilt, könnte der Overhead durch die Abtastung $\propto \|x\|_1^2$ die Verringerung des Trotter-Fehlers überwiegen. Weitere Informationen finden Sie im Leitfaden zur Auswahl von Trotter-Stufen.

Inhalt dieses Tutorials

Dieses Tutorial führt Sie in zwei Schritten durch einen MPF-Workflow von Anfang bis Ende. Zunächst zeigt ein Beispiel mit einem Simulator im kleinen Maßstab (10-Qubit-Heisenberg-Kette), wie man das Problem aufstellt, statische und dynamische MPF-Koeffizienten berechnet und die resultierenden Erwartungswerte mit der exakten Diagonalisierung vergleicht. Anschließend zeigt ein groß angelegtes Hardware-Beispiel (50-Qubit-XXZ-Kette), wie man den Code transpiliert, auf einer „ IBM Quantum “-Hardware mit Fehlerminderung ausführt und die Ergebnisse mithilfe der MPF-Koeffizienten nachbearbeitet. Wir verwenden das qiskit_addon_mpf Paket durchgehend zusammen mit den Standard-Tools von Qiskit.

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.22 oder später (pip install qiskit-ibm-runtime)
• Qiskit Aer-Simulator (pip install qiskit-aer)
• MPF-Qiskit-Add-on mit dem „ TeNPy “-Backend (pip install "qiskit-addon-mpf[tenpy]")
• Qiskit-Add-on-Dienstprogramme (pip install qiskit-addon-utils)
• SciPy (pip install scipy)

Einrichtung

Im Folgenden fassen wir alle in diesem Tutorial verwendeten Paketimporte in einer einzigen Zelle zusammen. XXPlusYYGateAußerdem definieren wir einen CollectAndCollapse Transpiler-Durchlauf, der benachbarte rxx und ryy Rotationen zu einer einzigen zusammenfasst. Dieser Durchlauf wird sowohl beim Aufbau der Schaltung in Schritt 1 (um die Anzahl der Gatter gering zu halten) als auch indirekt bei der Extraktion der Schichtstruktur für das dynamische MPF in Schritt 4 angewendet (die Funktion „ TeNPy “ erwartet Zwei-Qubit-Gatter und keine Paare von nicht fusionierten Rotationen).

import warnings

import numpy as np
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)


# --- Helper: collect XX + YY rotations into a single gate ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}


collect_function = partial(
    collect_using_filter_function,
    filter_function=filter_function,
    split_blocks=True,
    min_block_size=1,
)


def collapse_to_xx_plus_yy(block):
    param = 0.0
    for node in block.data:
        param += node.operation.params[0]
    return XXPlusYYGate(param)


collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

pm = PassManager()
pm.append(CollectAndCollapse(collect_function, collapse_function))

Beispiel für einen Simulator im kleinen Maßstab

Schritt 1: Klassische Eingaben einem Quantenproblem zuordnen

Wir beginnen mit einem 10-Qubit-Heisenberg-Modell auf einer Linie und verwenden dabei die Néel-Zustands $\vert 0101\ldots01 \rangle$ als Anfangszustand. Der Hamilton-Operator lautet:

wobei $J$ die Kopplungsstärke zwischen den nächsten Nachbarn ist. Wir messen den ZZ-Korrelator $Z_{L/2-1} Z_{L/2}$ an einem Qubit-Paar in der Mitte der Kette und verwenden Trotter-Schritte $k_j = [1, 2, 4]$ mit einer Produktformel zweiter Ordnung.

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)
print(hamiltonian)

Output:

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

Output:

SparsePauliOp(['IIIIZZIIII'],
              coeffs=[1.+0.j])

# MPF parameters
mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

trotter_times = np.arange(0.5, 1.55, 0.1)
exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)

Trotter-Schaltungen aufbauen

Wir erstellen die Schaltungen, die die approximativen Trotter-Zeitentwicklungen für jeden Zeitpunkt und jede Anzahl von Trotter-Schritten umsetzen. Der CollectAndCollapse im Abschnitt „Setup“ definierte Durchlauf fasst XX- und YY-Rotationen zu einzelnen XX+YY-Gates zusammen, um später eine effizientere Tensor-Netzwerk-Simulation zu ermöglichen.

# Initial Neel state preparation
initial_state_circ = QuantumCircuit(L)
initial_state_circ.x([i for i in range(L) if i % 2 != 0])


all_circs = []
for total_time in trotter_times:
    mpf_trotter_circs = [
        generate_time_evolution_circuit(
            hamiltonian,
            time=total_time,
            synthesis=SuzukiTrotter(reps=num_steps, order=order),
        )
        for num_steps in mpf_trotter_steps
    ]

    mpf_trotter_circs = pm.run(
        mpf_trotter_circs
    )  # Collect XX and YY into XX + YY

    mpf_circuits = [
        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs
    ]
    all_circs.append(mpf_circuits)

mpf_circuits[-1].draw("mpl", fold=-1)

Output:

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

Für das Beispiel im kleinen Maßstab nehmen wir den Aer-Simulator als Ziel. Bevor die Schaltkreise zur Ausführung bereit sind, finden zwei Umwandlungen statt:

1. Gate-Erfassung auf der Ebene der Hamilton-Simulation. XXPlusYYGateIn der Setup-Zelle haben wir einen CollectAndCollapse Durchgang erstellt, ryy der benachbarte rxx Rotationen zu einer einzigen zusammenführt. Diesen Durchlauf haben wir bereits bei der Erstellung der Trotter-Schaltungen in Schritt 1 (dem pm.run(...) Aufruf) angewendet. Dadurch wird sowohl die Anzahl der Zwei-Qubit-Gatter reduziert als auch eine Struktur geschaffen, die sich besser für die Tensor-Netzwerk-Simulation zur späteren Berechnung der dynamischen Koeffizienten eignet.

2. Anpassung an die ISA des Simulators. Im Folgenden führen wir den Qiskit-Preset-Pass-Manager aus, optimization_level=3 um jede Trotter-Schaltung an die Befehlssatzarchitektur (ISA) des Simulators anzupassen.

aer_sim = AerSimulator()
pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)

isa_circs_all_times = [
    pm_sim.run([deepcopy(c) for c in mpf_circuits])
    for mpf_circuits in all_circs
]

Schritt 3: Mit „ Qiskit primitives “ ausführen

Für das Beispiel in kleinem Maßstab lassen wir die ISA-optimierten Trotter-Schaltkreise durch die EstimatorV2 von Aer unterstützte Primitive laufen. Auf diese Weise erhalten wir für jedes $(k_j, t)$ -Paar einen rauschfreien Referenzwert – dies sind die $\langle A \rangle_{k_j}(t)$ -Werte, die der MPF in Schritt 4 kombinieren wird. Wir geben einen Überblick über die Entwicklungszeiträume, damit wir später die vollständige Zeitreihenkurve jeder einzelnen Produktformel und der MPF darstellen können.

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

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

In Schritt 4 wird das MPF tatsächlich erstellt. Auch wenn die Koeffizienten $x_j$ hier berechnet werden (und diese Berechnung bei der dynamischen Variante sehr rechenintensiv sein kann), stellen sie konzeptionell ein klassisches Verfahren dar, um die Quantenmessungen aus Schritt 3 zu einem einzigen korrigierten Erwartungswert zu kombinieren – daher betrachten wir den gesamten Arbeitsablauf rund um die Koeffizienten und die Kombination als Nachbearbeitung.

Um zu beurteilen, wie gut das MPF-Modell die tatsächliche Dynamik abbildet, berechnen wir zunächst die exakten, zeitabhängigen Erwartungswerte, indem wir den Hamilton-Operator direkt potenzieren. Dies ist nur deshalb machbar, weil $L = 10$ gilt; im folgenden Beispiel für groß angelegte Hardware müssen wir uns stattdessen auf Schätzungen mittels Tensor-Netzwerken stützen.

exact_expvals = []
for t in exact_evolution_times:
    exp_H = expm(-1j * t * hamiltonian.to_matrix())
    initial_state = Statevector(initial_state_circ).data
    time_evolved_state = exp_H @ initial_state

    exact_obs = (
        time_evolved_state.conj()
        @ observable.to_matrix()
        @ time_evolved_state
    ).real
    exact_expvals.append(exact_obs)

Statische MPF-Koeffizienten

Statische MPFs verwenden Koeffizienten $x_j$, die unabhängig von der Evolutionszeit, dem Hamilton-Operator und dem Anfangszustand sind. Wir stellen das im Hintergrund beschriebene lineare System $Ax = b$ auf und bestimmen die Koeffizienten. Die Matrix $A$ wird durch die Trotter-Schrittzahlen $k_j$, die Ordnung $\chi$ der Produktformel und die Frage bestimmt, ob die Formel symmetrisch ist (was die Exponenten $\eta_n$ bestimmt).

Für unser einfaches Beispiel verwenden wir $k_j = [1, 2, 4]$ mit einer asymmetrischen Suzuki-Trotter-Formel der Ordnung $2\chi=2$ (also $\chi=1$ und $\eta_n = 2 + n$, was $\eta_0 = 2,\, \eta_1 = 3$ ergibt). Das System lautet dann:

Die erste Zeile gewährleistet die Unvoreingenommenheit ( $\sum_j x_j = 1$ ); die zweite und dritte Zeile heben jeweils die führenden $1/k^2$ und die Trotter-Fehlerterme nächster Ordnung $1/k^3$ auf.

Richten Sie die LSE ein

Wir verwenden setup_static_lse von qiskit_addon_mpf.static bis, um die Matrix $A$ und den Vektor auf der rechten Seite $b$, wie oben beschrieben, zu bilden. Die Matrix $A$ hängt nicht nur von $k_j$ ab, sondern auch von der Wahl der Produktformel – insbesondere von deren Ordnung $\chi$ und davon, ob sie symmetrisch ist. Das symmetric Flag steuert das Exponentenmuster $\eta_n$ (symmetrische Formeln erzeugen ausschließlich Trotter-Fehlerterme mit gerader Potenz; siehe Ref. [1] ). Beachten Sie, dass, wie in [2] gezeigt, die Festlegung symmetric=True von nicht unbedingt erforderlich ist, selbst wenn die zugrunde liegende PF symmetrisch ist – die nicht-symmetrische LSE bleibt gültig (sie erzwingt zusätzliche, unnötige Einschränkungen).

Für unser Beispiel haben wir und symmetric = False bereits in Schritt 1 festgelegt order = 2 .

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)

Überprüfen Sie die erstellte Matrix $A$ und den Vektor $b$, um sicherzustellen, dass sie mit dem oben angegebenen System übereinstimmen.

lse.A

Output:

array([[1.      , 1.      , 1.      ],
       [1.      , 0.25    , 0.0625  ],
       [1.      , 0.125   , 0.015625]])

lse.b

Output:

array([1., 0., 0.])

Mit der LSE in der Hand bestimmen wir die statischen Koeffizienten $x_j$ über lse.solve() (dies ist die direkte Lösung nach $x = A^{-1}b$ ).

mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)

Output:

The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]

Optimieren Sie für „ $x$ “ mithilfe eines exakten Modells

Alternativ zur Berechnung von $x = A^{-1}b$ können Sie setup\_exact\_model verwenden, um eine cvxpy.Problem-Instanz zu erstellen, die die LSE als Nebenbedingungen verwendet und deren optimale Lösung $x$ ergibt.

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)

Output:

[ 0.04761905 -0.57142857  1.52380952]

print(
    "L1 norm of the exact coefficients:",
    np.linalg.norm(coeffs_exact.value, ord=1),
)

Output:

L1 norm of the exact coefficients: 2.1428571428556378

Optimieren Sie für „ $x$ “ mithilfe eines Näherungsmodells

Es kann vorkommen, dass die Norm „ $L_1$ “ für die gewählte Menge von „ $k_j$ “-Werten als zu hoch angesehen wird. Wenn dies der Fall ist und Sie keine anderen Werte für „ $k_j$ “ auswählen können, können Sie eine Näherungslösung verwenden, die die „ $L_1$ “-Norm auf einen gewählten Schwellenwert begrenzt und gleichzeitig „ $\|Ax - b\|$ “ minimiert. Lesen Sie dazu die Anleitung zur Verwendung des Näherungsmodells.

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=1.5
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)

Output:

[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

Dynamische MPF-Koeffizienten

Die statische MPF eliminiert Trotter-Fehlerterme auf eine vom Hamilton-Operator und vom Zustand unabhängige Weise, sodass sie für einen gegebenen Hamilton-Operator und einen gegebenen Anfangszustand nicht zwangsläufig den kleinstmöglichen Approximationsfehler erzeugt. Das dynamische MPF-Verfahren (Lit. [2], [3] ) ermittelt stattdessen zeitabhängige Koeffizienten $x_i(t)$, die den Abstand in der Frobenius-Norm $\|\rho(t) - \mu^D(t)\|_F^2$ zu jedem Zeitpunkt $t$ minimieren. Wie im Abschnitt „Hintergrund“ dargelegt, erfordert dies die Überlappungsmatrix $M_{ij}(t)$ zwischen den nach Trotter entwickelten Zuständen und die Überlappung $L_i(t)$ mit dem exakten Zustand – beides schätzen wir mithilfe von Tensor-Network-Backends ( TeNPy ) in qiskit_addon_mpf.

Um die dynamische LSE einzurichten, benötigen wir drei Komponenten:

1. Eine „Approximate Evolver Factory“, die das Add-on für jeden „ $k_j$ “ ausführt, um „ $\rho_{k_j}(t)$ “ als MPS/MPO zu erzeugen. Wir konstruieren ihn ausgehend von der Schichtstruktur der Ordnung $2$ Trotter-Schaltung (eine Schicht pro slice_by_depth), eingekapselt als mit den LayerwiseEvolver Trunkierungsparametern TeNPy.
2. Eine exakte Evolver-Fabrik, die ein hochgenaues Referenz $\rho(t)$ erzeugt. Wir verwenden eine Suzuki-Trotter-Schaltung vierter Ordnung mit kleinem Zeitschritt (dt=0.1, order=4) als Näherung für die exakte Evolution.
3. Eine Identitätsfabrik und ein MPS im Anfangszustand, die als Startdaten für die Simulation „ TeNPy “ dienen.

Die folgende Zelle erstellt die Fabrik für die approximativen Evolver.

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

dt=0.1Wir approximieren den exakten zeitabhängigen Zustand mit einer Suzuki-Trotter-Formel vierter Ordnung unter Verwendung eines kleinen Zeitschritts. Die Trunkierungsparameter von „ TeNPy “ können die Genauigkeit beeinflussen, daher ist es wichtig, verschiedene Werte zu testen.

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

Abschließend definieren wir ein identity_factory , das den anfänglichen MPO-Zustand liefert, und bereiten den Néel-Ausgangszustand als MPS vor, der mit dem vom geschichteten Trotter-Modell verwendeten Gitter übereinstimmt.

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)


mps_initial_state = MPS_neel_state(models[0].lat)

Nachdem die Fabriken eingerichtet sind, berechnen wir nun die dynamischen Koeffizienten zu jedem Evolutionszeitpunkt. Für jeden $t$ erstellt das Programm die entsprechenden Überlappungsmatrizen über TeNPy, und setup_frobenius_problem gibt einen cvxpy.Problem Wert zurück, setup_dynamic_lse der die Kosten der Frobenius-Norm minimiert. Der Solver gibt auf diesen Zeitpunkt zugeschnittene Koeffizienten $x_j(t)$ zurück; wir sammeln sie in mpf_dynamic_coeffs_list. Wenn der Solver bei einem bestimmten „ $t$ “ fehlschlägt, greifen wir auf Nullkoeffizienten zurück, sodass die Schleife fortgesetzt wird.

mpf_dynamic_coeffs_list = []
for t in trotter_times:
    print(f"Computing dynamic coefficients for time={t}")
    lse = setup_dynamic_lse(
        mpf_trotter_steps,
        t,
        identity_factory,
        exact_factory,
        approx_factory,
        mps_initial_state,
    )
    problem, coeffs = setup_frobenius_problem(lse)
    try:
        problem.solve()
        mpf_dynamic_coeffs_list.append(coeffs.value)
    except Exception as error:
        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))
        print(error, "Calculation Failed for time", t)
    print("")

Output:

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

Kombinieren Sie die Erwartungswerte nach Trotter mit den MPF-Koeffizienten

Nun berechnen wir für jeden Satz von Koeffizienten (static-exact, static-approximate und dynamic) den Wert von „ $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ “, übertragen die Standardfehler pro Schaltung und stellen die resultierenden Zeitreihen gegen die Kurve der exakten Diagonalisierung grafisch dar.

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    trotter_curve, trotter_curve_error = [], []
    for trotter_expvals, trotter_stds in zip(
        mpf_expvals_all_times, mpf_stds_all_times
    ):
        trotter_curve.append(trotter_expvals[k])
        trotter_curve_error.append(trotter_stds[k])

    plt.errorbar(
        trotter_times,
        trotter_curve,
        yerr=trotter_curve_error,
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

# Get expectation values at all times for the static MPF with exact coeffs
exact_mpf_curve, exact_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_exact.value, trotter_stds)
            ]
        )
    )
    exact_mpf_curve_error.append(mpf_std)
    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)

plt.errorbar(
    trotter_times,
    exact_mpf_curve,
    yerr=exact_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Exact",
    color="purple",
)


# Get expectation values at all times for the static MPF with approximate coeffs
approx_mpf_curve, approx_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_approx.value, trotter_stds)
            ]
        )
    )
    approx_mpf_curve_error.append(mpf_std)
    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)

plt.errorbar(
    trotter_times,
    approx_mpf_curve,
    yerr=approx_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Approx",
    color="orange",
)


# Get expectation values at all times for the dynamic MPF
dynamic_mpf_curve, dynamic_mpf_curve_error = [], []
for trotter_expvals, trotter_stds, dynamic_coeffs in zip(
    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(dynamic_coeffs, trotter_stds)
            ]
        )
    )
    dynamic_mpf_curve_error.append(mpf_std)
    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)

plt.errorbar(
    trotter_times,
    dynamic_mpf_curve,
    yerr=dynamic_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Dynamic MPF",
    color="pink",
)


# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

plt.title(f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ vs time")
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

Output:

Die obige Grafik veranschaulicht das Zusammenspiel zwischen Trotter-Fehler und Stichprobenfehler.

• Trotter-Fehler. Die einzelnen Produktformeln (graue Markierungen) weichen mit zunehmender Zeit immer stärker von der exakten Kurve ab. Der Regelkreis „ $k=1$ “ weist die größte Abweichung auf und ist am flachsten, befindet sich jedoch bereits in dem Bereich, in dem $t/k \gtrsim 1$ gilt, sodass der führende Fehlerterm $1/k^{2}$ groß ist. Die MPF-Kombinationen (farbige Markierungen) heben mehrere dieser führenden Trotter-Fehlerterme auf, sodass sie der exakten Kurve viel genauer folgen als jede einzelne „ $k_j$ “-Schaltung. Die verbleibende Lücke spiegelt die Trotter-Terme höherer Ordnung wider, die durch die MPF nicht eliminiert werden: Ein statisches MPF der Ordnung $2$, $r=3$ eliminiert lediglich die ersten beiden Fehlerordnungen, und bei großen $t/k_{\min}$ dominiert schließlich der nicht eliminierte Rest – daher garantiert die MPF nicht, dass sehr flache Schaltkreise zu beliebigen Zeitpunkten genau bleiben.

• Stichprobenfehler. Die breiteren Fehlerbalken bei den MPF-Kurven sind eine direkte Folge der linearen Kombination: Die Übertragung unabhängiger Standardfehler pro Schaltung $\sigma_{k_j}$ ergibt eine Gesamtvarianz $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Je größer also der „ $\|x\|_2$ “ (und in der Praxis der „ $\|x\|_1$ “, den wir steuern) ist, desto mehr Messungen sind erforderlich, um eine bestimmte Zielunsicherheit zu erreichen. Dies ist der Kompromiss, der hinter der Option „approximate-solver“ im Hintergrund steht: Wir begrenzen „ $\|x\|_1$ “, um diesen Mehraufwand überschaubar zu halten. Entscheidend ist, dass der Stichprobenfehler – anders als der Trotter-Fehler – mit $1/\sqrt{N_{\text{shots}}}$ abnimmt, sodass er durch den Einsatz weiterer Schüsse stets verringert werden kann.

Im folgenden Beispiel für eine groß angelegte Hardware-Implementierung tritt Hardware-Rauschen als zusätzliche Fehlerquelle bei jedem „ $\langle A \rangle_{k_j}$ “ auf, das ebenfalls durch die MPF-Koeffizienten verstärkt wird. In diesem Abschnitt werden wir sehen, wie die Fehlerminderung mit MPFs zusammenwirkt.

Beispiel für groß angelegte Hardware

In diesem Abschnitt erweitern wir das Problem über das hinaus, was sich exakt simulieren lässt. Wir reproduzieren einige der in Ref. [3] dargestellten Ergebnisse unter Verwendung einer 50-Qubit-XXZ-Kette zum Zeit $t = 3$. Wir folgen demselben vierstufigen Arbeitsablauf wie im Beispiel im kleinen Maßstab, richten uns nun jedoch auf echte Quantenhardware mit Fehlerminderung aus. Wie in der Vorlage ist jeder Schritt direkt im Code gekennzeichnet, und ein einzelner Schritt kann sich über mehrere Zellen erstrecken, wenn seine Zwischenergebnisse eine genauere Betrachtung wert sind.

Die Abbildung entspricht dem Beispiel im kleinen Maßstab: Definieren Sie einen Hamilton-Operator, wählen Sie Trotter-Parameter aus, berechnen Sie die MPF-Koeffizienten (statisch und dynamisch) und bauen Sie Schaltungen auf. Die wichtigsten Unterschiede sind:

• Ein XXZ-Hamilton-Operator an 50 Knoten mit zufälligen Kopplungen, entnommen aus $\mathcal{U}(0.5, 1.5)$ (Ref. [3] ).
• Eine symmetrische Trotter-Formel zweiter Ordnung mit $k_j = [3, 4, 6]$ (also $\chi=1$, symmetric=True).
• Ein einziger fester Evolutionszeit $t = 3$. Mit $k_{\min}=3$ ergibt sich daraus $t/k_{\min}=1$, wodurch die flachen Anteile im Trotter-Konvergenzbereich verbleiben, in dem das vom MPF verwendete Modell der führenden Fehlerquelle gültig ist.
• Ein zusätzlicher Einzelkreisvergleichslauf mit $k = 10$ Trotter-Schritten, der als Basislinie verwendet wird. Wir haben uns für „ $k = 10$ “ entschieden, da dessen 2-Qubit-Tiefe auf der Hardware größer ist als die des tiefsten MPF-Bestandteils ( $k_{\max}=6$ ) zuzüglich des Overheads für die Ausführung mehrerer MPF-Schaltungen – tief genug, um rauschbegrenzt zu sein. In diesem Bereich wird erwartet, dass die MPF-Kombination die Basislinie mit einer einzelnen Schaltung übertrifft. Es handelt sich um einen Vergleich mit einem „einzelnen tiefen Schaltkreis“ im Vergleich zur MPF-Kombination, nicht um einen Schaltkreis, der auf den effektiven Trotter-Fehler des MPF abzielt (was wesentlich mehr Schritte erfordern würde).

Beachten Sie, dass wir hier zwar noch in Schritt 1 sind (Abbildung und Schaltungsaufbau), in dieser Zelle jedoch neben den statischen auch die dynamischen Koeffizienten vorab berechnen. Die dynamischen Koeffizienten hängen von $H$ und $t$ ab, jedoch nicht von den Quantenmessungen; daher können sie jederzeit vor Schritt 4 berechnet werden. Wir machen das jetzt, um alle MPF-spezifischen Einstellungen an einem Ort zu bündeln.

# -------------------------Step 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# XXZ Hamiltonian with random couplings (Ref. [3])
np.random.seed(0)
even_edges = list(coupling_map.get_edges())[::2]
odd_edges = list(coupling_map.get_edges())[1::2]

Js = np.random.uniform(0.5, 1.5, size=L)
hamiltonian = SparsePauliOp(Pauli("I" * L))
for i, edge in enumerate(even_edges + odd_edges):
    hamiltonian += SparsePauliOp.from_sparse_list(
        [
            ("XX", (edge), 2 * Js[i]),
            ("YY", (edge), 2 * Js[i]),
            ("ZZ", (edge), 4 * Js[i]),
        ],
        num_qubits=L,
    )

observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)

total_time = 3
mpf_trotter_steps = [3, 4, 6]
order = 2
symmetric = True

# Static coefficients
lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(f"Static coefficients: {mpf_coeffs}")
print(f"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}")

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(f"Approximate coefficients: {coeffs_approx.value}")
print(f"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}")

# -------------------------Dynamic coefficients-------------------------
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)

layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 4,
    },
)

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 3,
    },
)


def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)


mps_initial_state = MPS_neel_state(models[0].lat)

print(f"Computing dynamic coefficients for time={total_time}")
lse_dyn = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs_dyn.value
except Exception as error:
    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))
    print(error, "Calculation Failed")

# -------------------------Step 1 (cont): Build circuits-------------------------
mpf_circuits = []
for k in mpf_trotter_steps:
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2])
    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(reps=k, order=order),
        time=total_time,
    )
    circuit.compose(trotter_circ, qubits=range(L), inplace=True)
    mpf_circuits.append(circuit)

# Baseline "single deep circuit" comparison run with k=10 Trotter steps.
# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus
# the overhead of running multiple circuits, pushing it into the noise-limited
# regime where MPF is expected to outperform. It does NOT target the MPF's effective
# Trotter error (which would require many more steps).
comp_circuit = QuantumCircuit(L)
comp_circuit.x([i for i in range(L) if i % 2])
trotter_circ = generate_time_evolution_circuit(
    hamiltonian,
    synthesis=SuzukiTrotter(reps=10, order=order),
    time=total_time,
)
comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)
mpf_circuits.append(comp_circuit)

Output:

Static coefficients: [ 0.42857143 -1.82857143  2.4       ]
L1 norm: 4.65714285714286
Approximate coefficients: [-0.4942491   0.40206845  1.09218065]
L1 norm (approx): 1.9884981979026675
Computing dynamic coefficients for time=3

Nun optimieren wir die Schaltungen für das ausgewählte Backend. Wir verwenden den Qiskit-Preset-Pass-Manager unter optimization_level=3, der automatisch eine geeignete Auswahl physikalischer Qubits auswählt und jede Schaltung auf die Gerätetopologie abbildet.

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

Output:

<IBMBackend('ibm_fez')>

Um komplexere Schaltkreise auf echter Hardware auszuführen, sind umfassende Maßnahmen zur Fehlerbehebung erforderlich. Wir ermöglichen dynamische Entkopplung, Gate- und Mess-Twirling, Messfehlerminderung sowie rauschfreie Extrapolation (ZNE). Beachten Sie, dass die hier verwendeten ZNE-Rauschfaktoren (1, 1.2, 1.4) kleiner sind als in einem Szenario mit flachen Schaltkreisen, da die tiefer liegenden MPF-Bestandteile bereits nahe an der Rauschschwelle liegen und starke Rauschverstärkungen sie über den Punkt hinausdrücken würden, ab dem eine ZNE-Extrapolation noch zuverlässig ist.

Wir reichen alle vier Schaltungen (drei MPF-Komponenten unter $k_j = [3, 4, 6]$ sowie die Basislinie „ $k = 10$ “) in einem einzigen Estimator-Auftrag ein.

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Error suppression/mitigation
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)
estimator.options.resilience.zne.extrapolator = "linear"

estimator.options.environment.job_tags = ["TUT_MPF"]

job_50 = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

Wir entnehmen dem Job-Ergebnis die Erwartungswerte und Standardabweichungen pro Schaltkreis und kombinieren diese dann mit den jeweiligen MPF-Koeffizienten genau wie im Beispiel für den kleinen Maßstab: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$, mit fortgeführter Varianz $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$.

# -------------------------Step 4-------------------------
result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

Output:

[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]
[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]

exact_mpf_std = np.sqrt(
    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])
)
print(
    "Exact static MPF expectation value: ",
    evs[:3] @ mpf_coeffs,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, std[:3])
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs[:3] @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs[:3] @ mpf_dynamic_coeffs,
    "+-",
    dynamic_mpf_std,
)

Output:

Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915
Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166
Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169

sym = {3: "^", 4: "s", 6: "p"}
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs[k],
        yerr=std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

plt.errorbar(
    3,
    evs[-1],
    yerr=std[-1],
    alpha=0.5,
    markersize=8,
    marker="x",
    color="blue",
    label="10 Trotter steps",
)

plt.errorbar(
    4,
    evs[:3] @ mpf_coeffs,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    5,
    evs[:3] @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    6,
    evs[:3] @ mpf_dynamic_coeffs,
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation
plt.axhline(
    y=exact_obs, linestyle="--", color="red", label="Exact time-evolution"
)

plt.title(
    f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ at time {total_time} for the different methods"
)
plt.xlabel("Method")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

Output:

Einige Anmerkungen zu den oben aufgeführten Hardware-Ergebnissen:

• Auf der Hardware ist eine tiefergehende Analyse nicht kostenlos. Die Basislinien der Einzelschaltungen sprechen eine deutliche Sprache: Die Schaltung „ $k = 6$ “ ist im Wesentlichen exakt ( $-0.256$ im Vergleich zur Referenz $-0.244$ ), während die Basislinie der tieferen Schaltung „ $k = 10$* “ schlechter* ist ( $-0.061$, Abweichung von $\sim 0.18$ ) und nicht besser. Sobald der Trotter-Fehler bereits gering ist, führt das Hinzufügen weiterer Schritte meist zu einer Vertiefung der Schaltung und zu einer Zunahme von Gate-Rauschen und Dekohärenz. Genau für diesen Anwendungsfall sind MPFs konzipiert: die Genauigkeit einer tiefen Schaltung allein mit flachen Bausteinen zu erreichen.

• Ein MPF mit kleiner Norm schlägt den Deep Single Circuit. Der annähernd statische MPF (begrenzt auf $\|x\|_1 \approx 2$ ) liegt bei $-0.259$, nur $\sim 0.015$ vom Referenzwert entfernt und damit weit näher als die $k = 10$ -Basislinie. Auch der dynamische MPF ( $-0.127$ ) übertrifft diese Basislinie deutlich. Beide kombinieren zwar nur die flachen „ $k_j = [3, 4, 6]$ “-Schaltkreise, liefern jedoch eine Antwort, die der tiefe Einzelschaltkreis nicht liefern konnte.

• Die Koeffizienten-Norm ist wichtiger als die mathematische Optimalität. Das exakt-statische MPF weist eine Norm von $\|x\|_1 = 4.66$ auf und ist der schlechteste Schätzer von allen ( $-0.567$, mit einer Abweichung von mehr als $0.3$ ): Die große Koeffizienten-Norm verstärkt das Rest-Gate-Rauschen, die Dekohärenz und den ZNE-Fehler bei jedem $\langle A \rangle_{k_j}$ um ungefähr denselben Faktor, wodurch die dadurch erzielte Aufhebung des Trotter-Fehlers zunichte gemacht wird. Durch die Begrenzung der Norm (des approximativ-statischen Lösers, $\|x\|_1 \approx 2$ ) wird diese Überlastung beseitigt und die beste Schätzung erzielt – auch wenn sich dessen Koeffizienten nicht mehr exakt mit dem führenden Trotter-Fehler aufheben.

• Auch einzelne Flachwasser-Wettkämpfe können nach wie vor wettbewerbsfähig sein. Die einzige Komponente „ $k = 6$ “ ( $-0.256$ ) ist hier im Grunde genommen exakt – bei diesem Durchlauf liegt sie sogar geringfügig näher als der annähernd statische MPF. Der Haken daran ist, dass man im Voraus nicht weiß*, welcher* einzelne $k$ genau den „Sweet Spot“ zwischen „konvergiert, aber noch nicht rauschbegrenzt“ trifft, und dass die auf den ersten Blick sichere Entscheidung, einfach tiefer zu gehen ( $k = 10$ ), um die Trotter-Konvergenz zu gewährleisten, genau die ist, die fehlschlägt. Das MPF bietet eine prinzipielle Kombination aus flachen Schaltkreisen, bei der die richtige Tiefe nicht erraten werden muss.

Die praktische Schlussfolgerung lautet, dass MPFs auf der Hardware mit einer starken Fehlerminderung für jeden einzelnen „ $\langle A \rangle_{k_j}$ “ kombiniert werden sollten, die Koeffizienten- $L_1$ -Norm moderat gehalten werden sollte (Verwendung des approximativen Lösers oder des dynamischen MPF) und die Trotter-Schritte $k_j$ so gewählt werden sollten, dass $t/k_{\min} \lesssim 1$ — hier liefert $k_{\min} = 3$ unter $t = 3$ das Ergebnis $t/k_{\min} = 1$, wodurch die Komponenten im konvergenten Bereich bleiben, in dem das Modell für den führenden Fehler, auf das sich der statische MPF stützt, gültig ist. Mit diesen Einstellungen erreichen die MPFs mit kleiner Norm hier die Leistung eines konvergierten Einzelkreis-Modells, während dies bei der naiven Basislinie „einfach tiefer gehen“ nicht der Fall ist, wodurch der in Ref. [3] gezeigte Vorteil von Tiefe gegenüber Genauigkeit wiederhergestellt wird. Beachten Sie außerdem, dass einzelne Durchläufe mit Rauschen behaftet sind – bei einer anderen Einreichung desselben Jobs (oder auf einem anderen Backend) kann sich die genaue Reihenfolge verschieben; die robusten Trends zeigen, dass kleine MPFs ( $\|x\|_1$ ) gut abschneiden, der große, exakt-statische MPF ( $\|x\|_1$ ) durch Hardware-Rauschen verstärkt wird und der übertiefe Einzelkreis durch Rauschen begrenzt ist.

Nächste Schritte

Referenzen

[1] Vázquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Gut konditionierte Mehrproduktformeln für eine hardwarefreundliche Hamilton-Simulation. Quantum, 7, 1067 (2023)

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Trotter-Fehlergrenzen und dynamische Mehrproduktformeln für die Hamilton-Simulation. Physical Review Research, 6(3), 033309 (2024)

[3] Robertson, N. F., et al. Durch Tensor-Netzwerke erweiterte dynamische Mehrproduktformeln. arXiv:2407.17405 (2024)