Multi-Produkt-Formeln zur Reduzierung des Trotter-Fehlers

Geschätzte QPU-Nutzung: Vier Minuten auf einem Heron-r2-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine tatsächliche Laufzeit kann abweichen.)

Hintergrund

Dieses Tutorial zeigt, wie du eine Multi-Produkt-Formel (MPF) einsetzt, um einen geringeren Trotter-Fehler beim Observable zu erzielen, als er durch den tiefsten Trotter-Schaltkreis entsteht, den du tatsächlich ausführst. MPFs reduzieren den Trotter-Fehler der Hamiltonschen Dynamik durch eine gewichtete Kombination mehrerer Schaltkreisausführungen. Betrachte die Aufgabe, Erwartungswerte von Observablen für den Quantenzustand $\rho(t)=e^{-i H t} \rho(0) e^{i H t}$ mit dem Hamiltonian $H$ zu berechnen. Man kann Produkt-Formeln (PFs) verwenden, um die Zeitentwicklung $e^{-i H t}$ wie folgt anzunähern:

• Schreibe den Hamiltonian $H$ als $H=\sum_{a=1}^d F_a,$ wobei $F_a$ hermitesche Operatoren sind, sodass jedes zugehörige Unitär effizient auf einem Quantenprozessor implementiert werden kann.
• Nähere die Terme $F_a$ an, die nicht miteinander kommutieren.

Die Produkt-Formel erster Ordnung (Lie-Trotter-Formel) lautet dann:

$S_1(t):=\prod_{a=1}^d e^{-i F_a t},$

mit dem quadratischen Fehlerterm $S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right)$. Es gibt auch Produkt-Formeln höherer Ordnung (Lie-Trotter-Suzuki-Formeln), die schneller konvergieren und rekursiv definiert sind:

$S_2(t):=\prod_{a=1}^d e^{-i F_a t/2}\prod_{a=1}^d e^{-i F_a t/2}$

$S_{2 \chi}(t):= S_{2 \chi -2}(s_{\chi}t)^2 S_{2 \chi -2}((1-4s_{\chi})t)S_{2 \chi -2}(s_{\chi}t)^2,$

wobei $\chi$ die Ordnung der symmetrischen PF und $s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}$ ist. Bei langen Zeitentwicklungen kann man das Zeitintervall $t$ in $k$ Abschnitte – sogenannte Trotter-Schritte – der Länge $t/k$ unterteilen und die Zeitentwicklung in jedem Abschnitt mit einer Produkt-Formel der Ordnung $\chi$, $S_{\chi}$, annähern. Die PF der Ordnung $\chi$ über $k$ Trotter-Schritte lautet damit:

$S_{\chi}^{k}(t) = \left[ S_{\chi} \left( \frac{t}{k} \right)\right]^k = e^{-i H t}+O\left(t \left( \frac{t}{k} \right)^{\chi} \right)$

wobei der Fehlerterm mit der Anzahl der Trotter-Schritte $k$ und der Ordnung $\chi$ der PF abnimmt.

Für eine ganze Zahl $k \geq 1$ und eine Produkt-Formel $S_{\chi}(t)$ erhält man den näherungsweise zeitentwickelten Zustand $\rho_k(t)$ aus $\rho_0$, indem man $k$ Iterationen der Produkt-Formel $S_{\chi}\left(\frac{t}{k}\right)$ anwendet.

$$\rho_k(t)=S_{\chi}\left(\frac{t}{k}\right)^k \rho_0 S_{\chi}\left(\frac{t}{k}\right)^{-k}$$

$\rho_k(t)$ ist eine Näherung für $\rho(t)$ mit dem Trotter-Approximationsfehler ||$\rho_k(t)-\rho(t) ||$. Betrachten wir eine Linearkombination von Trotter-Näherungen von $\rho(t)$:

$$\mu(t) = \sum_{j}^{l} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{verbleibender Trotter-Fehler} \, ,$$

wobei $x_j$ die Gewichtungskoeffizienten sind, $\rho^{k_j}_j$ die Dichtematrix des reinen Zustands ist, der durch Entwicklung des Anfangszustands mit der Produkt-Formel $S^{k_j}_{\chi}$ mit $k_j$ Trotter-Schritten entsteht, und $j \in {1, ..., l}$ die Anzahl der PFs indiziert, die die MPF bilden. Alle Terme in $\mu(t)$ verwenden dieselbe Produkt-Formel $S_{\chi}(t)$ als Basis. Das Ziel ist es, ||$\rho_k(t)-\rho(t) \|$ zu verbessern, indem man ein $\mu(t)$ mit noch kleinerem $\|\mu(t)-\rho(t)\|$ findet.

• $\mu(t)$ muss kein physikalischer Zustand sein, da $x_i$ nicht positiv sein müssen. Ziel ist es, den Fehler im Erwartungswert der Observablen zu minimieren, nicht einen physikalischen Ersatz für $\rho(t)$ zu finden.
• $k_j$ bestimmt sowohl die Schaltkreistiefe als auch das Niveau der Trotter-Näherung. Kleinere Werte von $k_j$ führen zu flacheren Schaltkreisen, die weniger Schaltkreisfehler verursachen, aber eine weniger genaue Näherung des gewünschten Zustands darstellen.

Der entscheidende Punkt ist, dass der verbleibende Trotter-Fehler von $\mu(t)$ kleiner ist als der Trotter-Fehler, den man durch einfache Verwendung des größten $k_j$-Werts erhalten würde.

Dieser Ansatz lässt sich aus zwei Perspektiven betrachten:

1. Bei einem festen Budget an Trotter-Schritten kannst du Ergebnisse mit insgesamt kleinerem Trotter-Fehler erzielen.
2. Wenn eine Zielanzahl von Trotter-Schritten zu groß ist, um sie auszuführen, kannst du die MPF nutzen, um eine Sammlung von Schaltkreisen mit geringerer Tiefe zu finden, die einen ähnlichen Trotter-Fehler ergeben.

Voraussetzungen

Stelle vor Beginn dieses Tutorials sicher, dass Folgendes installiert ist:

• Qiskit SDK v1.0 oder höher, mit Unterstützung für Visualisierung
• Qiskit Runtime v0.22 oder höher (pip install qiskit-ibm-runtime)
• MPF-Qiskit-Addons (pip install qiskit_addon_mpf)
• Qiskit-Addons-Utils (pip install qiskit_addon_utils)
• Quimb-Bibliothek (pip install quimb)
• Qiskit-Quimb-Bibliothek (pip install qiskit-quimb)
• Numpy v0.21 für paketübergreifende Kompatibilität (pip install numpy==0.21)

Teil I. Kleinskaliges Beispiel

Stabilität der MPF untersuchen

Es gibt keine offensichtliche Einschränkung bei der Wahl der Anzahl der Trotter-Schritte $k_j$, die den MPF-Zustand $\mu(t)$ bilden. Diese müssen jedoch sorgfältig gewählt werden, um Instabilitäten in den aus $\mu(t)$ berechneten Erwartungswerten zu vermeiden. Eine gute Faustregel ist, den kleinsten Trotter-Schritt $k_{\text{min}}$ so zu setzen, dass $t/k_{\text{min}} \lt 1$ gilt. Wenn du mehr darüber erfahren möchtest und wissen willst, wie du deine anderen $k_j$-Werte wählst, lies die Anleitung How to choose the Trotter steps for an MPF.

Im folgenden Beispiel untersuchen wir die Stabilität der MPF-Lösung, indem wir den Erwartungswert der Magnetisierung für einen Bereich von Zeiten mit verschiedenen zeitentwickelten Zuständen berechnen. Konkret vergleichen wir die Erwartungswerte der einzelnen näherungsweisen Zeitentwicklungen mit den entsprechenden Trotter-Schritten und der verschiedenen MPF-Modelle (statische und dynamische Koeffizienten) mit den exakten Werten des zeitentwickelten Observablen. Zunächst definieren wir die Parameter für die Trotter-Formeln und die Entwicklungszeiten:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-mpf qiskit-addon-utils qiskit-aer qiskit-ibm-runtime rustworkx scipy

import numpy as np

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)

Für dieses Beispiel verwenden wir den Néel-Zustand als Anfangszustand $\vert \text{Neel} \rangle = \vert 0101...01 \rangle$ und das Heisenberg-Modell auf einer Kette von 10 Gitterplätzen als Hamilton-Operator für die Zeitentwicklung:

$$\hat{\mathcal{H}}_{Heis} = J \sum_{i=1}^{L-1} \left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ Z_i Z_{(i+1)} \right) \, ,$$

wobei $J$ die Kopplungsstärke für nächste-Nachbar-Kanten ist.

from qiskit.transpiler import CouplingMap
from rustworkx.visualization import graphviz_draw
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
import numpy as np

L = 10

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

# Get a qubit operator describing the Heisenberg field model
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)

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

Das Observable, das wir messen werden, ist die Magnetisierung eines Qubit-Paars in der Mitte der Kette.

from qiskit.quantum_info import SparsePauliOp

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

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

Wir definieren einen Transpiler-Pass, der die XX- und YY-Rotationen im Schaltkreis als einzelnes XX+YY-Gate zusammenfasst. Damit können wir die Spinerhaltungseigenschaften von TeNPy während der MPO-Berechnung nutzen und die Berechnung erheblich beschleunigen.

from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)
from functools import partial

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

Anschließend erstellen wir die Schaltkreise, die die näherungsweisen Trotter-Zeitentwicklungen implementieren.

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

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

Als nächstes berechnen wir die zeitentwickelten Erwartungswerte aus den Trotter-Schaltkreisen.

from copy import deepcopy
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

aer_sim = AerSimulator()
estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for t, mpf_circuits in zip(trotter_times, all_circs):
    mpf_expvals = []
    circuits = [deepcopy(circuit) for circuit in mpf_circuits]
    pm_sim = generate_preset_pass_manager(
        backend=aer_sim, optimization_level=3
    )
    isa_circuits = pm_sim.run(circuits)
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals = [res.data.evs for res in result]
    mpf_stds = [res.data.stds for res in result]
    mpf_expvals_all_times.append(mpf_expvals)
    mpf_stds_all_times.append(mpf_stds)

Wir berechnen außerdem die exakten Erwartungswerte zum Vergleich.

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exact_expvals = []
for t in exact_evolution_times:
    # Exact expectation values
    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 sind solche, bei denen die $x_j$-Werte nicht von der Entwicklungszeit $t$ abhängen. Betrachten wir die PF der Ordnung $\chi = 1$ mit $k_j$ Trotter-Schritten, die sich wie folgt schreiben lässt:

$S_1^{k_j}\left( \frac{t}{k_j} \right)=e^{-i H t}+ \sum_{n=1}^{\infty} A_n \frac{t^{n+1}}{k_j^n}$

wobei $A_n$ Matrizen sind, die von den Kommutatoren der $F_a$-Terme in der Zerlegung des Hamiltonians abhängen. Wichtig ist, dass $A_n$ selbst unabhängig von der Zeit und der Anzahl der Trotter-Schritte $k_j$ sind. Daher lassen sich Fehlerterme niedrigerer Ordnung in $\mu(t)$ durch eine sorgfältige Wahl der Gewichte $x_j$ der Linearkombination aufheben. Um den Trotter-Fehler der ersten $l-1$ Terme (die den größten Beitrag liefern, da sie einer kleineren Anzahl von Trotter-Schritten entsprechen) im Ausdruck für $\mu(t)$ zu eliminieren, müssen die Koeffizienten $x_j$ folgende Gleichungen erfüllen:

$\sum_{j=1}^l x_j = 1$ $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{n}} = 0$

mit $n=0, ... l-2$. Die erste Gleichung stellt sicher, dass der konstruierte Zustand $\mu(t)$ keine Verzerrung aufweist, während die zweite Gleichung die Aufhebung der Trotter-Fehler garantiert. Für Produkt-Formeln höherer Ordnung wird die zweite Gleichung zu $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0$, wobei $\eta = \chi + 2n$ für symmetrische PFs und $\eta = \chi + n$ sonst gilt, mit $n=0, ..., l-2$. Der resultierende Fehler (Refs. [1],[2]) beträgt dann:

$\epsilon = \mathcal{O} \left( \frac{t^{l+1}}{k_1^l} \right).$

Die statischen MPF-Koeffizienten für eine gegebene Menge von $k_j$-Werten zu bestimmen, bedeutet, das lineare Gleichungssystem aus den beiden obigen Gleichungen für die Variablen $x_j$ zu lösen: $Ax=b$. Dabei sind $x$ die gesuchten Koeffizienten, $A$ eine Matrix, die von $k_j$ und der Art der verwendeten PF ($S$) abhängt, und $b$ ein Vektor von Randbedingungen. Konkret gilt:

$A_{0,j} = 1$ $A_{i>0,j} = k_{j}^{-(\chi + s(i-1))}$ $b_0 = 1$ $b_{i>0} = 0$

wobei $\chi$ die order ist, $s$ gleich $2$ ist, wenn symmetric gleich True ist, und sonst $1$, $k_{j}$ die trotter_steps sind und $x$ die zu lösenden Variablen sind. Die Indizes $i$ und $j$ beginnen bei $0$. In Matrixform lässt sich dies so darstellen:

$$A = \begin{bmatrix} A_{0,0} & A_{0,1} & A_{0,2} & ... \\ A_{1,0} & A_{1,1} & A_{1,2} & ... \\ A_{2,0} & A_{2,1} & A_{2,2} & ... \\ ... & ... & ... & ... \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & ... \\ k_{0}^{-(\chi + s(1-1))} & k_{1}^{-(\chi + s(1-1))} & k_{2}^{-(\chi + s(1-1))} & ... \\ k_{0}^{-(\chi + s(2-1))} & k_{1}^{-(\chi + s(2-1))} & k_{2}^{-(\chi + s(2-1))} & ... \\ ... & ... & ... & ... \end{bmatrix}$$

und

$$b = \begin{bmatrix} b_{0} \\ b_{1} \\ b_{2} \\ ... \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ ... \end{bmatrix}$$

Weitere Details findest du in der Dokumentation des Linearen Gleichungssystems (LSE).

Die analytische Lösung für $x$ lautet $x = A^{-1}b$; siehe zum Beispiel Refs. [1] oder [2]. Diese exakte Lösung kann jedoch „schlecht konditioniert" sein, was zu sehr großen L1-Normen unserer Koeffizienten $x$ führt und die MPF-Leistung verschlechtern kann. Alternativ kann man eine Näherungslösung berechnen, die die L1-Norm von $x$ minimiert, um das MPF-Verhalten zu optimieren.

Das LGS aufstellen

Nachdem wir unsere $k_j$-Werte gewählt haben, müssen wir zunächst das LGS $Ax=b$ wie oben beschrieben aufstellen. Die Matrix $A$ hängt nicht nur von $k_j$ ab, sondern auch von unserer Wahl der PF, insbesondere ihrer Ordnung. Zusätzlich kann man berücksichtigen, ob die PF symmetrisch ist oder nicht (siehe [1]), indem man symmetric=True/False setzt. Dies ist jedoch nicht erforderlich, wie Ref. [2] zeigt.

from qiskit_addon_mpf.static import setup_static_lse

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

Gehen wir die oben gewählten Werte durch, um die Matrix $A$ und den Vektor $b$ zu konstruieren. Mit $j=0,1,2$ Trotter-Schritten $k_j = [1, 2, 4]$, Ordnung $\chi = 2$ und nicht-symmetrischen Trotter-Schritten ($s=1$) werden die Matrixelemente von $A$ unterhalb der ersten Zeile durch den Ausdruck $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$ bestimmt, konkret:

$A_{0,0} = A_{0,1} = A_{0,2} = 1$ $A_{1,j} = k_{j}^{-1} \rightarrow A_{1,0} = \frac{1}{1^2}, \;, A_{1,1} = \frac{1}{2^2}, \;, A_{1,2} = \frac{1}{4^2}$ $A_{2,j} = k_{j}^{-2} \rightarrow A_{2,0} = \frac{1}{1^3}, \;, A_{2,1} = \frac{1}{2^3}, \;, A_{2,2} = \frac{1}{4^3}$

bzw. in Matrixform:

$$A = \begin{bmatrix} 1 & 1 & 1\\ 1 & \frac{1}{2^2} & \frac{1}{4^2} \\ 1 & \frac{1}{2^3} & \frac{1}{4^3} \\ \end{bmatrix}$$

Das lässt sich durch Inspektion des lse-Objekts überprüfen:

lse.A

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

Der Constraint-Vektor $b$ hat folgende Elemente: $b_{0} = 1$ $b_1 = b_2 = 0$

Also:

$$b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$$

Und entsprechend in lse:

lse.b

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

Das lse-Objekt besitzt Methoden zur Berechnung der statischen Koeffizienten $x_j$, die das Gleichungssystem erfüllen.

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

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

$x$ mit einem exakten Modell optimieren

Alternativ zur Berechnung von $x=A^{-1}b$ kannst du setup_exact_model verwenden, um eine cvxpy.Problem-Instanz zu konstruieren, die das LGS als Nebenbedingungen verwendet und deren optimale Lösung $x$ liefert.

Im nächsten Abschnitt wird deutlich, warum diese Schnittstelle existiert.

from qiskit_addon_mpf.costs import setup_exact_problem

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

[ 0.04761905 -0.57142857  1.52380952]

Als Indikator dafür, ob eine MPF mit diesen Koeffizienten gute Ergebnisse liefert, kann die L1-Norm verwendet werden (siehe auch Ref. [1]).

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

L1 norm of the exact coefficients: 2.1428571428556378

$x$ mit einem Näherungsmodell optimieren

Es kann vorkommen, dass die L1-Norm für den gewählten Satz von $k_j$-Werten als zu hoch eingestuft wird. Wenn das der Fall ist und du keinen anderen Satz von $k_j$-Werten wählen kannst, kannst du statt einer exakten Lösung eine Näherungslösung für das LGS verwenden.

Verwende dazu einfach setup_approximate_model, um eine andere cvxpy.Problem-Instanz zu konstruieren, die die L1-Norm auf einen gewählten Schwellenwert beschränkt und dabei die Differenz zwischen $Ax$ und $b$ minimiert.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

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

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

Beachte, dass du vollständige Freiheit darüber hast, wie du dieses Optimierungsproblem löst – du kannst den Optimierungs-Solver, seine Konvergenzschwellen und weiteres anpassen. Lies dazu die entsprechende Anleitung How to use the approximate model.

Dynamische MPF-Koeffizienten

Im vorherigen Abschnitt haben wir eine statische MPF eingeführt, die die Standard-Trotter-Näherung verbessert. Diese statische Variante minimiert den Approximationsfehler jedoch nicht notwendigerweise. Konkret ist die statische MPF $\mu^S(t)$ nicht die optimale Projektion von $\rho(t)$ auf den Unterraum, der von den Produkt-Formel-Zuständen $\{\rho_{k_i}(t)\}_{i=1}^r$ aufgespannt wird.

Um diesem Problem zu begegnen, betrachten wir eine dynamische MPF (eingeführt in Ref. [2] und experimentell demonstriert in Ref. [3]), die den Approximationsfehler in der Frobenius-Norm tatsächlich minimiert. Formal minimieren wir:

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr],$$

bezüglich einiger Koeffizienten $x_i(t)$ zu jedem Zeitpunkt $t$. Der optimale Projektor in der Frobenius-Norm ist dann $\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t)$, und wir nennen $\mu^D(t)$ die dynamische MPF. Setzt man die obigen Definitionen ein, ergibt sich:

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr] \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t) \right) \left( \rho(t) - \sum_{j=1}^r x_j(t)\,\rho_{k_j}(t) \right) \bigr] \;=\; \\ = 1 \;+\; \sum_{i,j=1}^r M_{i,j}(t)\,x_i(t)\,x_j(t) \;-\; 2 \sum_{i=1}^r L_i^{\mathrm{exact}}(t)\,x_i(t),$$

wobei $M_{i,j}(t)$ die Gram-Matrix ist, definiert durch:

$$M_{i,j}(t) \;=\; \mathrm{Tr}\bigl[\rho_{k_i}(t)\,\rho_{k_j}(t)\bigr] \;=\; \bigl|\langle \psi_{\mathrm{in}} \!\mid S\bigl(t/k_i\bigr)^{-k_i}\,S\bigl(t/k_j\bigr)^{k_j} \!\mid \psi_{\mathrm{in}} \rangle \bigr|^2.$$

und

$$L_i^{\mathrm{exact}}(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$$

die Überlappung zwischen dem exakten Zustand $\rho(t)$ und jeder Produkt-Formel-Näherung $\rho_{k_i}(t)$ darstellt. In der Praxis können diese Überlappungen aufgrund von Rauschen oder eingeschränktem Zugang zu $\rho(t)$ nur näherungsweise gemessen werden.

Hier ist $\lvert\psi_{\mathrm{in}}\rangle$ der Anfangszustand und $S(\cdot)$ die in der Produkt-Formel angewendete Operation. Indem man die Koeffizienten $x_i(t)$ wählt, die diesen Ausdruck minimieren (und dabei näherungsweise Überlappungsdaten berücksichtigt, wenn $\rho(t)$ nicht vollständig bekannt ist), erhält man die „beste" (im Sinne der Frobenius-Norm) dynamische Näherung von $\rho(t)$ im MPF-Unterraum. Die Größen $L_i(t)$ und $M_{i,j}(t)$ lassen sich effizient mit Tensornetzwerk-Methoden berechnen [3]. Das MPF-Qiskit-Addon bietet mehrere „Backends" für diese Berechnung. Das folgende Beispiel zeigt die flexibelste Vorgehensweise; die Dokumentation des TeNPy-lagerbasierten Backends erläutert dies ebenfalls ausführlich. Für diese Methode startest du vom Schaltkreis, der die gewünschte Zeitentwicklung implementiert, und erstellst Modelle, die diese Operationen anhand der Schichten des entsprechenden Schaltkreises darstellen. Schließlich wird ein Evolver-Objekt erzeugt, mit dem die zeitentwickelten Größen $M_{i,j}(t)$ und $L_i(t)$ berechnet werden können. Wir beginnen mit der Erstellung des Evolver-Objekts für die näherungsweise Zeitentwicklung (ApproxEvolverFactory), die durch die Schaltkreise implementiert wird. Achte dabei besonders auf die Variable order, damit die Werte übereinstimmen. Beachte, dass bei der Erzeugung der Schaltkreise für die näherungsweise Zeitentwicklung Platzhalterwerte für time = 1.0 und die Anzahl der Trotter-Schritte (reps=1) verwendet werden. Die korrekten Näherungsschaltkreise werden dann vom dynamischen Problemlöser in setup_dynamic_lse erzeugt.

from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.backends.tenpy_layers import LayerModel
from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver
from functools import partial

# 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,
    },
)

warnung

Die Optionen von LayerwiseEvolver, die die Details der Tensornetzwerk-Simulation bestimmen, müssen sorgfältig gewählt werden, um ein schlecht definiertes Optimierungsproblem zu vermeiden.

Dann richten wir den exakten Evolver ein (zum Beispiel ExactEvolverFactory), der ein Evolver-Objekt zurückgibt, das die wahre oder „Referenz"-Zeitentwicklung berechnet. Realistischerweise würde man die exakte Entwicklung durch eine Suzuki–Trotter-Formel höherer Ordnung oder eine andere zuverlässige Methode mit kleinem Zeitschritt annähern. Im Folgenden nähern wir den exakt zeitentwickelten Zustand mit einer Suzuki-Trotter-Formel vierter Ordnung mit kleinem Zeitschritt dt=0.1 an, was bedeutet, dass die Anzahl der Trotter-Schritte zur Zeit $t$ gleich $k=t/dt$ ist. Wir geben außerdem einige TeNPy-spezifische Trunkierungsoptionen an, um die maximale Bond-Dimension des zugrunde liegenden Tensornetzwerks sowie die minimalen Singulärwerte der aufgeteilten Tensornetzwerk-Bonds zu begrenzen. Diese Parameter können die Genauigkeit des mit den dynamischen MPF-Koeffizienten berechneten Erwartungswerts beeinflussen, daher ist es wichtig, einen Bereich von Werten zu erkunden, um die optimale Balance zwischen Rechenzeit und Genauigkeit zu finden. Beachte, dass die Berechnung der MPF-Koeffizienten nicht auf dem Erwartungswert der PF aus der Hardware-Ausführung basiert und daher in der Nachverarbeitung angepasst werden kann.

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,
    },
)

Als nächstes erstellst du den Anfangszustand des Systems in einem mit TeNPy kompatiblen Format (zum Beispiel MPS_neel_state=$\vert 0101...01 \rangle$). Damit wird der Vielteilchen-Wellenfunktionszustand $\lvert\psi_{\mathrm{in}}\rangle$, den du zeitlich entwickelst, als Tensor aufgesetzt.

from qiskit_addon_mpf.backends.tenpy_tebd import MPOState
from qiskit_addon_mpf.backends.tenpy_tebd import MPS_neel_state

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

mps_initial_state = MPS_neel_state(models[0].lat)

Für jeden Zeitschritt $t$ richten wir das dynamische lineare Gleichungssystem mit der Methode setup_dynamic_lse ein. Das zugehörige Objekt enthält die Informationen über das dynamische MPF-Problem: lse.A liefert die Gram-Matrix $M$, während lse.b die Überlappung $L$ liefert. Das LGS kann dann (sofern es nicht schlecht definiert ist) gelöst werden, um die dynamischen Koeffizienten mit setup_frobenius_problem zu finden. Wichtig ist, den Unterschied zu den statischen Koeffizienten zu beachten, die nur von den Details der verwendeten Produkt-Formel abhängen und unabhängig von den Details der Zeitentwicklung (Hamiltonian und Anfangszustand) sind.

from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import setup_frobenius_problem

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

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

Abschließend werden diese Erwartungswerte über die Entwicklungszeit hinweg geplottet.

import matplotlib.pyplot as plt

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
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",
    color="orange",
    label="Static MPF - Approximate",
)

# # 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",
    color="pink",
    label="Dynamic MPF",
)

plt.plot(
    exact_evolution_times,
    exact_expvals,
    lw=3,
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) as a function of time"
)
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend()
plt.grid()

In Fällen wie dem obigen Beispiel, in denen die PF mit $k=1$ zu allen Zeiten schlecht abschneidet, wird auch die Qualität der dynamischen MPF-Ergebnisse stark beeinträchtigt. In solchen Situationen ist es sinnvoll, die Möglichkeit zu untersuchen, einzelne PFs mit einer höheren Anzahl von Trotter-Schritten zu verwenden, um die Gesamtqualität der Ergebnisse zu verbessern. In diesen Simulationen sehen wir das Zusammenspiel verschiedener Fehlerarten: Fehler durch endliches Sampling und Trotter-Fehler durch die Produkt-Formeln. MPF hilft, den Trotter-Fehler der Produkt-Formeln zu reduzieren, verursacht im Vergleich zu den Produkt-Formeln jedoch einen höheren Sampling-Fehler. Das kann vorteilhaft sein, da Produkt-Formeln den Sampling-Fehler durch erhöhtes Sampling verringern können, der systematische Fehler durch die Trotter-Näherung jedoch bestehen bleibt.

Ein weiteres interessantes Verhalten, das wir im Plot beobachten können, ist, dass der Erwartungswert der PF für $k=1$ zu Zeiten, bei denen $t/k > 1$ gilt, unregelmäßig zu werden beginnt (zusätzlich dazu, keine gute Näherung für den exakten Wert zu sein), wie in der Anleitung zur Wahl der Anzahl der Trotter-Schritte erläutert.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Betrachten wir nun einen einzelnen Zeitpunkt $t=1.0$ und berechnen den Erwartungswert der Magnetisierung mit den verschiedenen Methoden auf einem QPU. Die Wahl von $t$ wurde so getroffen, dass der Unterschied zwischen den Methoden möglichst groß ist und ihre relative Wirksamkeit gut sichtbar wird. Um den Zeitbereich zu bestimmen, für den dynamisches MPF garantiert kleinere Fehler liefert als jede der einzelnen Trotter-Formeln im Multi-Produkt, können wir den „MPF-Test" anwenden – siehe Gleichung (17) und den umgebenden Text in [3].

Trotter-Circuits aufbauen

An diesem Punkt haben wir unsere Entwicklungskoeffizienten $x$ bestimmt. Was jetzt noch fehlt, ist die Erzeugung der Trotterisierten Quantencircuits. Auch hier hilft das Modul qiskit_addon_utils.problem_generators mit einer praktischen Funktion weiter:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

total_time = 1.0
mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial Neel state preparation
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2 != 0])

    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=order, reps=k),
        time=total_time,
    )

    circuit.compose(trotter_circ, inplace=True)

    mpf_circuits.append(pm.run(circuit))

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

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

Kehren wir zur Berechnung des Erwartungswerts für einen einzelnen Zeitpunkt zurück. Wir wählen ein Backend für die Ausführung des Experiments auf echter Hardware aus.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
print(backend)

qubits = list(range(backend.num_qubits))

Anschließend entfernen wir Ausreißer-Qubits aus der Coupling Map, damit die Layout-Phase des Transpilers diese nicht berücksichtigt. Im folgenden Code verwenden wir die im target-Objekt gespeicherten Backend-Eigenschaften und entfernen Qubits, deren Meßfehler oder Zwei-Qubit-Gate-Fehler einen bestimmten Schwellwert (max_meas_err, max_twoq_err) überschreitet oder deren $T_2$-Zeit (die den Kohärenzverlust bestimmt) unter einem bestimmten Schwellwert (min_t2) liegt.

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.012
min_t2 = 40
max_twoq_err = 0.005

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 10

Wir wählen max_meas_err, min_t2 und max_twoq_err so, dass wir eine ausreichend große Teilmenge von Qubits finden, auf der der Circuit ausgeführt werden kann. In unserem Fall reicht es, eine 1D-Kette aus 10 Qubits zu finden.

cust_cmap.draw()

Danach können wir den Circuit und den Observablen auf die physischen Qubits des Geräts abbilden.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

print(transpiled_circuits[-1].depth(lambda x: x.operation.num_qubits == 2))
print(transpiled_circuits[-1].count_ops())
transpiled_circuits[-1].draw("mpl", idle_wires=False, fold=False)

51
OrderedDict([('sx', 310), ('rz', 232), ('cz', 132), ('x', 19)])

Schritt 3: Ausführung mit Qiskit-Primitiven

Mit dem Estimator-Primitiv können wir den Erwartungswert vom QPU schätzen. Wir führen die optimierten AQC-Circuits mit zusätzlichen Fehlerunterdrückungs- und Fehlerminderungstechniken aus.

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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, 3, 5)
estimator.options.resilience.zne.extrapolator = ("exponential", "linear")

estimator.options.environment.job_tags = ["mpf small"]

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

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

Der einzige Nachverarbeitungsschritt besteht darin, die von den Qiskit-Runtime-Primitiven bei verschiedenen Trotter-Schritten erhaltenen Erwartungswerte mit den jeweiligen MPF-Koeffizienten zu kombinieren. Für einen Observablen $A$ gilt:

$\langle A \rangle_{\text{mpf}} = \text{Tr} [A \mu(t)] = \sum_{j} x_j \text{Tr} [A \rho_{k_j}] = \sum_{j} x_j \langle A \rangle_j$

Zuerst extrahieren wir die einzelnen Erwartungswerte, die für jeden der Trotter-Circuits berechnet wurden:

result_exp = job.result()
evs_exp = [res.data.evs for res in result_exp]
evs_std = [res.data.stds for res in result_exp]

print(evs_exp)

[array(-0.06361607), array(-0.23820448), array(-0.50271805)]

Anschließend kombinieren wir diese Werte einfach mit unseren MPF-Koeffizienten, um die Gesamt-Erwartungswerte des MPF zu erhalten. Wir tun dies für jede der verschiedenen Arten, mit denen wir $x$ berechnet haben.

exact_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_exact.value, evs_std)
        ]
    )
)
print(
    "Exact static MPF expectation value: ",
    evs_exp @ coeffs_exact.value,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, evs_std)
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs_exp @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs_list[7], evs_std)
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs_exp @ mpf_dynamic_coeffs_list[7],
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.6329590442738475 +- 0.012798249760406036
Approximate static MPF expectation value:  -0.5690390035339492 +- 0.010459559917168473
Dynamic MPF expectation value:  -0.4655579758795695 +- 0.007639139186720507

Für dieses kleine Problem können wir den exakten Referenzwert schließlich mit scipy.linalg.expm wie folgt berechnen:

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exp_H = expm(-1j * total_time * hamiltonian.to_matrix())

initial_state_circuit = QuantumCircuit(L)
initial_state_circuit.x([i for i in range(L) if i % 2 != 0])
initial_state = Statevector(initial_state_circuit).data

time_evolved_state = exp_H @ initial_state

exact_obs = (
    time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
)
print("Exact expectation value ", exact_obs.real)

Exact expectation value  -0.39909900734489434

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs_exp[k],
        yerr=evs_std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )  # , , )

plt.errorbar(
    3,
    evs_exp @ coeffs_exact.value,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    4,
    evs_exp @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    5,
    evs_exp @ mpf_dynamic_coeffs_list[7],
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

plt.axhline(
    y=exact_obs.real,
    linestyle="--",
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

Im obigen Beispiel schneidet die dynamische MPF-Methode beim Erwartungswert am besten ab und übertrifft das Ergebnis, das wir allein mit der höchsten Anzahl an Trotter-Schritten erhalten hätten. Obwohl die verschiedenen MPF-Techniken nicht immer einen verbesserten Erwartungswert gegenüber der höchsten Anzahl von Trotter-Schritten erzielen (wie das exakte und das approximative Modell im obigen Diagramm zeigen), erfasst die Standardabweichung dieser Werte gut die erhöhte Varianz, die durch die MPF-Technik entsteht. Dies verdeutlicht die Unsicherheit rund um den erhaltenen Erwartungswert, der stets den Erwartungswert einschließt, den wir von einer exakten Zeitentwicklung des Systems erwarten würden. Demgegenüber gelingt es den mit einer geringeren Anzahl von Trotter-Schritten berechneten Erwartungswerten nicht, den exakten Wert innerhalb ihrer Unsicherheit zu erfassen – sie liefern also mit hoher Konfidenz das falsche Ergebnis.

def relative_error(ev, exact_ev):
    return abs(ev - exact_ev)

relative_error_k = [relative_error(ev, exact_obs.real) for ev in evs_exp]
relative_error_mpf = relative_error(evs_exp @ mpf_coeffs, exact_obs.real)
relative_error_approx_mpf = relative_error(
    evs_exp @ coeffs_approx.value, exact_obs.real
)
relative_error_dynamic_mpf = relative_error(
    evs_exp @ mpf_dynamic_coeffs_list[7], exact_obs.real
)

print("relative error for each trotter steps", relative_error_k)
print("relative error with MPF exact coeffs", relative_error_mpf)
print("relative error with MPF approx coeffs", relative_error_approx_mpf)
print("relative error with MPF dynamic coeffs", relative_error_dynamic_mpf)

relative error for each trotter steps [0.33548293650112293, 0.16089452939226306, 0.10361904247828346]
relative error with MPF exact coeffs 0.2338600369291003
relative error with MPF approx coeffs 0.16993999618905486
relative error with MPF dynamic coeffs 0.06645896853467514

Teil II: Größer skalieren

Lass uns das Problem über das hinaus skalieren, was sich exakt simulieren lässt. In diesem Abschnitt konzentrieren wir uns darauf, einige der Ergebnisse aus Ref. [3] nachzuvollziehen.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Hamiltonian

Für das großskalige Beispiel verwenden wir das XXZ-Modell auf einer Kette von 50 Gitterplätzen:

$$\hat{\mathcal{H}}_{XXZ} = \sum_{i=1}^{L-1} J_{i,(i+1)}\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\cdot Z_i Z_{(i+1)} \right) \, ,$$

wobei $J_{i,(i+1)}$ ein zufälliger Koeffizient für die Kante $(i, i+1)$ ist. Dies ist der Hamiltonian, der in der Demonstration aus Ref. [3] verwendet wird.

L = 50
# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

import numpy as np
from qiskit.quantum_info import SparsePauliOp, Pauli

# Generate random coefficients for our XXZ Hamiltonian
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,
    )

print(hamiltonian)

SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'XXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'YYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
              coeffs=[1.        +0.j, 2.09762701+0.j, 2.09762701+0.j, 4.19525402+0.j,
 2.43037873+0.j, 2.43037873+0.j, 4.86075747+0.j, 2.20552675+0.j,
 2.20552675+0.j, 4.4110535 +0.j, 2.08976637+0.j, 2.08976637+0.j,
 4.17953273+0.j, 1.8473096 +0.j, 1.8473096 +0.j, 3.6946192 +0.j,
 2.29178823+0.j, 2.29178823+0.j, 4.58357645+0.j, 1.87517442+0.j,
 1.87517442+0.j, 3.75034885+0.j, 2.783546  +0.j, 2.783546  +0.j,
 5.567092  +0.j, 2.92732552+0.j, 2.92732552+0.j, 5.85465104+0.j,
 1.76688304+0.j, 1.76688304+0.j, 3.53376608+0.j, 2.58345008+0.j,
 2.58345008+0.j, 5.16690015+0.j, 2.05778984+0.j, 2.05778984+0.j,
 4.11557968+0.j, 2.13608912+0.j, 2.13608912+0.j, 4.27217824+0.j,
 2.85119328+0.j, 2.85119328+0.j, 5.70238655+0.j, 1.14207212+0.j,
 1.14207212+0.j, 2.28414423+0.j, 1.1742586 +0.j, 1.1742586 +0.j,
 2.3485172 +0.j, 1.04043679+0.j, 1.04043679+0.j, 2.08087359+0.j,
 2.66523969+0.j, 2.66523969+0.j, 5.33047938+0.j, 2.5563135 +0.j,
 2.5563135 +0.j, 5.112627  +0.j, 2.7400243 +0.j, 2.7400243 +0.j,
 5.48004859+0.j, 2.95723668+0.j, 2.95723668+0.j, 5.91447337+0.j,
 2.59831713+0.j, 2.59831713+0.j, 5.19663426+0.j, 1.92295872+0.j,
 1.92295872+0.j, 3.84591745+0.j, 2.56105835+0.j, 2.56105835+0.j,
 5.12211671+0.j, 1.23654885+0.j, 1.23654885+0.j, 2.4730977 +0.j,
 2.27984204+0.j, 2.27984204+0.j, 4.55968409+0.j, 1.28670657+0.j,
 1.28670657+0.j, 2.57341315+0.j, 2.88933783+0.j, 2.88933783+0.j,
 5.77867567+0.j, 2.04369664+0.j, 2.04369664+0.j, 4.08739329+0.j,
 1.82932388+0.j, 1.82932388+0.j, 3.65864776+0.j, 1.52911122+0.j,
 1.52911122+0.j, 3.05822245+0.j, 2.54846738+0.j, 2.54846738+0.j,
 5.09693476+0.j, 1.91230066+0.j, 1.91230066+0.j, 3.82460133+0.j,
 2.1368679 +0.j, 2.1368679 +0.j, 4.2737358 +0.j, 1.0375796 +0.j,
 1.0375796 +0.j, 2.0751592 +0.j, 2.23527099+0.j, 2.23527099+0.j,
 4.47054199+0.j, 2.22419145+0.j, 2.22419145+0.j, 4.44838289+0.j,
 2.23386799+0.j, 2.23386799+0.j, 4.46773599+0.j, 2.88749616+0.j,
 2.88749616+0.j, 5.77499231+0.j, 2.3636406 +0.j, 2.3636406 +0.j,
 4.7272812 +0.j, 1.7190158 +0.j, 1.7190158 +0.j, 3.4380316 +0.j,
 1.87406391+0.j, 1.87406391+0.j, 3.74812782+0.j, 2.39526239+0.j,
 2.39526239+0.j, 4.79052478+0.j, 1.12045094+0.j, 1.12045094+0.j,
 2.24090189+0.j, 2.33353343+0.j, 2.33353343+0.j, 4.66706686+0.j,
 2.34127574+0.j, 2.34127574+0.j, 4.68255148+0.j, 1.42076512+0.j,
 1.42076512+0.j, 2.84153024+0.j, 1.2578526 +0.j, 1.2578526 +0.j,
 2.51570519+0.j, 1.6308567 +0.j, 1.6308567 +0.j, 3.2617134 +0.j])

Als Observable wählen wir $Z_{24}Z_{25}$, wie im unteren Bereich von Abb. 5 in Ref. [3] zu sehen.

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

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

Trotter-Schritte wählen

Das in Abb. 4 von Ref. [3] gezeigte Experiment verwendet $k_j = [2, 3, 4]$ symmetrische Trotter-Schritte der Ordnung $2$. Wir konzentrieren uns auf die Ergebnisse für den Zeitpunkt $t=3$, bei dem die MPF und eine PF mit einer höheren Anzahl von Trotter-Schritten (hier 6) denselben Trotter-Fehler aufweisen. Der MPF-Erwartungswert wird jedoch aus Circuits berechnet, die den niedrigeren Trotter-Schritten entsprechen und daher flacher sind. In der Praxis erwarten wir, dass der aus den MPF-Circuits berechnete experimentelle Erwartungswert näher am theoretischen liegt — auch wenn MPF und der tiefere Trotter-Circuit denselben Trotter-Fehler haben —, weil die flacheren Circuits weniger anfällig für Hardware-Rauschen sind als der Circuit mit der höheren Trotter-Schritt-PF.

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

Das LGS aufstellen

Hier betrachten wir die statischen MPF-Koeffizienten für dieses Problem.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)
print("L1 norm:", np.linalg.norm(mpf_coeffs, ord=1))

The static coefficients associated with the ansatze are: [ 0.26666667 -2.31428571  3.04761905]
L1 norm: 5.628571428571431

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

[-0.24255546 -0.25744454  1.5       ]
L1 norm of the approximate coefficients: 2.0

Dynamische Koeffizienten

# 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": 4,
    },
)

# Create exact time-evolution circuits
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)
]

# Create the time-evolution object
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)

lse = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs = setup_frobenius_problem(lse)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs.value
except Exception as error:
    print(error, "Calculation Failed for time", total_time)
print("")

Die Trotter-Circuits unserer MPF-Zerlegung aufbauen

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial state preparation |1010..>
    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)

Trotter-Circuit mit vergleichbarem Trotter-Fehler wie die MPF aufbauen

k = 6

# Initial state preparation |1010..>
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=k, order=order),
    time=total_time,
)

comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)

mpf_circuits.append(comp_circuit)

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

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.055
min_t2 = 30
max_twoq_err = 0.01

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 73

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

Schritt 3: Mit Qiskit Primitives ausführen

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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 = ["mpf large"]

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

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

result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

[array(-0.08034071), array(-0.00605026), array(-0.15345759), array(-0.18127293)]
[array(0.04482517), array(0.03438413), array(0.21540776), array(0.21520829)]

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

Exact static MPF expectation value:  -0.47510243192011536 +- 0.6613940032465087
Approximate static MPF expectation value:  -0.20914170384216998 +- 0.32341567460419135
Dynamic MPF expectation value:  -0.07994951978722761 +- 0.07423091963310202

sym = {2: "^", 3: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
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="6 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"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

Bei der Ausführung von Circuits auf echter Hardware können zusätzliche Herausforderungen auftreten, wenn es darum geht, genaue Erwartungswerte zu erhalten — bedingt durch Hardware-Rauschen. Dieses ist im MPF-Formalismus nicht berücksichtigt und kann sich nachteilig auf die MPF-Lösung auswirken. Das könnte zum Beispiel der Grund dafür sein, dass die dynamischen Koeffizienten im Plot keinen besseren Schätzwert liefern als die approximativen statischen Koeffizienten. Das heißt: Der approximative Evolver, der den approximativen Circuit simuliert, spiegelt die Ergebnisse der tatsächlichen Circuit-Ausführung unter Hardware-Rauschen nicht akkurat wider. Aus diesen Gründen empfiehlt es sich, verschiedene Fehlerminderungstechniken zu kombinieren, um Ergebnisse zu erzielen, die so nah wie möglich an den idealen Werten jeder Produktformel liegen. Das macht die Vorteile des MPF-Ansatzes konsistent sichtbar.

Insgesamt liefern die approximativen statischen Koeffizienten nach wie vor eine genauere Lösung als die Produktformel mit einer höheren Anzahl von Trotter-Schritten bei gleichem Trotter-Fehler in der rauschfreien Umgebung.

Es ist außerdem wichtig zu beachten, dass im Beispiel, das das Experiment aus Ref. [3] nachvollzieht, der Zeitpunkt $t=3$ jenseits der Grenze liegt, ab der erwartet wird, dass die PF mit $k=2$ gut funktioniert — nämlich $t/k>1$, wie in diesem Leitfaden erläutert.

Referenzen

[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. (2023). Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 7, 1067.

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. (2024). Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 6(3), 033309.

[3] Robertson, N. F., et al. (2024). Tensor network enhanced dynamic multiproduct formulas. arXiv preprint arXiv:2407.17405.