Multi-product formulas pro snížení Trotterovy chyby

Odhadované využití: čtyři minuty na procesoru Heron r2 (POZNÁMKA: Jedná se pouze o odhad. Skutečná doba běhu se může lišit.)

Výsledky učení

• Jak vícesložkové vzorce (MPF) snižují Trotterovu chybu při simulaci Hamiltoniánu kombinací středních hodnot z více mělkých obvodů

• Kdy jsou MPF výhodné oproti standardním produktovým vzorcům a kdy nejsou tím správným nástrojem

• Jak vypočítat statické a dynamické koeficienty MPF pomocí balíčku qiskit_addon_mpf

• Jak spustit celý postup MPF na hardwaru IBM Quantum®, včetně transpilace, mitigace chyb a následného zpracování

Předpoklady

• Metody kompilace pro obvody simulace Hamiltoniánu — představuje Trotterovy (produktové vzorce) obvody v Qiskitu.

• Produktové vzorce v Qiskitu, konkrétně syntetizační třídy SuzukiTrotter a LieTrotter.

• Primitivy Qiskit a rozhraní Estimator.

Základy

Co jsou vícesložkové vzorce?

Při simulaci kvantových systémů na kvantovém počítači je ústředním úkolem aproximovat operátor časového vývoje $e^{-iHt}$ pro Hamiltonián $H$. Standardní přístup používá produktové vzorce (PF), známé také jako Trotter-Suzukiho rozklady. Ty rozkládají $H = \sum_{a=1}^d F_a$ na členy, jejichž jednotlivé unitáry $e^{-iF_a t}$ lze efektivně implementovat, a poté aproximují celou evoluci jako uspořádaný součin těchto jednodušších unitár.

Produktový vzorec prvního řádu (Lie-Trotter) je:

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

který zavádí kvadratickou chybu: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. Symetrické vzorce vyššího řádu $S_{2\chi}(t)$, kde $\chi$ značí řád symetrického produktového vzorce (viz odkaz [1]), konvergují rychleji jako $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, ale za cenu hlubších obvodů na krok.

Aby se snížila chyba při pevném řádu $\chi$, obvykle se celkový čas evoluce $t$ rozdělí na $k$ menších Trotterových kroků. Každý krok aproximuje $e^{-iHt/k}$ produktovým vzorcem a kroky se zřetězí:

$$e^{-iHt} \approx \left[S_{2\chi}(t/k)\right]^k.$$

Pro symetrický vzorec $2\chi$-tého řádu pak zbytková Trotterova chyba škáluje jako $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Zvyšování $k$ tedy rychle potlačuje Trotterovu chybu — ale také lineárně prohlubuje obvod, a na zašuměném hardwaru to znamená více nahromaděného šumu hradel. Toto napětí mezi Trotterovou chybou (upřednostňuje větší $k$) a šumem hardwaru (upřednostňuje menší $k$) je přesně to, co mají vícesložkové vzorce vyřešit. Všimni si, že MPF kombinují výsledky z různých voleb $k$ při pevném řádu $\chi$ — neměnní řád podkladového produktového vzorce.

Vícesložkové vzorce (MPF) [1] sestavují váženou lineární kombinaci středních hodnot získaných z několika mělčích Trotterových obvodů, z nichž každý používá jiný počet Trotterových kroků $k_1, k_2, \ldots, k_r$ (množinu $r$ počtů kroků):

$$\langle A \rangle_{\text{MPF}}(t) = \sum_{j=1}^r x_j \, \langle A \rangle_{k_j}(t),$$

kde $\langle A \rangle_{k_j}(t)$ je očekávaná hodnota pozorovatelné $A$ v čase $t$ odhadnutá z Trotterova Circuit s $k_j$ kroky a koeficienty $\{x_j\}_{j=1}^r$ jsou zvoleny tak, aby se v kombinaci vzájemně zrušily vedoucí členy Trotterovy chyby. K tomuto výrazu se vrátíme v kroku 4, kde jej explicitně vyhodnotíme, abychom zkombinovali naše Trotterovy výsledky. Klíčové praktické zjištění je, že nejhlubší Circuit v MPF potřebuje pouze $k_{\max}$ kroků, což je mnohem méně než jediné $k$, které by bylo potřeba k dosažení stejné efektivní Trotterovy chyby přímo. Mělčí Circuity dělají přístup MPF vhodnějším pro zašuměný hardware.

Jak se určují koeficienty?

Existují dvě rodiny koeficientů MPF:

Statické koeficienty nezávisí na Hamiltoniánu, počátečním stavu ani na čase evoluce. Zjišťují se řešením lineární soustavy $Ax = b$, která vynucuje zrušení vedoucích členů Trotterovy chyby. Pro množinu Trotterových kroků $\{k_j\}_{j=1}^r$ použitých s symetrickým součinovým vzorcem řádu $2\chi$ vede rozvoj Trotterovy chyby v inverzních mocninách $k_j$ k omezujícím rovnicím tvaru:

$$\sum_{j=1}^r x_j = 1, \quad \sum_{j=1}^r \frac{x_j}{k_j^{\eta_n}} = 0 \quad (n = 0, \ldots, r-2),$$

kde celočíselné exponenty $\{\eta_n\}$ jsou řády postupných členů Trotterovy chyby pro zvolený součinový vzorec. Pro symetrický PF řádu $2\chi$ škáluje vedoucí chyba v $\left[S_{2\chi}(t/k)\right]^k$ jako $1/k^{2\chi}$, s dalšími korekcemi na $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — exponenty jsou tedy $\eta_n = 2\chi + 2n$. Pro nesymetrické PF přispívají jak liché, tak sudé mocniny a $\eta_n = 2\chi + n$. Úplné odvození viz odkaz [1]. První rovnice v soustavě výše zajišťuje nezkreslenost (MPF v limitě $k_j \to \infty$ reprodukuje přesnou očekávanou hodnotu) a zbývajících $r-1$ rovnic postupně ruší prvních $r-1$ členů Trotterovy chyby. Když je výsledná $L_1$-norma $\|x\|_1$ příliš velká (což zesiluje šum ze vzorkování), můžeš místo toho vyřešit přibližnou optimalizaci, která omezí $\|x\|_1$ a zároveň minimalizuje $\|Ax - b\|$.

Dynamické koeficienty [2], [3] navíc závisí na Hamiltoniánu, počátečním stavu a čase evoluce $t$. Minimalizují vzdálenost ve Frobeniově normě mezi skutečným časově vyvíjeným stavem a MPF aproximací:

$$\|\rho(t) - \mu^D(t)\|_F^2 = 1 + \sum_{i,j} M_{ij}(t)\, x_i(t)\, x_j(t) - 2\sum_i L_i(t)\, x_i(t),$$

kde $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ je Gramova matice překryvů mezi Trotterem vyvinutými stavy pro různé počty kroků $k_i, k_j$ a $L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ měří překryv s (přibližně) přesným stavem. V tomto tutoriálu jsou tyto veličiny počítány efektivně pomocí metod tenzorových sítí, konkrétně backendů založených na TeNPy v qiskit_addon_mpf.

Kdy použít MPF

MPF jsou nejvíce prospěšné, když:

• Hloubka Circuit je omezujícím faktorem. Pokud hardwarový šum omezuje hloubku, kterou můžeš spustit, použij MPF k dosažení vyšší efektivní přesnosti Trottera z mělčích Circuitů.

• Potřebuješ přesné očekávané hodnoty, ne úplnou přípravu stavu. MPF pracují na úrovni očekávaných hodnot — kombinují klasická čísla, ne kvantové stavy. Jsou proto ideální pro odhad pozorovatelných při použití primitiva Estimator.

• Kombinuješ mírný počet hodnot Trotterových kroků. Obvykle stačí zkombinovat $r = 3$–$5$ různých počtů kroků $k_j$, aby se zrušilo několik vedoucích členů Trotterovy chyby při zachování zvladatelné hodnoty $\|x\|_1$.

Kdy MPF nemusí pomoci

• Velmi krátké časy evoluce. Když je $t$ dostatečně malé na to, aby byl přesný už jediný Trotterův vzorec nízkého řádu, je režie spojená se spouštěním více Circuitů zbytečná.

• Úlohy přípravy stavu. MPF vytvářejí opravenou očekávanou hodnotu, ne opravený kvantový stav. Pokud potřebuješ skutečný časově vyvinutý stav (například jako vstup do jiné kvantové subrutiny), MPF se nepoužijí.

• Počty Trotterových kroků, které porušují režim konvergence. Odvození statických koeficientů rozvíjí každý jednotlivý výraz $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ jako řadu v $t/k_j$; tento rozvoj dobře konverguje pouze tehdy, když $t/k_{\min} \lesssim 1$. Pokud je $k_{\min}$ zvoleno příliš malé pro dané $t$, nejmělčí Circuit se ocitá daleko mimo perturbativní režim, členy chyby vyššího řádu, které MPF nezruší, se stanou velkými a zrušení může vyžadovat velké koeficienty. $L_1$-norma $\|x\|_1$ je praktickou diagnostikou: když $\|x\|_1 \gg 1$, může režie vzorkování $\propto \|x\|_1^2$ převážit nad snížením Trotterovy chyby. Podrobnosti najdeš v průvodci výběrem Trotterových kroků.

Co tento tutoriál pokrývá

Tento tutoriál provádí kompletním pracovním postupem MPF ve dvou fázích. Nejprve příklad v malém měřítku pomocí simulátoru (10qubitový Heisenbergův řetězec) ukazuje, jak nastavit úlohu, vypočítat statické a dynamické koeficienty MPF a porovnat výsledné očekávané hodnoty s přesnou diagonalizací. Poté příklad ve velkém měřítku na hardwaru (50qubitový XXZ řetězec) ukazuje, jak provést transpilaci, spuštění na hardwaru IBM Quantum s potlačením chyb a následné zpracování výsledků pomocí koeficientů MPF. V celém tutoriálu používáme balíček qiskit_addon_mpf spolu se standardními nástroji Qiskit.

Požadavky

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalováno následující:

• Qiskit SDK v2.0 nebo novější, s podporou vizualizace

• Qiskit Runtime v0.22 nebo novější (pip install qiskit-ibm-runtime)

• Simulátor Qiskit Aer (pip install qiskit-aer)

• MPF Qiskit addon s backendem TeNPy (pip install "qiskit-addon-mpf[tenpy]")

• Qiskit addon utilities (pip install qiskit-addon-utils)

• SciPy (pip install scipy)

Nastavení

Níže shromažďujeme všechny importy balíčků použité v tomto tutoriálu do jedné buňky. Také definujeme transpiler pass CollectAndCollapse, který slučuje sousední rotace rxx a ryy do jediného XXPlusYYGate. Tento pass se používá jak během konstrukce Circuit v kroku 1 (aby se udržel nízký počet gate), tak nepřímo, když v kroku 4 extrahujeme vrstvenou strukturu pro dynamický MPF (TeNPy očekává dvouqubitové gate, ne dvojice nesloučených rotací).

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

import warnings

import numpy as np
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)

# --- Helper: collect XX + YY rotations into a single gate ---
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))

Příklad v malém měřítku pomocí simulátoru

Krok 1: Mapování klasických vstupů na kvantový problém

Začínáme s 10qubitovým Heisenbergovým modelem na přímce, přičemž jako počáteční stav používáme Néelův stav $\vert 0101\ldots01 \rangle$. Hamiltonián je:

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

kde $J$ je síla vazby mezi nejbližšími sousedy. Měříme ZZ korelátor $Z_{L/2-1} Z_{L/2}$ na dvojici qubitů uprostřed řetězce a používáme Trotterovy kroky $k_j = [1, 2, 4]$ se součinovým vzorcem druhého řádu.

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

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

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

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

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

Sestavení Trotterových Circuit

Vytváříme Circuity implementující přibližné Trotterovy časové evoluce pro každý časový bod a každý počet Trotterových kroků. Pass CollectAndCollapse definovaný v sekci Nastavení shromažďuje rotace XX a YY do jednotlivých gate XX+YY, aby se připravil na efektivnější simulaci pomocí tenzorových sítí později.

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

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru

Pro příklad v malém měřítku cílíme na simulátor Aer. Před tím, než jsou Circuity připraveny ke spuštění, proběhnou dvě transformace:

1. Slučování gate na úrovni simulace Hamiltoniánu. V buňce Nastavení jsme sestavili pass CollectAndCollapse, který slučuje sousední rotace rxx a ryy do jediného XXPlusYYGate. Tento pass jsme již použili při sestavování Trotterových Circuit v kroku 1 (volání pm.run(...)). To jednak snižuje počet dvouqubitových gate, jednak vytváří strukturu, která je vhodnější pro simulaci pomocí tenzorových sítí při pozdějším výpočtu dynamických koeficientů.

2. Snížení na ISA simulátoru. Níže spouštíme přednastavený pass manager Qiskit s optimization_level=3, abychom snížili každý Trotterův Circuit na architekturu instrukční sady (ISA) simulátoru.

aer_sim = AerSimulator()
pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)

isa_circs_all_times = [
    pm_sim.run([deepcopy(c) for c in mpf_circuits])
    for mpf_circuits in all_circs
]

Krok 3: Spuštění pomocí Qiskit primitiv

Pro příklad v malém měřítku spouštíme ISA-snížené Trotterovy Circuity pomocí primitiva EstimatorV2 s backendem Aer. Tím získáme nezašuměnou referenční hodnotu pro každou dvojici $(k_j, t)$ — jsou to hodnoty $\langle A \rangle_{k_j}(t)$, které MPF v kroku 4 zkombinuje. Procházíme rozsah časů evoluce, abychom mohli později vykreslit úplnou časovou křivku každého jednotlivého součinového vzorce i MPF.

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

Krok 4: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu

Krok 4 je místem, kde je MPF skutečně sestaven. I když jsou koeficienty $x_j$ zde vypočítány (a u dynamické varianty může být tento výpočet náročný), koncepčně jde o klasický recept pro kombinování kvantových měření z kroku 3 do jediné opravené očekávané hodnoty — proto celý pracovní postup výpočtu koeficientů a kombinace považujeme za post-processing.

Abychom posoudili, jak dobře MPF sleduje skutečnou dynamiku, nejprve vypočítáme přesné časově vyvinuté očekávané hodnoty přímou exponenciací Hamiltoniánu. To je proveditelné pouze díky tomu, že $L = 10$; v příkladu ve velkém měřítku na hardwaru níže se místo toho budeme muset spolehnout na odhady pomocí tenzorových sítí.

exact_expvals = []
for t in exact_evolution_times:
    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)

Statické MPF koeficienty

Statické MPF používají koeficienty $x_j$, které jsou nezávislé na čase evoluce, Hamiltoniánu a počátečním stavu. Sestavíme lineární soustavu $Ax = b$ popsanou v sekci Kontext a vyřešíme ji pro koeficienty. Matice $A$ je určena počty Trotterových kroků $k_j$, řádem $\chi$ součinového vzorce a tím, zda je vzorec symetrický (což řídí exponenty $\eta_n$).

Pro náš příklad v malém měřítku používáme $k_j = [1, 2, 4]$ s nesymetrickým Suzuki-Trotterovým vzorcem řádu $2\chi=2$ (tedy $\chi=1$ a $\eta_n = 2 + n$, což dává $\eta_0 = 2,\, \eta_1 = 3$). Soustava má tvar:

$$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}, \quad b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}.$$

První řádek vynucuje nezkreslenost ($\sum_j x_j = 1$); druhý a třetí řádek ruší vedoucí člen $1/k^2$ a člen $1/k^3$ dalšího řádu Trotterovy chyby.

Sestavení LSE

Používáme setup_static_lse z qiskit_addon_mpf.static k sestavení matice $A$ a vektoru pravé strany $b$ popsaných výše. Matice $A$ nezávisí pouze na $k_j$, ale i na naší volbě součinového vzorce — konkrétně na jeho řádu $\chi$ a na tom, zda je symetrický. Příznak symmetric řídí vzor exponentů $\eta_n$ (symetrické vzorce produkují pouze členy Trotterovy chyby se sudými mocninami; viz odkaz [1]). Všimni si, že jak ukazuje odkaz [2], nastavení symmetric=True není striktně nutné, ani když je základní PF symetrický — nesymetrická LSE zůstává platná (vynucuje nadbytečná, nepotřebná omezení).

Pro náš příklad jsme již v kroku 1 nastavili order = 2 a symmetric = False.

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

Prohlédni si sestrojenou matici $A$ a vektor $b$, abys potvrdil, že odpovídají soustavě uvedené výše.

lse.A

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

lse.b

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

S hotovou LSE vyřešíme statické koeficienty $x_j$ pomocí lse.solve() (jedná se o přímé řešení $x = A^{-1}b$).

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]

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

Alternativně k výpočtu $x = A^{-1}b$ lze použít setup_exact_model pro sestavení instance cvxpy.Problem, která používá LSE jako omezení a jejíž optimální řešení dá $x$.

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

[ 0.04761905 -0.57142857  1.52380952]

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

L1 norm of the exact coefficients: 2.1428571428556378

Optimalizace pro $x$ pomocí přibližného modelu

Může se stát, že $L_1$-norma pro zvolenou sadu hodnot $k_j$ je považována za příliš vysokou. Pokud tomu tak je a nemůžeš zvolit jinou sadu hodnot $k_j$, můžeš použít přibližné řešení, které omezí $L_1$-normu na zvolený práh při současné minimalizaci $\|Ax - b\|$. Podívej se na průvodce Jak použít přibližný model.

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

Dynamické MPF koeficienty

Statický MPF ruší členy Trotterovy chyby způsobem nezávislým na Hamiltoniánu a stavu, takže nemusí nutně produkovat nejmenší možnou chybu aproximace pro daný Hamiltonián a počáteční stav. Dynamický MPF (odkazy [2], [3]) místo toho nalézá časově závislé koeficienty $x_i(t)$, které minimalizují vzdálenost ve Frobeniově normě $\|\rho(t) - \mu^D(t)\|_F^2$ v každém čase $t$. Jak bylo ukázáno v Kontextu, to vyžaduje matici překryvů $M_{ij}(t)$ mezi Trotterem vyvinutými stavy a překryv $L_i(t)$ s přesným stavem — obě tyto veličiny odhadujeme pomocí backendů tenzorových sítí (TeNPy) v qiskit_addon_mpf.

K sestavení dynamické LSE potřebujeme tři složky:

1. Přibližnou továrnu na evolvery, kterou addon spustí pro každé $k_j$, aby vytvořil $\rho_{k_j}(t)$ jako MPS/MPO. Sestavíme ji z vrstvené struktury Trotterova Circuit druhého řádu (jedna vrstva na slice_by_depth), zabalené jako LayerwiseEvolver s parametry ořezávání TeNPy.

2. Přesnou továrnu na evolvery, která produkuje vysoce přesnou referenci $\rho(t)$. Jako náhradu za přesnou evoluci používáme Suzuki-Trotterův Circuit čtvrtého řádu s malým časovým krokem (dt=0.1, order=4).

3. Továrnu na identitu a MPS počátečního stavu, které inicializují simulaci TeNPy.

Buňka níže sestrojuje přibližnou továrnu na evolvery.

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

varování

Možnosti LayerwiseEvolver, které určují podrobnosti simulace tenzorové sítě, je třeba volit pečlivě, aby se předešlo nastavení špatně definovaného optimalizačního problému.

Přesně časově vyvinutý stav aproximujeme pomocí Suzuki-Trotterova vzorce čtvrtého řádu s malým časovým krokem dt=0.1. Parametry ořezávání (truncation) TeNPy mohou ovlivnit přesnost, proto je důležité prozkoumat rozsah hodnot.

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

Nakonec definujeme identity_factory, který vytváří počáteční stav MPO, a připravíme počáteční Néelův stav jako MPS odpovídající mřížce použité ve vrstveném Trotterově modelu.

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

mps_initial_state = MPS_neel_state(models[0].lat)

S hotovými továrnami nyní vypočítáme dynamické koeficienty v každém čase evoluce. Pro každé $t$ setup_dynamic_lse sestaví příslušné matice překryvů pomocí TeNPy a setup_frobenius_problem vrátí cvxpy.Problem, který minimalizuje náklady ve Frobeniově normě. Řešič vrátí koeficienty $x_j(t)$ přizpůsobené danému času; shromažďujeme je v mpf_dynamic_coeffs_list. Pokud řešič pro dané $t$ selže, vrátíme se k nulovým koeficientům, aby smyčka mohla pokračovat.

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

Kombinace Trotterových očekávaných hodnot s koeficienty MPF

Nyní vyhodnotíme $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ pro každou sadu koeficientů (statická-přesná, statická-přibližná a dynamická), propagujeme standardní chyby jednotlivých Circuit a vykreslíme výslednou časovou řadu proti křivce přesné diagonalizace.

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

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

# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

plt.title(f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ vs time")
plt.xlabel("Time")
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()

Graf výše ilustruje vzájemné působení Trotterovy chyby a chyby vzorkování.

• Trotterova chyba. Jednotlivé součinové vzorce (šedé značky) se s rostoucím časem stále více odchylují od přesné křivky. Circuit s $k=1$ má největší odchylku a je nejmělčí, ale zároveň se už nachází v režimu, kde $t/k \gtrsim 1$, takže vedoucí člen chyby $1/k^{2}$ je velký. Kombinace MPF (barevné značky) ruší několik z těchto vedoucích členů Trotterovy chyby, takže sledují přesnou křivku mnohem těsněji než kterýkoli jednotlivý Circuit s daným $k_j$. Zbývající rozdíl odráží členy Trotterovy chyby vyššího řádu, které MPF neruší: statický MPF řádu $2$ s $r=3$ ruší pouze první dva řády chyby a při velkém $t/k_{\min}$ nakonec dominuje nezrušený zbytek — MPF tedy nezaručuje, že velmi mělké Circuity zůstanou přesné v libovolných časech.

• Chyba vzorkování. Širší chybové úsečky u křivek MPF jsou přímým důsledkem lineární kombinace: propagace nezávislých standardních chyb jednotlivých Circuit $\sigma_{k_j}$ dává celkový rozptyl $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Čím větší je tedy $\|x\|_2$ (a v praxi $\|x\|_1$, kterou kontrolujeme), tím více snímků je potřeba k dosažení dané cílové nejistoty. To je kompromis, na kterém je založena možnost přibližného řešiče v Kontextu: omezujeme $\|x\|_1$, abychom udrželi tuto režii zvladatelnou. Na rozdíl od Trotterovy chyby se chyba vzorkování zásadně zmenšuje s $1/\sqrt{N_{\text{shots}}}$, takže ji lze vždy snížit spotřebováním více snímků.

V příkladu ve velkém měřítku na hardwaru níže vstupuje hardwarový šum jako další zdroj chyby u každé hodnoty $\langle A \rangle_{k_j}$, což je podobně zesíleno koeficienty MPF. V dané sekci uvidíme, jak potlačení chyb interaguje s MPF.

Příklad ve velkém měřítku na hardwaru

V této sekci rozšíříme úlohu nad rámec toho, co lze přesně simulovat. Reprodukujeme některé výsledky uvedené v odkazu [3], a to pomocí 50qubitového XXZ řetězce v čase $t = 3$. Postupujeme podle stejného čtyřkrokového pracovního postupu jako u příkladu v malém měřítku, tentokrát cílíme na skutečný kvantový hardware s potlačením chyb. Stejně jako v šabloně je každý krok označen přímo v kódu a jeden krok může zahrnovat více buněk, pokud stojí za to prozkoumat jeho mezivýsledky. Mapování odpovídá příkladu v malém měřítku: definice Hamiltoniánu, volba parametrů Trottera, výpočet koeficientů MPF (statických i dynamických) a sestavení Circuit. Klíčové rozdíly jsou:

• XXZ Hamiltonián na 50 místech s náhodnými vazbami taženými z $\mathcal{U}(0.5, 1.5)$ (odkaz [3]).

• Symetrický Trotterův vzorec druhého řádu s $k_j = [3, 4, 6]$ (tedy $\chi=1$, symmetric=True).

• Jediný pevný čas evoluce $t = 3$. S $k_{\min}=3$ to dává $t/k_{\min}=1$, čímž se mělké složky udržují v Trotterově konvergenčním režimu, kde je model vedoucí chyby, na kterém MPF spoléhá, platný.

• Dodatečné srovnávací spuštění jediného Circuit s $k = 10$ Trotterovými kroky, použité jako základní hodnota (baseline). Zvolili jsme $k = 10$, protože jeho dvouqubitová hloubka na hardwaru je hlubší než nejhlubší složka MPF ($k_{\max}=6$) plus režie spuštění více Circuit MPF — dostatečně hluboká na to, aby byla omezena šumem, což je právě režim, v němž se očekává, že kombinace MPF překoná základní hodnotu jediného Circuit. Jde o srovnání typu „jeden hluboký Circuit“ proti kombinaci MPF, nikoli o Circuit cílící na efektivní Trotterovu chybu MPF (což by vyžadovalo mnohem více kroků).

Všimni si, že i když se stále nacházíme v kroku 1 (mapování a konstrukce Circuit), v této buňce také předpočítáme dynamické koeficienty spolu se statickými. Dynamické koeficienty závisí na $H$ a $t$, ale ne na kvantových měřeních, takže je lze vypočítat kdykoli před krokem 4. Děláme to nyní, abychom udrželi veškeré nastavení specifické pro MPF na jednom místě.

# -------------------------Step 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# XXZ Hamiltonian with random couplings (Ref. [3])
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,
    )

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

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

# Static coefficients
lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(f"Static coefficients: {mpf_coeffs}")
print(f"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}")

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(f"Approximate coefficients: {coeffs_approx.value}")
print(f"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}")

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

layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

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

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

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

mps_initial_state = MPS_neel_state(models[0].lat)

print(f"Computing dynamic coefficients for time={total_time}")
lse_dyn = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs_dyn.value
except Exception as error:
    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))
    print(error, "Calculation Failed")

# -------------------------Step 1 (cont): Build circuits-------------------------
mpf_circuits = []
for k in mpf_trotter_steps:
    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)

# Baseline "single deep circuit" comparison run with k=10 Trotter steps.
# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus
# the overhead of running multiple circuits, pushing it into the noise-limited
# regime where MPF is expected to outperform. It does NOT target the MPF's effective
# Trotter error (which would require many more steps).
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=10, order=order),
    time=total_time,
)
comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)
mpf_circuits.append(comp_circuit)

Static coefficients: [ 0.42857143 -1.82857143  2.4       ]
L1 norm: 4.65714285714286
Approximate coefficients: [-0.4942491   0.40206845  1.09218065]
L1 norm (approx): 1.9884981979026675
Computing dynamic coefficients for time=3

Nyní optimalizujeme Circuity pro zvolený Backend. Používáme přednastavený pass manager Qiskit s optimization_level=3, který automaticky vybere vhodnou sadu fyzických qubitů a namapuje každý Circuit na topologii zařízení.

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

<IBMBackend('ibm_fez')>

Spouštění hlubších Circuit na skutečném hardwaru vyžaduje agresivní potlačení chyb. Povolujeme dynamické oddělování (dynamical decoupling), twirling gate a měření, potlačení chyb měření a extrapolaci k nulovému šumu (ZNE). Všimni si, že faktory šumu ZNE, které zde používáme (1, 1.2, 1.4), jsou menší než ve scénáři s mělkým Circuit, protože hlubší složky MPF jsou už blízko prahu šumu a velká zesílení šumu by je posunula za bod, kde je extrapolace ZNE spolehlivá.

Odesíláme všechny čtyři Circuity (tři složky MPF při $k_j = [3, 4, 6]$ plus základní hodnotu $k = 10$) v jediné úloze Estimator.

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Error suppression/mitigation
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 = ["TUT_MPF"]

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

Z výsledku úlohy získáme očekávané hodnoty a směrodatné odchylky jednotlivých Circuit, poté je zkombinujeme s každou sadou koeficientů MPF přesně stejně jako v příkladu v malém měřítku: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$, s propagovaným rozptylem $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$.

# -------------------------Step 4-------------------------
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.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]
[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]

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.5665938395816946 +- 0.3925273058119915
Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166
Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169

sym = {3: "^", 4: "s", 6: "p"}
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="10 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"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ 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()

Několik postřehů k výsledkům na hardwaru výše:

• Jít hlouběji na hardwaru nic nestojí zdarma. Základní hodnoty jednotlivých Circuit to říkají přímo: Circuit $k = 6$ je v podstatě přesný ($-0,256$ oproti referenci $-0,244$), zatímco hlubší základní hodnota $k = 10$ je horší ($-0,061$, odchylka přibližně $0,18$), nikoli lepší. Jakmile je Trotterova chyba už malá, přidávání kroků hlavně prohlubuje Circuit a hromadí více šumu gate a dekoherence. Přesně toto je režim, pro který jsou MPF navrženy: dosáhnout přesnosti hlubokého Circuit pouze pomocí mělkých složek.

• Malonormová MPF poráží hluboký jednotlivý Circuit. Přibližná statická MPF (omezená na $\|x\|_1 \approx 2$) dosahuje hodnoty $-0,259$, v rozmezí přibližně $0,015$ od reference a mnohem blíže než základní hodnota $k = 10$. Dynamická MPF ($-0,127$) také pohodlně poráží tuto základní hodnotu. Obě kombinují pouze mělké Circuity $k_j = [3, 4, 6]$, a přesto obnovují odpověď, kterou hluboký jednotlivý Circuit nedokázal.

• Norma koeficientů je důležitější než matematická optimalita. Přesná statická MPF má $\|x\|_1 = 4,66$ a je nejhorším odhadem ze všech ($-0,567$, odchylka o více než $0,3$): velká norma koeficientů zesiluje zbytkový šum gate, dekoherenci a chybu ZNE u každé hodnoty $\langle A \rangle_{k_j}$ přibližně stejným faktorem, čímž přehlušuje zrušení Trotterovy chyby, které tím získává. Omezení normy (přibližný statický řešič, $\|x\|_1 \approx 2$) toto přehlušení odstraňuje a dává nejlepší odhad — i když jeho koeficienty už přesně neruší vedoucí Trotterovu chybu.

• Jednotlivé mělké Circuity mohou být stále konkurenceschopné. Samotná složka $k = 6$ ($-0,256$) je zde sama o sobě v podstatě přesná — v tomto běhu je dokonce nepatrně blíže než přibližná statická MPF. Háček je v tom, že předem nevíš, které jednotlivé $k$ leží v ideálním bodě „konvergováno, ale ještě ne omezeno šumem“, a zdánlivě bezpečná volba jednoduše jít hlouběji ($k = 10$), aby se zaručila konvergence Trottera, je právě ta, která selhává. MPF poskytuje principiální kombinaci mělkých Circuit, která nevyžaduje odhadování správné hloubky.

Praktickým poučením je, že na hardwaru by MPF měly být párovány se silným potlačením chyb u každé jednotlivé hodnoty $\langle A \rangle_{k_j}$, $L_1$-norma koeficientů by měla být udržována mírná (použij přibližný řešič nebo dynamický MPF) a Trotterovy kroky $k_j$ by měly být zvoleny tak, aby $t/k_{\min} \lesssim 1$ — zde $k_{\min} = 3$ při $t = 3$ dává $t/k_{\min} = 1$, čímž se složky udržují v konvergentním režimu, kde je model vedoucí chyby, na kterém statická MPF spoléhá, platný. Při těchto volbách se malonormové MPF zde vyrovnají konvergovanému jednotlivému Circuit, zatímco naivní základní hodnota „prostě jdi hlouběji“ nikoli, čímž se obnovuje výhoda hloubky versus přesnosti uvedená v odkazu [3]. Všimni si také, že jednotlivé běhy jsou zašuměné — při jiném odeslání stejné úlohy (nebo na jiném Backend) se přesné pořadí může posunout; robustní trendy jsou, že MPF s malým $\|x\|_1$ fungují dobře, přesná statická MPF s velkým $\|x\|_1$ je zesílena hardwarovým šumem a nadměrně hluboký jednotlivý Circuit je omezen šumem.

Další kroky

Doporučení

Pokud tě tato práce zaujala, mohl by tě zajímat následující materiál:

• Jak zvolit Trotterovy kroky pro MPF — praktický návod, jak zvolit hodnoty $k_j$ tak, aby se předešlo nestabilitám
• Jak použít přibližný model — ladění omezení $L_1$-normy a možností řešiče pro přibližnou statickou MPF
• Referenční dokumentace API qiskit-addon-mpf — úplná dokumentace pro statické, dynamické a backendové moduly

Reference

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

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

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