Formuły wieloproduktowe do redukcji błędu Trottera

Szacowane użycie QPU: cztery minuty na procesorze Heron r2 (UWAGA: to tylko szacunek. Twój czas wykonania może się różnić.)

Podstawy

Ten samouczek pokazuje, jak używać Formuły Wieloproduktowej (MPF), aby uzyskać mniejszy błąd Trottera dla obserwabli w porównaniu z błędem generowanym przez najgłębszy obwód Trottera, który faktycznie wykonamy. MPF redukuje błąd Trottera dynamiki Hamiltonowskiej poprzez ważoną kombinację kilku wykonań obwodów. Rozważmy zadanie wyznaczenia wartości oczekiwanych obserwabli dla stanu kwantowego $\rho(t)=e^{-i H t} \rho(0) e^{i H t}$ z hamiltonianem $H$. Można użyć Formuł Produktowych (PF) do przybliżenia ewolucji czasowej $e^{-i H t}$ w następujący sposób:

• Zapisz Hamiltonian $H$ jako $H=\sum_{a=1}^d F_a,$ gdzie $F_a$ są operatorami hermitowskimi takimi, że każdy odpowiadający im operator unitarny może być efektywnie zaimplementowany na urządzeniu kwantowym.
• Przybliż wyrazy $F_a$, które nie komutują ze sobą.

Wtedy formuła produktowa pierwszego rzędu (formuła Lie-Trottera) to:

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

która ma kwadratowy człon błędu $S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right)$. Można też używać formuł produktowych wyższego rzędu (formuły Lie-Trotter-Suzuki), które zbiegają szybciej i są zdefiniowane rekurencyjnie jako:

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

gdzie $\chi$ jest rzędem symetrycznej PF, a $s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}$. Dla długich ewolucji czasowych można podzielić przedział czasu $t$ na $k$ podprzedziałów, zwanych krokami Trottera, o długości $t/k$ i przybliżyć ewolucję czasową w każdym podprzedziale formułą produktową rzędu $\chi$, $S_{\chi}$. Tym samym PF rzędu $\chi$ dla operatora ewolucji czasowej po $k$ krokach Trottera wynosi:

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

gdzie człon błędu maleje wraz z liczbą kroków Trottera $k$ i rzędem $\chi$ PF.

Dla danej liczby całkowitej $k \geq 1$ i formuły produktowej $S_{\chi}(t)$ przybliżony stan po ewolucji czasowej $\rho_k(t)$ można uzyskać z $\rho_0$ przez zastosowanie $k$ iteracji formuły produktowej $S_{\chi}\left(\frac{t}{k}\right)$.

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

$\rho_k(t)$ jest przybliżeniem $\rho(t)$ z błędem przybliżenia Trottera ||$\rho_k(t)-\rho(t) ||$. Jeśli rozważymy liniową kombinację przybliżeń Trottera $\rho(t)$:

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

gdzie $x_j$ są współczynnikami wagowymi, $\rho^{k_j}_j$ jest macierzą gęstości odpowiadającą czystemu stanowi uzyskanemu przez ewolucję stanu początkowego formułą produktową $S^{k_j}_{\chi}$ z $k_j$ krokami Trottera, a $j \in {1, ..., l}$ indeksuje liczbę PF składających się na MPF. Wszystkie wyrazy w $\mu(t)$ używają tej samej formuły produktowej $S_{\chi}(t)$ jako podstawy. Celem jest poprawa ||$\rho_k(t)-\rho(t) \|$ przez znalezienie $\mu(t)$ z jeszcze mniejszym $\|\mu(t)-\rho(t)\|$.

• $\mu(t)$ nie musi być stanem fizycznym, ponieważ $x_i$ nie musi być dodatnie. Celem jest minimalizacja błędu wartości oczekiwanej obserwabli, a nie znalezienie fizycznego zastępstwa dla $\rho(t)$.
• $k_j$ wyznacza zarówno głębokość obwodu, jak i poziom przybliżenia Trottera. Mniejsze wartości $k_j$ prowadzą do krótszych obwodów, które mają mniej błędów, ale są mniej dokładnym przybliżeniem pożądanego stanu.

Kluczem jest to, że pozostały błąd Trottera dany przez $\mu(t)$ jest mniejszy niż błąd Trottera, który uzyskano by po prostu używając największej wartości $k_j$.

Możesz spojrzeć na użyteczność tej metody z dwóch perspektyw:

1. Przy ustalonym budżecie kroków Trottera, które możesz wykonać, możesz uzyskać wyniki z błędem Trottera mniejszym w sumie.
2. Dla pewnej docelowej liczby kroków Trottera, która jest zbyt duża do wykonania, możesz użyć MPF, aby znaleźć zestaw obwodów o mniejszej głębokości, których uruchomienie prowadzi do podobnego błędu Trottera.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane:

• Qiskit SDK v1.0 lub nowszy, z obsługą wizualizacji
• Qiskit Runtime v0.22 lub nowszy (pip install qiskit-ibm-runtime)
• MPF Qiskit addons (pip install qiskit_addon_mpf)
• Qiskit addons utils (pip install qiskit_addon_utils)
• Biblioteka Quimb (pip install quimb)
• Biblioteka Qiskit Quimb (pip install qiskit-quimb)
• Numpy v0.21 dla zgodności między pakietami (pip install numpy==0.21)

Część I. Przykład małej skali

Badanie stabilności MPF

Nie ma oczywistych ograniczeń co do wyboru liczby kroków Trottera $k_j$ składających się na stan MPF $\mu(t)$. Jednak należy je dobierać ostrożnie, aby uniknąć niestabilności wynikowych wartości oczekiwanych obliczonych z $\mu(t)$. Dobrą ogólną zasadą jest ustawienie najmniejszego kroku Trottera $k_{\text{min}}$ tak, aby $t/k_{\text{min}} \lt 1$. Jeśli chcesz dowiedzieć się więcej o tym i jak wybierać pozostałe wartości $k_j$, zapoznaj się z przewodnikiem Jak wybrać kroki Trottera dla MPF.

W poniższym przykładzie badamy stabilność rozwiązania MPF, obliczając wartość oczekiwaną magnetyzacji dla zakresu czasów przy użyciu różnych stanów po ewolucji czasowej. Konkretnie porównujemy wartości oczekiwane obliczone z każdej przybliżonej ewolucji czasowej zaimplementowanej z odpowiednimi krokami Trottera oraz różnymi modelami MPF (współczynniki statyczne i dynamiczne) z dokładnymi wartościami obserwabli po ewolucji czasowej. Najpierw zdefiniujmy parametry dla formuł Trottera i czasów ewolucji

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

W tym przykładzie użyjemy stanu Neela jako stanu początkowego $\vert \text{Neel} \rangle = \vert 0101...01 \rangle$ oraz modelu Heisenberga na łańcuchu 10 węzłów dla Hamiltonianu rządzącego ewolucją czasową

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

gdzie $J$ jest stałą sprzężenia dla krawędzi najbliższych sąsiadów.

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

Obserwabla, którą będziemy mierzyć, to magnetyzacja pary qubitów w środku łańcucha.

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

Definiujemy przebieg transpilera zbierający rotacje XX i YY w obwodzie jako pojedynczą bramkę XX+YY. Pozwoli to nam wykorzystać właściwości zachowania spinu TeNPy podczas obliczania MPO, co znacznie przyspieszy obliczenia.

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

Następnie tworzymy obwody implementujące przybliżone ewolucje czasowe Trottera.

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)

Następnie obliczamy wartości oczekiwane po ewolucji czasowej z obwodów Trottera.

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)

Obliczamy też dokładne wartości oczekiwane do porównania.

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)

Statyczne współczynniki MPF

Statyczne MPF to takie, w których wartości $x_j$ nie zależą od czasu ewolucji $t$. Rozważmy PF rzędu $\chi = 1$ z $k_j$ krokami Trottera, co można zapisać jako:

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

gdzie $A_n$ są macierzami zależnymi od komutatorów wyrazów $F_a$ w dekompozycji Hamiltonianu. Ważne jest, że $A_n$ same w sobie są niezależne od czasu i liczby kroków Trottera $k_j$. Dlatego możliwe jest anulowanie wyrazów błędu niższego rzędu wnoszących wkład do $\mu(t)$ przez staranny dobór wag $x_j$ kombinacji liniowej. Aby anulować błąd Trottera dla pierwszych $l-1$ wyrazów (które dają największy wkład, bo odpowiadają mniejszej liczbie kroków Trottera) w wyrażeniu na $\mu(t)$, współczynniki $x_j$ muszą spełniać następujące równania:

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

z $n=0, ... l-2$. Pierwsze równanie gwarantuje brak obciążenia w konstruowanym stanie $\mu(t)$, drugie zapewnia anulowanie błędów Trottera. Dla PF wyższego rzędu drugie równanie ma postać $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0$, gdzie $\eta = \chi + 2n$ dla symetrycznych PF i $\eta = \chi + n$ w przeciwnym razie, z $n=0, ..., l-2$. Wynikowy błąd (Ref. [1], [2]) wynosi wtedy

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

Wyznaczenie statycznych współczynników MPF dla danego zestawu wartości $k_j$ sprowadza się do rozwiązania układu równań liniowych zdefiniowanego przez powyższe dwa równania względem zmiennych $x_j$: $Ax=b$. Gdzie $x$ są szukanymi współczynnikami, $A$ jest macierzą zależną od $k_j$ i typu używanej PF ($S$), a $b$ jest wektorem ograniczeń. Konkretnie:

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

gdzie $\chi$ to order, $s$ wynosi $2$ jeśli symmetric jest True, a $1$ w przeciwnym razie, $k_{j}$ to trotter_steps, a $x$ to zmienne do wyznaczenia. Indeksy $i$ i $j$ zaczynają się od $0$. Można to też zobrazować w postaci macierzowej:

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

oraz

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

Więcej szczegółów znajdziesz w dokumentacji Układu Równań Liniowych (LSE).

Możemy znaleźć rozwiązanie dla $x$ analitycznie jako $x = A^{-1}b$; patrz np. Ref. [1] lub [2]. Jednak to dokładne rozwiązanie może być „źle uwarunkowane", co prowadzi do bardzo dużych norm L1 naszych współczynników $x$, co może skutkować słabą wydajnością MPF. Zamiast tego można uzyskać przybliżone rozwiązanie minimalizujące normę L1 parametru $x$, aby spróbować zoptymalizować zachowanie MPF.

Konfiguracja LSE

Teraz, gdy wybraliśmy wartości $k_j$, musimy najpierw skonstruować LSE, $Ax=b$ zgodnie z powyższym opisem. Macierz $A$ zależy nie tylko od $k_j$, ale też od naszego wyboru PF, w szczególności jej rzędu. Dodatkowo możesz uwzględnić, czy PF jest symetryczna czy nie (patrz [1]), ustawiając symmetric=True/False. Nie jest to jednak wymagane, jak pokazuje Ref. [2].

from qiskit_addon_mpf.static import setup_static_lse

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

Przeanalizujmy wartości wybrane powyżej, aby skonstruować macierz $A$ i wektor $b$. Dla $j=0,1, 2$ kroków Trottera $k_j = [1, 2, 4]$, rzędu $\chi = 2$ i wyboru niesymetrycznych kroków Trottera ($s=1$), elementy macierzy $A$ poniżej pierwszego wiersza wyznacza wyrażenie $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$, konkretnie:

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

lub w postaci macierzowej:

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

Można to sprawdzić, inspekcjonując obiekt lse:

lse.A

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

Natomiast wektor ograniczeń $b$ ma następujące elementy: $b_{0} = 1$ $b_1 = b_2 = 0$

Zatem

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

I analogicznie w lse:

lse.b

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

Obiekt lse ma metody do wyznaczania statycznych współczynników $x_j$ spełniających układ równań.

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]

Optymalizacja $x$ za pomocą modelu dokładnego

Alternatywnie do obliczania $x=A^{-1}b$, możesz też użyć setup_exact_model, aby skonstruować instancję cvxpy.Problem, która używa LSE jako ograniczeń, a jej optymalne rozwiązanie daje $x$.

W następnej sekcji stanie się jasne, dlaczego ten interfejs istnieje.

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]

Jako wskaźnik tego, czy MPF skonstruowane z tymi współczynnikami da dobre wyniki, możemy użyć normy L1 (patrz też 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

Optymalizacja $x$ za pomocą modelu przybliżonego

Może się zdarzyć, że norma L1 dla wybranego zestawu wartości $k_j$ zostanie uznana za zbyt wysoką. Jeśli tak jest i nie możesz wybrać innego zestawu wartości $k_j$, możesz użyć przybliżonego rozwiązania LSE zamiast dokładnego.

W tym celu wystarczy użyć setup_approximate_model, aby skonstruować inną instancję cvxpy.Problem, która ogranicza normę L1 do wybranego progu, minimalizując jednocześnie różnicę $Ax$ i $b$.

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

Zauważ, że masz pełną swobodę w sposobie rozwiązywania tego problemu optymalizacyjnego, co oznacza, że możesz zmieniać solver optymalizacyjny, jego progi zbieżności i tak dalej. Sprawdź odpowiedni przewodnik na temat Jak używać modelu przybliżonego.

Dynamiczne współczynniki MPF

W poprzedniej sekcji wprowadziliśmy statyczne MPF poprawiające standardowe przybliżenie Trottera. Jednak ta statyczna wersja niekoniecznie minimalizuje błąd przybliżenia. Konkretnie, statyczne MPF, oznaczone $\mu^S(t)$, nie jest optymalną projekcją $\rho(t)$ na podprzestrzeń rozpiętą przez stany formuły produktowej $\{\rho_{k_i}(t)\}_{i=1}^r$.

Aby temu zaradzić, rozważamy dynamiczne MPF (wprowadzone w Ref. [2] i eksperymentalnie zademonstrowane w Ref. [3]), które minimalizuje błąd przybliżenia w normie Frobeniusa. Formalnie koncentrujemy się na minimalizacji

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

względem pewnych współczynników $x_i(t)$ w każdej chwili $t$. Optymalny projektor w normie Frobeniusa to $\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t)$, i nazywamy $\mu^D(t)$ dynamicznym MPF. Podstawiając powyższe definicje:

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

gdzie $M_{i,j}(t)$ jest macierzą Grama, zdefiniowaną przez

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

oraz

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

reprezentuje nakładanie się dokładnego stanu $\rho(t)$ z każdym przybliżeniem formuły produktowej $\rho_{k_i}(t)$. W praktycznych scenariuszach te nakładania można mierzyć tylko w przybliżeniu ze względu na szum lub częściowy dostęp do $\rho(t)$.

Tutaj $\lvert\psi_{\mathrm{in}}\rangle$ jest stanem początkowym, a $S(\cdot)$ to operacja stosowana w formule produktowej. Wybierając współczynniki $x_i(t)$ minimalizujące to wyrażenie (i obsługując przybliżone dane o nakładaniu, gdy $\rho(t)$ nie jest w pełni znany), uzyskujemy „najlepsze" (w sensie normy Frobeniusa) dynamiczne przybliżenie $\rho(t)$ w podprzestrzeni MPF. Wielkości $L_i(t)$ i $M_{i,j}(t)$ można obliczać efektywnie metodami sieci tensorowych [3]. Addon MPF dla Qiskit udostępnia kilka „backendów" do tego obliczenia. Poniższy przykład pokazuje najbardziej elastyczny sposób jego realizacji, a dokumentacja backendu warstwowego TeNPy szczegółowo to wyjaśnia. Aby użyć tej metody, zacznij od obwodu implementującego pożądaną ewolucję czasową i utwórz modele reprezentujące te operacje z warstw odpowiedniego obwodu. Na końcu tworzony jest obiekt Evolver, który można użyć do wygenerowania wielkości $M_{i,j}(t)$ i $L_i(t)$ po ewolucji czasowej. Zaczynamy od stworzenia obiektu Evolver odpowiadającego przybliżonej ewolucji czasowej (ApproxEvolverFactory) zaimplementowanej przez obwody. Zwróć szczególną uwagę na zmienną order, aby wartości się zgadzały. Pamiętaj, że przy generowaniu obwodów odpowiadających przybliżonej ewolucji czasowej używamy wartości zastępczych time = 1.0 i liczby kroków Trottera (reps=1). Właściwe obwody przybliżające są następnie generowane przez solver dynamicznego problemu w setup_dynamic_lse.

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

ostrzeżenie

Opcje LayerwiseEvolver określające szczegóły symulacji sieci tensorowej muszą być dobrane ostrożnie, aby uniknąć sformułowania źle zdefiniowanego problemu optymalizacyjnego.

Następnie konfigurujemy dokładny evolver (na przykład ExactEvolverFactory), który zwraca obiekt Evolver obliczający prawdziwą lub „referencyjną" ewolucję czasową. W praktyce przybliżalibyśmy dokładną ewolucję, używając formuły Suzuki-Trottera wyższego rzędu lub innej wiarygodnej metody z małym krokiem czasowym. Poniżej przybliżamy dokładny stan po ewolucji formułą Suzuki-Trottera czwartego rzędu z małym krokiem czasowym dt=0.1, co oznacza, że liczba kroków Trottera używana w chwili $t$ wynosi $k=t/dt$. Określamy też opcje obcinania specyficzne dla TeNPy, aby ograniczyć maksymalny wymiar więzi bazowej sieci tensorowej oraz minimalne wartości osobliwe rozciętych więzi sieci. Parametry te mogą wpływać na dokładność wartości oczekiwanej obliczonej za pomocą dynamicznych współczynników MPF, dlatego ważne jest badanie zakresu wartości, aby znaleźć optymalny balans między czasem obliczeniowym a dokładnością. Pamiętaj, że obliczanie współczynników MPF nie polega na wartości oczekiwanej PF uzyskanej z wykonania sprzętowego, więc można je dostrajać w postprocesingu.

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

Następnie utwórz stan początkowy systemu w formacie zgodnym z TeNPy (np. MPS_neel_state=$\vert 0101...01 \rangle$). Konfiguruje to wielociałową funkcję falową, którą będziesz ewoluować w czasie $\lvert\psi_{\mathrm{in}}\rangle$ jako tensor.

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)

Dla każdego kroku czasowego $t$ konfigurujemy dynamiczny układ równań liniowych metodą setup_dynamic_lse. Odpowiadający mu obiekt zawiera informacje o dynamicznym problemie MPF: lse.A daje macierz Grama $M$, natomiast lse.b daje nakładanie $L$. Możemy wtedy rozwiązać LSE (gdy nie jest źle zdefiniowany), aby znaleźć dynamiczne współczynniki za pomocą setup_frobenius_problem. Ważne jest odróżnienie ich od współczynników statycznych, które zależą tylko od szczegółów używanej formuły produktowej i są niezależne od szczegółów ewolucji czasowej (Hamiltoniana i stanu początkowego).

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

Na koniec wykreśl te wartości oczekiwane w funkcji czasu ewolucji.

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

W przypadkach takich jak powyższy przykład, gdzie PF $k=1$ zachowuje się słabo przez cały czas, jakość wyników dynamicznego MPF jest również mocno zaburzona. W takich sytuacjach warto zbadać możliwość użycia poszczególnych PF z wyższą liczbą kroków Trottera, aby poprawić ogólną jakość wyników. W tych symulacjach obserwujemy współzależność różnych rodzajów błędów: błędu z skończonego próbkowania oraz błędu Trottera z formuł produktowych. MPF pomaga redukować błąd Trottera pochodzący z formuł produktowych, ale wiąże się z wyższym błędem próbkowania w porównaniu z formułami produktowymi. Może to być korzystne, ponieważ formuły produktowe mogą redukować błąd próbkowania zwiększonym próbkowaniem, ale systematyczny błąd z przybliżenia Trottera pozostaje niezmieniony.

Innym interesującym zachowaniem, które możemy zaobserwować na wykresie, jest to, że wartość oczekiwana PF dla $k=1$ zaczyna zachowywać się chaotycznie (poza tym, że nie jest dobrym przybliżeniem dokładnej) w czasach, dla których $t/k > 1$, zgodnie z wyjaśnieniami w przewodniku dotyczącym wyboru liczby kroków Trottera.

Krok 1: Odwzorowanie danych wejściowych na problem kwantowy

Rozważmy teraz pojedynczy czas $t=1.0$ i obliczmy wartość oczekiwaną magnetyzacji różnymi metodami przy użyciu jednego QPU. Konkretny wybór $t$ został dokonany tak, aby zmaksymalizować różnicę między poszczególnymi metodami i zaobserwować ich względną skuteczność. Aby wyznaczyć przedział czasu, dla którego dynamiczny MPF gwarantuje produkcję obserwowalnych o mniejszym błędzie niż którakolwiek z indywidualnych formuł Trottera wchodzących w skład multi-produktu, możemy zastosować „test MPF" — patrz równanie (17) i otaczający tekst w [3].

Przygotowanie obwodów Trottera

W tym momencie znaleźliśmy współczynniki rozwinięcia $x$ i pozostaje tylko wygenerować kwantowe obwody Trottera. Moduł qiskit_addon_utils.problem_generators ponownie przychodzi z pomocą, oferując użyteczną funkcję do tego celu:

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)

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Wróćmy do obliczania wartości oczekiwanej dla pojedynczego punktu czasowego. Wybierzemy Backend do wykonania eksperymentu na sprzęcie.

from qiskit_ibm_runtime import QiskitRuntimeService

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

qubits = list(range(backend.num_qubits))

Następnie usuwamy odstające qubity z mapy sprzężeń, aby zapewnić, że etap układania Transpilera nie uwzględni ich. Poniżej korzystamy z zapisanych właściwości Backendu przechowywanych w obiekcie target i usuwamy qubity, których błąd pomiaru lub bramki dwu-qubitowej przekracza określony próg (max_meas_err, max_twoq_err) albo czas $T_2$ (określający utratę koherencji) jest poniżej określonego progu (min_t2).

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

Chcemy ustawić max_meas_err, min_t2 oraz max_twoq_err tak, aby znaleźć wystarczająco duży podzbiór qubitów obsługujący uruchamiany obwód. W naszym przypadku wystarczy znaleźć jednowymiarowy łańcuch 10 qubitów.

cust_cmap.draw()

Możemy następnie odwzorować obwód i obserwowalną na fizyczne qubity urządzenia.

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

Krok 3: Wykonanie przy użyciu prymitywów Qiskit

Przy użyciu prymitywu Estimator możemy uzyskać oszacowanie wartości oczekiwanej z QPU. Wykonujemy zoptymalizowane obwody AQC z dodatkowymi technikami łagodzenia i tłumienia błędów.

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

Krok 4: Post-processing i zwrócenie wyniku w pożądanym formacie klasycznym

Jedynym krokiem post-processingu jest połączenie wartości oczekiwanych uzyskanych z prymitywów Qiskit Runtime przy różnych liczbach kroków Trottera przy użyciu odpowiednich współczynników MPF. Dla obserwowalnej $A$ mamy:

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

Najpierw wyodrębniamy indywidualne wartości oczekiwane uzyskane dla każdego z obwodów Trottera:

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

Następnie łączymy je ze współczynnikami MPF, aby uzyskać całkowite wartości oczekiwane MPF. Poniżej robimy to dla każdego z różnych sposobów, przy użyciu których obliczyliśmy $x$.

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

Na koniec, dla tego małego problemu możemy obliczyć dokładną wartość referencyjną przy użyciu scipy.linalg.expm w następujący sposób:

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

W powyższym przykładzie metoda dynamicznego MPF wypada najlepiej pod względem wartości oczekiwanej, poprawiając wynik w porównaniu do tego, który uzyskalibyśmy stosując samą najwyższą liczbę kroków Trottera. Choć różne techniki MPF nie zawsze osiągają poprawioną wartość oczekiwaną w porównaniu z najwyższą liczbą kroków Trottera (jak modele dokładny i przybliżony na powyższym wykresie), odchylenia standardowe tych wartości dobrze oddają zwiększoną wariancję wynikającą ze stosowania techniki MPF. Podkreśla to niepewność wokół uzyskanej wartości oczekiwanej, która zawsze obejmuje wartość oczekiwaną, jakiej spodziewalibyśmy się po dokładnej ewolucji czasowej układu. Z drugiej strony wartości oczekiwane obliczone przy mniejszej liczbie kroków Trottera nie są w stanie uchwycić dokładnej wartości oczekiwanej w granicach swojej niepewności, zwracając z wysoką pewnością błędny wynik.

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

Część II: skalowanie problemu

Przeskalujmy problem poza zakres dokładnej symulacji. W tej sekcji skupimy się na odtworzeniu niektórych wyników przedstawionych w Ref. [3].

Krok 1: Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Hamiltonian

Dla przykładu na dużą skalę używamy modelu XXZ na linii 50 węzłów:

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

gdzie $J_{i,(i+1)}$ jest losowym współczynnikiem odpowiadającym krawędzi $(i, i+1)$. Jest to Hamiltonian rozważany w demonstracji przedstawionej w Ref. [3].

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

Jako obserwowalną wybieramy $Z_{24}Z_{25}$, zgodnie z dolnym panelem Rys. 5 z Ref. [3].

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

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

Wybór kroków Trottera

Eksperyment przedstawiony na Rys. 4 z Ref. [3] używa $k_j = [2, 3, 4]$ symetrycznych kroków Trottera rzędu $2$. Skupiamy się na wynikach dla czasu $t=3$, gdzie MPF i PF z większą liczbą kroków Trottera (w tym przypadku 6) mają ten sam błąd Trottera. Jednak wartość oczekiwana MPF jest obliczana z obwodów odpowiadających mniejszej liczbie kroków Trottera, a zatem jest płytsza. W praktyce, nawet jeśli MPF i obwód z głębszymi krokami Trottera mają ten sam błąd Trottera, spodziewamy się, że eksperymentalna wartość oczekiwana obliczona z obwodów MPF będzie bliższa wartości teoretycznej, ponieważ wymaga uruchomienia płytszych obwodów, które są mniej narażone na szumy sprzętowe w porównaniu z obwodem odpowiadającym PF z większą liczbą kroków.

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

Konfiguracja LSE

Przyjrzyjmy się tutaj statycznym współczynnikom MPF dla tego problemu.

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

Współczynniki dynamiczne

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

Konstruowanie każdego z obwodów Trottera w naszym rozkładzie MPF

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)

Konstruowanie obwodu Trottera o porównywalnym błędzie Trottera do MPF

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)

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

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
]

Krok 3: Wykonanie przy użyciu prymitywów Qiskit

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

Krok 4: Przetwarzanie końcowe i zwrócenie wyniku w żądanym formacie klasycznym

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

Wykonując obwody na sprzęcie, możemy napotkać dodatkowe trudności w uzyskaniu dokładnych wartości oczekiwanych ze względu na obecność szumów sprzętowych. Nie jest to uwzględnione w formalizmie MPF i może działać na niekorzyść rozwiązania MPF. Na przykład może to być przyczyną tego, że dynamiczne współczynniki nie zapewniają lepszego oszacowania wartości oczekiwanej w porównaniu z przybliżonym współczynnikiem statycznym na wykresie. Oznacza to, że przybliżony ewolutor, który symuluje przybliżony obwód, nie odzwierciedla dokładnie wyników uzyskanych przez wykonanie przybliżonych obwodów w obecności szumów sprzętowych. Z tych powodów zaleca się łączenie różnych technik łagodzenia błędów, aby uzyskać wyniki jak najbliższe wartościom idealnym dla każdej z formuł produktowych. Pozwoli to wykazać spójne korzyści płynące z podejścia MPF.

Ogólnie rzecz biorąc, przybliżone współczynniki statyczne wciąż dają dokładniejsze rozwiązanie niż formuła produktowa z większą liczbą kroków Trottera przy tej samej wielkości błędu Trottera w scenariuszu bez szumów.

Ważne jest również, aby zwrócić uwagę, że w przykładzie odtwarzającym eksperyment z Ref. [3] punkt czasowy $t=3$ wykracza poza granicę, przy której oczekuje się dobrego zachowania PF z $k=2$, co wynosi $t/k>1$, zgodnie z omówieniem w tym przewodniku.

Referencje

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