Redukcja błędu Trottera w dynamice hamiltonowskiej za pomocą formuł wieloproduktowych

W tym notebooku nauczysz się, jak używać Formuły Wieloproduktowej (MPF), aby uzyskać mniejszy błąd Trottera na naszej obserwabli w porównaniu z błędem wprowadzanym przez najgłębszy Circuit Trottera, który faktycznie wykonamy. Osiągniesz to, przechodząc przez kolejne kroki wzorca Qiskit:

• Krok 1: Odwzorowanie na problem kwantowy
  • Inicjalizacja Hamiltonianu naszego problemu
  • Użycie MPF do generowania obwodów ewolucji czasowej metodą Trottera
• Krok 2: Optymalizacja problemu
  • Tutaj transpilujemy nasze Circuit dla GenericBackendV2
• Krok 3: Wykonanie eksperymentów
  • Użycie StatevectorEstimator dla uproszczenia w tym notebooku
• Krok 4: Rekonstrukcja wyników
  • Obliczenie wartości oczekiwanej MPF

Krok 1: Odwzorowanie na problem kwantowy

1a: Konfiguracja Hamiltonianu

Używamy modelu Isinga na łańcuchu 10 węzłów:

$$\hat{\mathcal{H}}_{\text{Ising}} = \sum_{i=1}^{9} J_{i,(i+1)} Z_i Z_{(i+1)} + \sum_{i=1}^{10} h_i X_i \, ,$$

gdzie $J$ to siła sprzężenia między dwoma węzłami, a $h$ to zewnętrzne pole magnetyczne. Pakiet qiskit_addon_utils dostarcza pewne wielokrotnie używane funkcje do różnych celów.

Moduł qiskit_addon_utils.problem_generators udostępnia funkcje do generowania hamiltonianów podobnych do modelu Heisenberga na zadanym grafie połączeń. Ten graf może być albo rustworkx.PyGraph, albo CouplingMap, co ułatwia jego użycie w przepływach pracy zorientowanych na Qiskit.

Poniżej tworzymy prosty łańcuch 10 Qubitów przy użyciu metody CouplingMap.from_line.

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

from qiskit.transpiler import CouplingMap

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

from rustworkx.visualization import graphviz_draw

graphviz_draw(coupling_map.graph, method="circo")

Następnie generujemy SparsePauliOp na podanej topologii połączeń z żądanymi stałymi.

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

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

Obserwabla, którą będziemy mierzyć, to całkowite namagnesowanie, które możemy po prostu skonstruować jak pokazano poniżej:

from qiskit.quantum_info import SparsePauliOp

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list([("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L)
print(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])

1b: Formuły Wieloproduktowe

MPF redukują błąd Trottera w dynamice hamiltonowskiej poprzez ważoną kombinację kilku wykonań Circuit.

Aby to ująć bardziej konkretnie, definiujemy MPF jako:

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

gdzie $x_j$ to nasze współczynniki wagowe, $\rho^{k_j}_j$ to macierz gęstości odpowiadająca czystemu stanowi uzyskanemu przez ewolucję stanu początkowego za pomocą formuły iloczynowej, $S^{k_j}$, obejmującej $k_j$ kroków Trottera, a $j$ indeksuje liczbę formuł iloczynowych składających się na MPF.

Kluczowe jest to, że pozostały błąd Trottera jest mniejszy niż błąd Trottera, który uzyskano by, stosując po prostu największą wartość $k_j$!

Użyteczność MPF możesz rozważać z dwóch perspektyw:

1. Przy ustalonym budżecie kroków Trottera, które jesteś w stanie wykonać, możesz uzyskać wyniki z łącznie mniejszym błędem Trottera.
2. Przy liczbie kroków Trottera skutkującej głębokimi Circuit, możesz użyć MPF, aby znaleźć kilka Circuit o mniejszej głębokości do uruchomienia, które dają podobny błąd Trottera.

Wprowadzenie do statycznych MPF

Statyczne MPF to takie, w których wartości $x_j$ NIE zależą od czasu ewolucji, $t$.

Wyznaczenie statycznych współczynników MPF dla danego zbioru wartości $k_j$ sprowadza się do rozwiązania układu równań liniowych: $Ax=b$, gdzie $x$ to nasze poszukiwane współczynniki, $A$ to macierz zależna od $k_j$ i rodzaju używanej formuły iloczynowej (PF) -- w szczególności jej rzędu, a $b$ to wektor ograniczeń. Dla zwięzłości nie będziemy tutaj wchodzić w więcej szczegółów i zamiast tego odsyłamy do dokumentacji LSE.

Możemy znaleźć rozwiązanie analityczne dla $x$ jako $x = A^{-1}b$, zob. np. Carrera Vazquez et al., 2023 lub Zhuk et al., 2023. Jednak to dokładne rozwiązanie może być „źle uwarunkowane", co skutkuje bardzo dużymi normami L1 naszych współczynników $x$, co może prowadzić do słabej wydajności MPF. Zamiast tego można też uzyskać rozwiązanie przybliżone, które minimalizuje normę L1 $x$, aby spróbować zoptymalizować zachowanie MPF.

Poniżej nauczysz się, jak to wszystko zrobić.

Wybór $k_j$

Wybór $k_j$ należy do użytkownika końcowego. Zasadniczo można wybrać dowolne wartości, ale niektóre $k_j$ będą prowadzić do większego wzmocnienia szumu na rzeczywistych urządzeniach niż inne wybory. Dlatego ważne jest, aby próbować znaleźć „dobre" wartości $k_j$.

Tutaj po prostu wybierzemy kilka ustalonych wartości $k_j$. Najmniejsza wartość jest motywowana docelowym czasem ewolucji $t=8.0$, który normalnie wskazuje na spełnienie $t/k_{\text{min}} \lt 1$, ale empirycznie wiemy, że ustawienie go równego $1$ zazwyczaj też działa. Jeśli chcesz dowiedzieć się więcej na ten temat i jak dobierać inne wartości $k_j$, zapoznaj się z odpowiednim przewodnikiem: How to choose the Trotter steps for an MPF.

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

Konfiguracja LSE

Teraz, gdy wybraliśmy nasze $k_j$, musimy najpierw skonstruować LSE, $Ax=b$, jak wyjaśniono powyżej. Macierz $A$ zależy nie tylko od $k_j$, ale też od naszego wyboru formuły iloczynowej (PF) -- w szczególności jej rzędu. Dodatkowo można wziąć pod uwagę, czy PF jest symetryczna, czy nie (zob. Carrera Vazquez et al., 2023), ustawiając symmetric=True. Nie jest to jednak wymagane, jak pokazano przez Zhuk et al., 2023.

Tutaj użyjemy formuły Suzuki-Trotter drugiego rzędu, co daje order=2, i ustawimy symmetric=True.

from qiskit_addon_mpf.static import setup_static_lse

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

Analityczne wyznaczenie $x$

Jak wspomniano wcześniej, możemy znaleźć $x$ analitycznie:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

Optymalizacja $x$ przy użyciu modelu dokładnego

Alternatywnie do obliczenia $x=A^{-1}b$, możesz też użyć setup_exact_problem, aby skonstruować instancję cvxpy.Problem, która używa LSE jako ograniczeń, a której optymalne rozwiązanie da $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.17239057 -1.19447005  2.02207947]

Jako wskaźnik tego, czy MPF zbudowane z tymi współczynnikami da dobre wyniki, możemy użyć normy L1 (zob. również Carrera Vazquez et al., 2023).

print(np.linalg.norm(coeffs_exact.value, ord=1))

3.3889400921655914

Optymalizacja $x$ przy użyciu modelu przybliżonego

Może się zdarzyć, że norma L1 dla wybranego zbioru wartości $k_j$ zostanie uznana za zbyt wysoką. W takim przypadku, gdy 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 użyj po prostu setup_sum_of_squares_problem, aby skonstruować inną instancję cvxpy.Problem, która ogranicza normę L1 do wybranego progu, minimalizując jednocześnie różnicę między $Ax$ a $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=3.0)
model_approx.solve()
print(coeffs_approx.value)
print(np.linalg.norm(coeffs_approx.value, ord=1))

[-0.40454257  0.57553173  0.8290123 ]
1.8090865903790838

Zwróć uwagę, ż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. Zapoznaj się z odpowiednim przewodnikiem How to use the approximate model.

1c: Konfiguracja Circuit Trottera

W tym momencie znaleźliśmy nasze współczynniki rozwinięcia, $x$, i pozostało już tylko wygenerowanie kwantowych Circuit przy użyciu metody Trottera. Ponownie moduł qiskit_addon_utils.problem_generators przychodzi z pomocą, aby to zrobić:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit

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)

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

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

Krok 2: Optymalizacja problemu

Normalnie jest to krok we wzorcu, podczas którego optymalizujesz swoje Circuit pod kątem wykonania na sprzęcie. Tutaj, ponieważ używamy jedynie bezszumowego symulatora, po prostu transpilujemy nasz Circuit dla GenericBackendV2.

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import generate_preset_pass_manager

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

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

Krok 3: Wykonanie eksperymentów kwantowych

Jak wyjaśniono na samym początku, pominiemy krok optymalizacji 2, ponieważ po prostu obliczymy wartości oczekiwane naszej docelowej obserwabli przy użyciu symulatora bez szumu, czyli StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

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

Krok 4: Rekonstrukcja wyników

Najpierw wyodrębniamy poszczególne wartości oczekiwane uzyskane dla każdego z Circuit Trottera:

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

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

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

print("Analytical    solution:", evs @ coeffs_analytical)
print("Exact model   solution:", evs @ coeffs_exact.value)
print("Approx. model solution:", evs @ coeffs_approx.value)

Analytical    solution: 0.3954847855980006
Exact model   solution: 0.39548478559800204
Approx. model solution: 0.42991214253489807

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

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(exact_obs.real)

0.40060242487899755

Wyraźnie widać, że MPF zredukowało błąd Trottera w porównaniu z błędem uzyskanym przy użyciu najgłębszej indywidualnej PF z $k_j=19$. Widzimy jednak również, że model przybliżony nie jest doskonały, ponieważ w rzeczywistości dał wartość oczekiwaną gorszą niż dokładne rozwiązanie. Pokazuje to znaczenie stosowania ścisłych kryteriów zbieżności w modelu przybliżonym, czego nauczysz się w przewodniku How to use the approximate model.