Redukce Trotterovy chyby hamiltonovy dynamiky pomocí multi-produktových formulí

V tomto notebooku se naučíš, jak použít Multi-Product Formula (MPF) k dosažení nižší Trotterovy chyby na naší observabli ve srovnání s chybou způsobenou nejhlubším Trotterovým Circuit em, který skutečně spustíme. Budeš to dělat tak, že projdeš kroky Qiskit vzoru:

• Krok 1: Mapování na kvantový problém
  • Inicializace Hamiltoniánu našeho problému
  • Použití MPF ke generování Trotterizovaných Circuit ů časové evoluce
• Krok 2: Optimalizace problému
  • Zde transpilujeme naše Circuit y pro GenericBackendV2
• Krok 3: Spuštění experimentů
  • Použití StatevectorEstimator pro jednoduchost v tomto notebooku
• Krok 4: Rekonstrukce výsledků
  • Výpočet střední hodnoty MPF

Krok 1: Mapování na kvantový problém

1a: Nastavení Hamiltoniánu

Používáme Isingův model na řadě 10 uzlů:

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

kde $J$ je síla vazby mezi dvěma uzly a $h$ je vnější magnetické pole. Balíček qiskit_addon_utils poskytuje znovupoužitelné funkce pro různé účely.

Jeho modul qiskit_addon_utils.problem_generators poskytuje funkce pro generování Heisenbergových Hamiltonianů na daném grafu propojení. Tento graf může být buď rustworkx.PyGraph nebo CouplingMap, což usnadňuje použití v pracovních postupech zaměřených na Qiskit.

V následujícím příkladu vytvoříme jednoduchou řadu 10 Qubitů pomocí 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")

Dále vygenerujeme SparsePauliOp na daném propojení s požadovanými konstantami.

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

Observabla, kterou budeme měřit, je celková magnetizace, kterou jednoduše sestrojíme takto:

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: Multi-Product Formulas

MPF snižují Trotterovu chybu hamiltonovy dynamiky pomocí váženého součtu několika spuštění Circuit ů.

Abychom to konkretizovali, definujeme MPF jako:

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

kde $x_j$ jsou naše váhové koeficienty, $\rho^{k_j}_j$ je matice hustoty odpovídající čistému stavu získanému evolucí počátečního stavu pomocí produktové formule $S^{k_j}$ zahrnující $k_j$ Trotterových kroků, a $j$ indexuje počet produktových formulí tvořících MPF.

Klíčem je, že zbývající Trotterova chyba je menší než chyba, kterou by člověk získal pouhým použitím největší hodnoty $k_j$!

Užitečnost MPF lze nahlížet ze dvou perspektiv:

1. Pro pevný rozpočet Trotterových kroků, které jsi schopen spustit, lze získat výsledky s celkově menší Trotterovou chybou.
2. Pro počet Trotterových kroků vedoucí k hlubokým Circuit ům lze pomocí MPF najít několik Circuit ů s kratší hloubkou, jejichž spuštění povede k podobné Trotterově chybě.

Úvod do statických MPF

Statické MPF jsou ty, kde hodnoty $x_j$ NEZÁVISÍ na čase evoluce $t$.

Určení statických MPF koeficientů pro danou sadu hodnot $k_j$ spočívá v řešení soustavy lineárních rovnic: $Ax=b$, kde $x$ jsou naše hledané koeficienty, $A$ je matice závisející na $k_j$ a druhu PF, který používáme ($S$), a $b$ je vektor omezení. Pro stručnost zde nebudeme zacházet do dalších podrobností a místo toho tě odkážeme na dokumentaci LSE.

Řešení pro $x$ lze nalézt analyticky jako $x = A^{-1}b$, viz např. Carrera Vazquez et al., 2023 nebo Zhuk et al., 2023. Toto přesné řešení však může být „špatně podmíněné", což vede k velmi velkým L1-normám našich koeficientů $x$, a tím ke špatnému výkonu MPF. Místo toho lze také získat přibližné řešení, které minimalizuje L1-normu $x$, aby se pokusilo optimalizovat chování MPF.

V následujícím textu se naučíš, jak to vše provést.

Volba $k_j$

Volba $k_j$ je na koncovém uživateli. V principu lze zvolit libovolné hodnoty, ale některé hodnoty $k_j$ povedou k většímu zesílení šumu na reálných zařízeních než jiné. Proto je důležité snažit se najít „dobré" hodnoty $k_j$.

Zde jednoduše zvolíme několik pevných hodnot pro $k_j$. Nejmenší hodnota je motivována cílovým časem evoluce $t=8.0$, který nám obecně říká, že máme splnit $t/k_{\text{min}} \lt 1$, ale empiricky víme, že nastavení rovnosti $1$ obvykle také funguje. Pokud se chceš dozvědět více o tom, jak vybrat ostatní hodnoty $k_j$, přečti si příslušný návod: How to choose the Trotter steps for an MPF.

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

Nastavení LSE

Nyní, když jsme zvolili naše $k_j$, musíme nejprve zkonstruovat LSE, $Ax=b$, jak bylo vysvětleno výše. Matice $A$ závisí nejen na $k_j$, ale také na volbě produktové formule (PF) -- konkrétně na jejím řádu. Navíc lze zohlednit, zda je PF symetrická nebo ne (viz Carrera Vazquez et al., 2023), nastavením symmetric=True. To však není povinné, jak ukazuje Zhuk et al., 2023.

Zde budeme používat Suzuki-Trotterovu formuli druhého řádu, což dává order=2, a nastavíme 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.]))

Analytické řešení pro $x$

Jak bylo zmíněno, $x$ lze nalézt analyticky:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

Optimalizace $x$ pomocí přesného modelu

Alternativou k výpočtu $x=A^{-1}b$ je použití setup_exact_problem ke konstrukci instance cvxpy.Problem, která používá LSE jako omezení a jejíž optimální řešení dá $x$.

V další části bude jasné, proč toto rozhraní existuje.

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 ukazatel toho, zda MPF sestavené s těmito koeficienty přinese dobré výsledky, můžeme použít L1-normu (viz také Carrera Vazquez et al., 2023).

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

3.3889400921655914

Optimalizace $x$ pomocí aproximativního modelu

Může se stát, že L1-norma pro zvolenou sadu hodnot $k_j$ je považována za příliš vysokou. Pokud tomu tak je a nelze zvolit jinou sadu hodnot $k_j$, lze místo přesného řešení LSE použít aproximativní řešení.

K tomu jednoduše použij setup_sum_of_squares_problem ke konstrukci jiné instance cvxpy.Problem, která omezuje L1-normu na zvolenou prahovou hodnotu a zároveň minimalizuje rozdíl $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

Všimni si, že máš úplnou volnost v tom, jak tento optimalizační problém řešit, což znamená, že můžeš změnit optimalizační solver, jeho konvergenční prahy a podobně. Podívej se na příslušný návod How to use the approximate model.

1c: Nastavení Trotterových Circuit

V tuto chvíli jsme nalezli koeficienty rozkladu, $x$, a zbývá nám pouze vygenerovat Trotterizované kvantové Circuit. Opět nám k tomu poslouží modul qiskit_addon_utils.problem_generators:

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: Optimalizace problému

Obvykle se v tomto kroku vzoru optimalizují Circuit pro spuštění na hardwaru. Zde, protože používáme pouze simulátor bez šumu, jednoduše transpilujeme náš Circuit pro 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: Spuštění kvantových experimentů

Jak bylo vysvětleno na samém začátku, krok optimalizace 2 přeskočíme, protože očekávané hodnoty naší cílové observace jednoduše vypočítáme pomocí bezšumového simulátoru, konkrétně StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

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

Krok 4: Rekonstrukce výsledků

Nejprve extrahujeme jednotlivé očekávané hodnoty získané pro každý z Trotterových Circuit:

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

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

Dále je jednoduše zkombinujeme s koeficienty MPF, čímž získáme celkové očekávané hodnoty MPF. Níže to provedeme pro každý ze způsobů, jakými jsme vypočítali $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

Nakonec pro tento malý problém můžeme vypočítat přesnou referenční hodnotu pomocí scipy.linalg.expm takto:

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

Jasně vidíme, že MPF snížila Trotterovu chybu v porovnání s tou, kterou bychom získali při použití nejhlubšího jednotlivého PF s $k_j=19$. Zároveň vidíme, že přibližný model není bezchybný, protože ve skutečnosti vedl k očekávané hodnotě horší než přesné řešení. To ukazuje důležitost použití přísných konvergenčních kritérií u přibližného modelu, jak se dozvíš v průvodci How to use the approximate model.