Erste Schritte mit Multi-Produkt-Formeln (MPF)

Erste Schritte mit Multi-Produkt-Formeln (MPFs)

Paketversionen

Der Code auf dieser Seite wurde mit den folgenden Anforderungen entwickelt. Wir empfehlen, diese oder neuere Versionen zu verwenden.

qiskit[all]~=2.3.0
qiskit-addon-utils~=0.3.0
qiskit-addon-mpf~=0.3.0
scipy~=1.16.3

Dieser Leitfaden zeigt, wie du das Paket qiskit-addon-mpf verwendest, und nutzt dabei die Zeitentwicklung eines Ising-Modells als Beispiel. Mit diesem Paket kannst du eine Multi-Produkt-Formel (MPF) erstellen, die einen geringeren Trotter-Fehler bei Observable-Messungen erreicht. Die bereitgestellten Werkzeuge ermöglichen es dir, die Gewichte einer gewählten MPF zu bestimmen, mit denen die geschätzten Erwartungswerte aus mehreren Zeitentwicklungs-Circuits – jeweils mit einer unterschiedlichen Anzahl von Trotter-Schritten – neu kombiniert werden können.

Betrachte zunächst den Hamiltonian eines Ising-Modells mit 10 Gitterstellen:

$H_{\text{Ising}} = \sum_{i=1}^9 J_{i,(i+1)}Z_iZ_{i+1} + \sum_{i=1}^{10} h_i X_i$

wobei $J_{i,(i+1)}$ die Kopplungsstärke und $h_i$ eine externe Magnetfeldstärke ist. Für die Problemstellung wird die Gesamtmagnetisierung des Systems als zu messende Observable gewählt:

$\langle M \rangle = \sum_{i=1}^{10} \langle Z_i \rangle.$

Der folgende Code-Ausschnitt bereitet den Hamiltonian der Ising-Kette mithilfe des Pakets qiskit-addon-utils vor und definiert die zu messende Observable.

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

from qiskit.transpiler import CouplingMap
from qiskit.synthesis import SuzukiTrotter
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import StatevectorEstimator
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import generate_preset_pass_manager
from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
)
from qiskit_addon_mpf.static import setup_static_lse

from scipy.linalg import expm
import numpy as np

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(10, bidirectional=False)

# Get a qubit operator describing the Ising field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(0.0, 0.0, 1.0),
    ext_magnetic_field=(0.4, 0.0, 0.0),
)
print(f"Hamiltonian:\n {hamiltonian}\n")

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L
)
print(f"Observable:\n {observable}")

Hamiltonian:
 SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],
              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, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,
 0.4+0.j, 0.4+0.j, 0.4+0.j])

Observable:
 SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j,
 0.05+0.j, 0.05+0.j, 0.05+0.j])

Als Nächstes bereitest du die MPF vor. Die erste Entscheidung ist, ob die Koeffizienten statisch (zeitunabhängig) oder dynamisch sein sollen; dieses Tutorial verwendet eine statische MPF. Die nächste Entscheidung ist die Wahl der $k_j$ -Werte. Diese legt fest, wie viele Terme die MPF enthält und wie viele Trotter-Schritte jeder Term zur Simulation der Zeitentwicklung verwendet. Die Wahl der $k_j$ -Werte ist heuristischer Natur, daher musst du eine eigene Menge „guter" $k_j$ -Werte ermitteln. Hinweise zur Findung einer guten Wertemenge findest du am Ende der Einstiegsseite.

Sind die $k_j$ -Werte bestimmt, kann das Gleichungssystem $Ax=b$ aufgestellt werden, das gelöst werden soll. Die Matrix $A$ wird ebenfalls durch die zu verwendende Produktformel bestimmt. Die Wahlmöglichkeiten sind ihre Ordnung – hier auf $2$ gesetzt – sowie ob die Produktformel symmetrisch sein soll – hier auf True gesetzt. Der folgende Code-Ausschnitt wählt eine Gesamtzeit für die Systementwicklung, die zu verwendenden $k_j$ -Werte sowie das zu lösende Gleichungssystem mithilfe der Methode qiskit_addon_mpf.static.setup_static_lse.

time = 8.0
trotter_steps = (8, 12, 19)

lse = setup_static_lse(trotter_steps, order=2, symmetric=True)
print(lse)

LSE(A=array([[1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
       [1.56250000e-02, 6.94444444e-03, 2.77008310e-03],
       [2.44140625e-04, 4.82253086e-05, 7.67336039e-06]]), b=array([1., 0., 0.]))

Sobald das lineare Gleichungssystem instanziiert wurde, kann es entweder exakt oder über ein Näherungsmodell mittels einer Summe von Quadraten (bzw. der Frobenius-Norm für dynamische Koeffizienten; weitere Informationen findest du in der API-Referenz) gelöst werden. Die Wahl eines Näherungsmodells ergibt sich typischerweise, wenn die Norm der Koeffizienten für die gewählte Menge von $k_j$ -Werten als zu hoch eingestuft wird und keine andere Menge von $k_j$ -Werten gewählt werden kann. Dieser Leitfaden demonstriert beide Ansätze zum Vergleich der Ergebnisse.

model_exact, coeffs_exact = setup_exact_problem(lse)
model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=3.0
)
model_exact.solve()
model_approx.solve()
print(f"Exact solution: {coeffs_exact.value}")
print(f"Approximate solution: {coeffs_approx.value}")

Exact solution: [ 0.17239057 -1.19447005  2.02207947]
Approximate solution: [-0.40454257  0.57553173  0.8290123 ]

hinweis

Das LSE-Objekt besitzt auch eine Methode LSE.solve(), die das Gleichungssystem exakt löst. Der Grund, warum in diesem Leitfaden setup_exact_problem() verwendet wird, ist die Demonstration der Schnittstelle, die von den anderen Näherungsmethoden bereitgestellt wird.

Trotter-Circuits einrichten und ausführen

Nachdem die Koeffizienten $x_j$ ermittelt wurden, besteht der letzte Schritt darin, die Zeitentwicklungs-Circuits für die Ordnung und die gewählte Schrittzahl $k_j$ der MPF zu generieren. Das Paket qiskit-addon-utils kann diesen Prozess beschleunigen.

circuits = []
for k in trotter_steps:
    circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=2, reps=k),
        time=time,
    )
    circuits.append(circ)

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

Ausgabe der vorherigen Code-Zelle

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

Ausgabe der vorherigen Code-Zelle

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

Ausgabe der vorherigen Code-Zelle

Sobald diese Circuits konstruiert wurden, kannst du sie transpilieren und mit einem QPU ausführen. Für dieses Beispiel verwenden wir einen der rauschfreien Simulatoren, um zu zeigen, wie der Trotter-Fehler reduziert wird.

backend = GenericBackendV2(num_qubits=10)
transpiler = generate_preset_pass_manager(
    optimization_level=2, backend=backend
)

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

estimator = StatevectorEstimator()
job = estimator.run([(circ, observable) for circ in transpiled_circuits])
result = job.result()

mpf_evs = [res.data.evs for res in result]
print(mpf_evs)

[array(0.23799162), array(0.35754312), array(0.38649906)]

Ergebnisse rekonstruieren

Nachdem die Circuits ausgeführt wurden, ist die Rekonstruktion der Ergebnisse recht unkompliziert. Wie auf der MPF-Übersichtsseite beschrieben, wird unsere Observable durch die gewichtete Summe rekonstruiert:

$\langle M \rangle = \sum_j x_j \langle M_j \rangle.$

wobei $x_j$ die ermittelten Koeffizienten und $\langle M_j \rangle$ die Schätzung der Observable $\sum_i \langle Z_i \rangle$ für jeden Circuit sind. Anschließend können wir die erhaltenen Ergebnisse mit dem exakten Wert mithilfe des Pakets scipy.linalg vergleichen.

exp_H = expm(-1j * time * hamiltonian.to_matrix())
initial_state = np.zeros(exp_H.shape[0])
initial_state[0] = 1.0

time_evolved_state = exp_H @ initial_state
exact_obs = (
    time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
)

# Print out the different observable measurements
print(f"Exact value: {exact_obs.real}")
print(f"PF with 19 steps: {mpf_evs[-1]}")
print(f"MPF using exact solution: {mpf_evs @ coeffs_exact.value}")
print(f"MPF using approximate solution: {mpf_evs @ coeffs_approx.value}")

Exact value: 0.4006024248789992
PF with 19 steps: 0.3864990619977402
MPF using exact solution: 0.3954847855979902
MPF using approximate solution: 0.4299121425348959

Hier ist zu sehen, dass die MPF den Trotter-Fehler im Vergleich zu dem reduziert hat, der mit einer einzigen PF mit $k_j=19$ erzielt wurde. Das Näherungsmodell lieferte jedoch einen Erwartungswert, der schlechter als das exakte Modell ist. Dies verdeutlicht, wie wichtig es ist, enge Konvergenzkriterien für das Näherungsmodell zu verwenden und eine „gute" Menge von $k_j$ -Werten zu finden.