Formule multi-prodotto per ridurre l'errore di Trotter

Utilizzo stimato: Quattro minuti su un processore Heron r2 (NOTA: Questa è solo una stima. Il tempo di esecuzione potrebbe variare.)

Risultati di apprendimento

• Come le formule multi-prodotto (MPF) riducano l'errore di Trotter nella simulazione Hamiltoniana combinando valori di aspettazione da più circuiti poco profondi

• Quando le MPF sono vantaggiose rispetto alle formule prodotto standard e quando non sono lo strumento giusto

• Come calcolare i coefficienti MPF statici e dinamici usando il pacchetto qiskit_addon_mpf

• Come eseguire un flusso di lavoro MPF end-to-end su hardware IBM Quantum®, inclusa la transpilazione, la mitigazione degli errori e il post-processing

Prerequisiti

• Metodi di compilazione per circuiti di simulazione Hamiltoniana — introduce i circuiti Trotter (formula prodotto) in Qiskit.

• Formule prodotto in Qiskit, in particolare le classi di sintesi SuzukiTrotter e LieTrotter.

• Primitive Qiskit e l'interfaccia Estimator.

Contesto

Cosa sono le formule multi-prodotto?

Quando si simulano sistemi quantistici su un computer quantistico, un compito centrale è approssimare l'operatore di evoluzione temporale $e^{-iHt}$ per un'Hamiltoniana $H$. L'approccio standard usa formule prodotto (PF), note anche come decomposizioni di Trotter-Suzuki. Queste decompongono $H = \sum_{a=1}^d F_a$ in termini le cui unitarie individuali $e^{-iF_a t}$ sono efficienti da implementare, e poi approssimano l'evoluzione completa come un prodotto ordinato di queste unitarie più semplici.

La formula prodotto del primo ordine (Lie-Trotter) è:

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

che comporta un errore quadratico: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. Formule simmetriche di ordine superiore $S_{2\chi}(t)$, dove $\chi$ indica l'ordine della formula prodotto simmetrica (vedi Rif. [1]), convergono più rapidamente come $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, ma al costo di circuiti più profondi per passo.

Per ridurre l'errore a un ordine $\chi$ fisso, di solito si divide il tempo totale di evoluzione $t$ in $k$ passi di Trotter più piccoli. Ogni passo approssima $e^{-iHt/k}$ tramite una formula prodotto e i passi vengono concatenati:

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

Per una formula simmetrica di ordine $2\chi$, l'errore residuo di Trotter scala quindi come $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Quindi aumentare $k$ sopprime rapidamente l'errore di Trotter — ma approfondisce anche linearmente il circuito, e su hardware rumoroso ciò significa più rumore delle porte accumulato. Questa tensione tra errore di Trotter (favorisce $k$ maggiore) e rumore hardware (favorisce $k$ minore) è esattamente ciò che le formule multi-prodotto sono progettate per risolvere. Nota che le MPF riguardano la combinazione di risultati da scelte diverse di $k$ a un ordine $\chi$ fisso — non cambiano l'ordine della formula prodotto sottostante.

Le formule multi-prodotto (MPF) [1] costruiscono una combinazione lineare pesata di valori di aspettazione ottenuti da diversi circuiti di Trotter meno profondi, ciascuno usando un numero diverso di passi di Trotter $k_1, k_2, \ldots, k_r$ (un insieme di $r$ conteggi di passi):

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

dove $\langle A \rangle_{k_j}(t)$ è il valore di aspettazione di un osservabile $A$ al tempo $t$ stimato da un circuito di Trotter con $k_j$ passi, e i coefficienti $\{x_j\}_{j=1}^r$ sono scelti in modo che i termini di errore di Trotter dominanti nella combinazione si cancellino. Rivisiteremo questa espressione nel Passo 4, dove la valutiamo esplicitamente per combinare i nostri risultati di Trotter. Il punto pratico chiave è che il circuito più profondo nella MPF richiede solo $k_{\max}$ passi, molto meno del singolo $k$ che sarebbe necessario per raggiungere direttamente lo stesso errore di Trotter effettivo. I circuiti meno profondi rendono l'approccio MPF più adatto all'hardware rumoroso.

Come vengono determinati i coefficienti?

Esistono due famiglie di coefficienti MPF:

I coefficienti statici sono indipendenti dall'Hamiltoniana, dallo stato iniziale e dal tempo di evoluzione. Sono trovati risolvendo un sistema lineare $Ax = b$ che impone la cancellazione dei termini di errore di Trotter dominanti. Per un insieme di passi di Trotter $\{k_j\}_{j=1}^r$ usati con una formula prodotto simmetrica di ordine $2\chi$, espandere l'errore di Trotter in potenze inverse di $k_j$ porta a equazioni di vincolo della forma:

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

dove gli esponenti interi $\{\eta_n\}$ sono gli ordini dei successivi termini di errore di Trotter per la formula prodotto scelta. Per una PF simmetrica di ordine $2\chi$, l'errore principale in $\left[S_{2\chi}(t/k)\right]^k$ scala come $1/k^{2\chi}$, con correzioni successive a $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — quindi gli esponenti sono $\eta_n = 2\chi + 2n$. Per PF non simmetriche, contribuiscono sia le potenze dispari che quelle pari e $\eta_n = 2\chi + n$. Vedi il Rif. [1] per la derivazione completa. La prima equazione nel sistema sopra garantisce l'assenza di bias (la MPF riproduce il valore di aspettazione esatto nel limite $k_j \to \infty$), e le restanti $r-1$ equazioni cancellano successivamente i primi $r-1$ termini di errore di Trotter. Quando la norma $L_1$ risultante $\|x\|_1$ è troppo grande (il che amplifica il rumore di campionamento), puoi invece risolvere un'ottimizzazione approssimata che limita $\|x\|_1$ minimizzando $\|Ax - b\|$.

I coefficienti dinamici [2], [3] dipendono inoltre dall'Hamiltoniana, dallo stato iniziale e dal tempo di evoluzione $t$. Essi minimizzano la distanza in norma di Frobenius tra il vero stato evoluto nel tempo e l'approssimazione MPF:

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

dove $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ è la matrice di Gram delle sovrapposizioni tra stati evoluti da Trotter per diversi conteggi di passi $k_i, k_j$, e $L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ misura la sovrapposizione con lo stato esatto (approssimato). In questo tutorial queste quantità vengono calcolate in modo efficiente utilizzando metodi di tensor-network, in particolare i backend basati su TeNPy in qiskit_addon_mpf.

Quando usare le MPF

Le MPF sono più vantaggiose quando:

• La profondità del circuito è il collo di bottiglia. Se il rumore hardware limita la profondità che puoi eseguire, usa le MPF per ottenere una maggiore accuratezza di Trotter effettiva partendo da circuiti più superficiali.

• Hai bisogno di valori di aspettazione accurati, non della preparazione dello stato completo. Le MPF operano a livello dei valori di aspettazione — combinano numeri classici, non stati quantistici. Sono quindi ideali per la stima di osservabili quando si utilizza la primitiva Estimator.

• Combini un numero modesto di conteggi di passi di Trotter. Tipicamente combinare $r = 3$–$5$ diversi conteggi di passi $k_j$ è sufficiente per cancellare diversi termini di errore di Trotter principali mantenendo $\|x\|_1$ gestibile.

Quando le MPF potrebbero non aiutare

• Tempi di evoluzione molto brevi. Quando $t$ è abbastanza piccolo che una singola formula di Trotter di ordine basso è già accurata, l'overhead di eseguire più circuiti non è necessario.

• Compiti di preparazione dello stato. Le MPF producono un valore di aspettazione corretto, non uno stato quantistico corretto. Se hai bisogno dello stato effettivo evoluto nel tempo (ad esempio, come input per un'altra subroutine quantistica), le MPF non si applicano.

• Conteggi di passi di Trotter che violano il regime di convergenza. La derivazione a coefficienti statici espande ogni singolo $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ come serie in $t/k_j$; questa espansione converge bene solo quando $t/k_{\min} \lesssim 1$. Se $k_{\min}$ è scelto troppo piccolo per il dato $t$, il circuito più superficiale è ben al di fuori del regime perturbativo, i termini di errore di ordine superiore che l'MPF lascia non cancellati diventano grandi, e la cancellazione può richiedere coefficienti elevati. La norma $L_1$ $\|x\|_1$ è il diagnostico pratico: quando $\|x\|_1 \gg 1$, l'overhead di campionamento $\propto \|x\|_1^2$ potrebbe superare la riduzione dell'errore di Trotter. Vedi la guida sulla scelta dei passi di Trotter per i dettagli.

Cosa copre questo tutorial

Questo tutorial illustra un workflow MPF end-to-end in due fasi. Prima, un esempio con simulatore su piccola scala (catena di Heisenberg a 10 qubit) dimostra come impostare il problema, calcolare i coefficienti MPF statici e dinamici, e confrontare i valori di aspettazione risultanti con la diagonalizzazione esatta. Poi, un esempio hardware su larga scala (catena XXZ a 50 qubit) mostra come effettuare il transpiling, eseguire su hardware IBM Quantum con mitigazione degli errori, e post-elaborare i risultati utilizzando i coefficienti MPF. Per tutto il tutorial, utilizziamo il pacchetto qiskit_addon_mpf insieme agli strumenti Qiskit standard.

Requisiti

Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:

• Qiskit SDK v2.0 o successivo, con supporto per la visualizzazione

• Qiskit Runtime v0.22 o successivo (pip install qiskit-ibm-runtime)

• Simulatore Qiskit Aer (pip install qiskit-aer)

• Addon MPF Qiskit con il backend TeNPy (pip install "qiskit-addon-mpf[tenpy]")

• Utilità di componenti aggiuntivi Qiskit (pip install qiskit-addon-utils)

• SciPy (pip install scipy)

Configurazione

Di seguito raccogliamo tutti gli import dei pacchetti utilizzati in questo tutorial in un'unica cella. Definiamo inoltre un pass del transpiler CollectAndCollapse che fonde le rotazioni rxx e ryy adiacenti in un unico XXPlusYYGate. Questo pass viene applicato sia durante la costruzione del circuito nel Passo 1 (per mantenere basso il conteggio dei gate) sia indirettamente quando estraiamo la struttura a strati per la MPF dinamica nel Passo 4 (TeNPy si aspetta gate a due qubit, non coppie di rotazioni non fuse).

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_mpf": "qiskit-addon-mpf", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy", "tenpy": "physics-tenpy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

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

Esempio con simulatore su piccola scala

Passo 1: Mappare gli input classici a un problema quantistico

Iniziamo con un modello di Heisenberg a 10 qubit su una linea, utilizzando lo stato di Néel $\vert 0101\ldots01 \rangle$ come stato iniziale. L'Hamiltoniana è:

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

dove $J$ è la forza di accoppiamento tra primi vicini. Misuriamo il correlatore ZZ $Z_{L/2-1} Z_{L/2}$ su una coppia di qubit al centro della catena, e utilizziamo passi di Trotter $k_j = [1, 2, 4]$ con una formula prodotto del secondo ordine.

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)

Costruire i circuiti di Trotter

Creiamo i circuiti che implementano le evoluzioni temporali approssimate di Trotter per ogni punto temporale e ogni conteggio di passi di Trotter. Il pass CollectAndCollapse definito nella sezione Configurazione raccoglie le rotazioni XX e YY in singoli gate XX+YY, per prepararsi a una simulazione tensor-network più efficiente in seguito.

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

Step 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Per l'esempio su piccola scala puntiamo al simulatore Aer. Due trasformazioni avvengono prima che i circuiti siano pronti per l'esecuzione:

1. Raccolta dei gate a livello di simulazione dell'Hamiltoniana. Nella cella di Configurazione abbiamo costruito un pass CollectAndCollapse che fonde le rotazioni rxx e ryy adiacenti in un unico XXPlusYYGate. Abbiamo già applicato questo pass quando abbiamo costruito i circuiti di Trotter nel Passo 1 (la chiamata pm.run(...)). Questo riduce sia il conteggio dei gate a due qubit sia produce una struttura più adatta alla simulazione tensor-network per il calcolo dei coefficienti dinamici in seguito.

2. Abbassamento all'ISA del simulatore. Di seguito eseguiamo il pass manager preimpostato di Qiskit con optimization_level=3 per abbassare ogni circuito di Trotter all'architettura del set di istruzioni (ISA) del simulatore.

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
]

Step 3: Eseguire utilizzando le primitive di Qiskit

Per l'esempio su piccola scala eseguiamo i circuiti di Trotter abbassati all'ISA attraverso la primitiva EstimatorV2 supportata da Aer. Facendo così otteniamo un valore di riferimento senza rumore per ogni coppia $(k_j, t)$ — questi sono i valori $\langle A \rangle_{k_j}(t)$ che la MPF combinerà nel Passo 4. Effettuiamo uno sweep sui tempi di evoluzione in modo da poter tracciare successivamente la curva completa della serie temporale di ogni singola formula prodotto e della 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])

Passo 4: Post-elaborare e restituire il risultato nel formato classico desiderato

Il Passo 4 è dove l'MPF viene effettivamente costruito. Anche se i coefficienti $x_j$ vengono calcolati qui (e, per la variante dinamica, questo calcolo può essere intensivo), concettualmente sono una ricetta classica per combinare le misurazioni quantistiche del Passo 3 in un unico valore di aspettazione corretto — quindi trattiamo l'intero workflow di coefficienti e combinazione come post-elaborazione.

Per valutare quanto bene la MPF segue la dinamica reale, prima calcoliamo i valori di aspettazione esatti evoluti nel tempo esponenziando direttamente l'Hamiltoniana. Questo è trattabile solo perché $L = 10$; nell'esempio hardware su larga scala qui sotto dovremo invece affidarci a stime tensor-network.

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)

Coefficienti MPF statici

Le MPF statiche utilizzano coefficienti $x_j$ che sono indipendenti dal tempo di evoluzione, dall'Hamiltoniana e dallo stato iniziale. Impostiamo il sistema lineare $Ax = b$ descritto nella sezione Contesto e risolviamo per i coefficienti. La matrice $A$ è determinata dai conteggi di passi di Trotter $k_j$, dall'ordine $\chi$ della formula prodotto, e da se la formula è simmetrica (il che controlla gli esponenti $\eta_n$).

Per il nostro esempio su piccola scala utilizziamo $k_j = [1, 2, 4]$ con una formula di Suzuki-Trotter non simmetrica di ordine $2\chi=2$ (quindi $\chi=1$ e $\eta_n = 2 + n$, dando $\eta_0 = 2,\, \eta_1 = 3$). Il sistema diventa:

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

La prima riga garantisce l'assenza di bias ($\sum_j x_j = 1$); la seconda e la terza riga cancellano rispettivamente i termini di errore di Trotter principale $1/k^2$ e di ordine successivo $1/k^3$.

Impostare l'LSE

Utilizziamo setup_static_lse da qiskit_addon_mpf.static per assemblare la matrice $A$ e il vettore del lato destro $b$ descritti sopra. La matrice $A$ dipende non solo da $k_j$ ma anche dalla nostra scelta della formula prodotto — in particolare dal suo ordine $\chi$ e da se è simmetrica. Il flag symmetric controlla lo schema di esponenti $\eta_n$ (le formule simmetriche producono solo termini di errore di Trotter con potenze pari; vedi il Rif. [1]). Nota che, come mostrato nel Rif. [2], impostare symmetric=True non è strettamente necessario anche quando la PF sottostante è simmetrica — l'LSE non simmetrico rimane valido (impone vincoli aggiuntivi non necessari).

Per il nostro esempio abbiamo già impostato order = 2 e symmetric = False nel Passo 1.

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

Ispeziona la matrice $A$ e il vettore $b$ costruiti per confermare che corrispondano al sistema scritto sopra.

lse.A

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

lse.b

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

Con l'LSE a disposizione, risolviamo per i coefficienti statici $x_j$ tramite lse.solve() (questa è la soluzione diretta $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]

Ottimizzare per $x$ usando un modello esatto

Alternativamente al calcolo di $x = A^{-1}b$, potete utilizzare setup_exact_model per costruire un'istanza di cvxpy.Problem che utilizza l'LSE come vincoli e la cui soluzione ottimale produrrà $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

Ottimizzare per $x$ usando un modello approssimato

Potrebbe accadere che la norma $L_1$ per l'insieme scelto di valori $k_j$ sia ritenuta troppo alta. Se questo è il caso e non puoi scegliere un diverso insieme di valori $k_j$, puoi utilizzare una soluzione approssimata che limita la norma $L_1$ a una soglia scelta minimizzando $\|Ax - b\|$. Consulta la guida su Come utilizzare il modello approssimato.

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

Coefficienti MPF dinamici

La MPF statica cancella i termini di errore di Trotter in modo agnostico rispetto all'Hamiltoniana e allo stato, quindi non produce necessariamente il minimo errore di approssimazione possibile per una data Hamiltoniana e stato iniziale. La MPF dinamica (Rif. [2], [3]) trova invece coefficienti dipendenti dal tempo $x_i(t)$ che minimizzano la distanza in norma di Frobenius $\|\rho(t) - \mu^D(t)\|_F^2$ a ogni tempo $t$. Come mostrato nella sezione Contesto, questo richiede la matrice di sovrapposizione $M_{ij}(t)$ tra stati evoluti da Trotter e la sovrapposizione $L_i(t)$ con lo stato esatto — entrambe stimate utilizzando backend tensor-network (TeNPy) in qiskit_addon_mpf.

Per impostare l'LSE dinamico abbiamo bisogno di tre ingredienti:

1. Una factory dell'evolutore approssimato che l'addon eseguirà per ogni $k_j$ per produrre $\rho_{k_j}(t)$ come MPS/MPO. La costruiamo dalla struttura a strati del circuito di Trotter di ordine $2$ (uno strato per slice_by_depth), avvolta come un LayerwiseEvolver con parametri di troncamento TeNPy.

2. Una factory dell'evolutore esatto che produce un riferimento $\rho(t)$ ad alta precisione. Utilizziamo un circuito di Suzuki-Trotter del quarto ordine con passo temporale piccolo (dt=0.1, order=4) come proxy dell'evoluzione esatta.

3. Una factory di identità e un MPS dello stato iniziale che avviano la simulazione TeNPy.

La cella sottostante costruisce la factory dell'evolutore approssimato.

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

avviso

Le opzioni di LayerwiseEvolver che determinano i dettagli della simulazione della rete tensoriale devono essere scelte con attenzione per evitare di impostare un problema di ottimizzazione mal definito.

Approssimiamo lo stato esatto evoluto nel tempo con una formula di Suzuki-Trotter del quarto ordine usando un piccolo passo temporale dt=0.1. I parametri di troncamento di TeNPy possono influire sull'accuratezza, quindi è importante esplorare un intervallo di valori.

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

Infine, definiamo una identity_factory che produce lo stato MPO iniziale e prepariamo lo stato iniziale di Néel come un MPS che corrisponde al reticolo utilizzato dal modello di Trotter a strati.

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

mps_initial_state = MPS_neel_state(models[0].lat)

Con le factory pronte, ora calcoliamo i coefficienti dinamici a ogni tempo di evoluzione. Per ogni $t$, setup_dynamic_lse costruisce le matrici di sovrapposizione rilevanti tramite TeNPy, e setup_frobenius_problem restituisce un cvxpy.Problem che minimizza il costo in norma di Frobenius. Il solver restituisce coefficienti $x_j(t)$ adattati a quel tempo; li raccogliamo in mpf_dynamic_coeffs_list. Se il solver fallisce per un dato $t$, ricorriamo a coefficienti nulli in modo che il ciclo continui.

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

Combinare i valori di aspettazione di Trotter con i coefficienti MPF

Ora valutiamo $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ per ogni insieme di coefficienti (statico-esatto, statico-approssimato, e dinamico), propaghiamo gli errori standard per singolo circuito, e tracciamo la serie temporale risultante rispetto alla curva della diagonalizzazione esatta.

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

Il grafico sopra illustra l'interazione tra errore di Trotter ed errore di campionamento.

• Errore di Trotter. Le singole formule prodotto (marcatori grigi) si discostano sempre più dalla curva esatta man mano che il tempo cresce. Il circuito $k=1$ ha la deviazione maggiore ed è il più superficiale, ma è anche già nel regime in cui $t/k \gtrsim 1$, quindi il termine di errore principale $1/k^{2}$ è grande. Le combinazioni MPF (marcatori colorati) cancellano diversi di questi termini di errore di Trotter principali, quindi seguono la curva esatta molto più da vicino di qualsiasi singolo circuito $k_j$. Il divario rimanente riflette i termini di Trotter di ordine superiore che la MPF non cancella: una MPF statica di ordine $2$, $r=3$ elimina solo i primi due ordini di errore, e a $t/k_{\min}$ elevato la coda non cancellata alla fine domina — quindi la MPF non garantisce che i circuiti molto superficiali rimangano accurati a tempi arbitrari.

• Errore di campionamento. Le barre di errore più ampie sulle curve MPF sono una conseguenza diretta della combinazione lineare: propagare errori standard indipendenti per singolo circuito $\sigma_{k_j}$ dà una varianza totale $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Pertanto, più grande è $\|x\|_2$ (e in pratica $\|x\|_1$, che è ciò che controlliamo), più shot sono richiesti per raggiungere una data incertezza obiettivo. Questo è il compromesso alla base dell'opzione del solver approssimato nella sezione Contesto: limitiamo $\|x\|_1$ per mantenere questo overhead gestibile. Fondamentalmente, a differenza dell'errore di Trotter, l'errore di campionamento si riduce con $1/\sqrt{N_{\text{shots}}}$, quindi può sempre essere ridotto spendendo più shot.

Nell'esempio hardware su larga scala qui sotto, il rumore hardware entra come fonte di errore aggiuntiva su ogni $\langle A \rangle_{k_j}$, che viene similmente amplificata dai coefficienti MPF. Vedremo come la mitigazione degli errori interagisce con le MPF in quella sezione.

Esempio hardware su larga scala

In questa sezione scaliamo il problema oltre ciò che è possibile simulare esattamente. Riproduciamo alcuni dei risultati mostrati nel Rif. [3], utilizzando una catena XXZ a 50 qubit al tempo $t = 3$. Seguiamo lo stesso workflow a quattro passi dell'esempio su piccola scala, mirando ora a hardware quantistico reale con mitigazione degli errori. Come nel template, ogni passo è contrassegnato inline nel codice, e un singolo passo può estendersi su più celle quando i suoi output intermedi meritano di essere ispezionati. La mappatura rispecchia l'esempio su piccola scala: definire un'Hamiltoniana, scegliere i parametri di Trotter, calcolare i coefficienti MPF (statici e dinamici), e costruire i circuiti. Le differenze chiave sono:

• Un'Hamiltoniana XXZ su 50 siti con accoppiamenti casuali estratti da $\mathcal{U}(0.5, 1.5)$ (Rif. [3]).

• Una formula di Trotter del secondo ordine simmetrica con $k_j = [3, 4, 6]$ (quindi $\chi=1$, symmetric=True).

• Un singolo tempo di evoluzione fisso $t = 3$. Con $k_{\min}=3$ questo dà $t/k_{\min}=1$, mantenendo i costituenti superficiali all'interno del regime di convergenza di Trotter in cui il modello di errore principale su cui si basa l'MPF è valido.

• Un'ulteriore esecuzione di confronto a singolo circuito con $k = 10$ passi di Trotter, utilizzata come baseline. Abbiamo scelto $k = 10$ perché la sua profondità a due qubit sull'hardware è maggiore del costituente MPF più profondo ($k_{\max}=6$) più l'overhead di eseguire più circuiti MPF — abbastanza profonda da essere limitata dal rumore, che è il regime in cui la combinazione MPF dovrebbe superare la baseline a singolo circuito. È un confronto tra un "singolo circuito profondo" e la combinazione MPF, non un circuito che punta all'errore di Trotter effettivo dell'MPF (che richiederebbe molti più passi).

Nota che anche se siamo ancora nel Passo 1 qui (mappatura e costruzione del circuito), pre-calcoliamo anche i coefficienti dinamici insieme a quelli statici in questa cella. I coefficienti dinamici dipendono da $H$ e $t$ ma non dalle misurazioni quantistiche, quindi possono essere calcolati in qualsiasi momento prima del Passo 4. Lo facciamo ora per mantenere tutta la configurazione specifica per MPF in un unico posto.

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

Ora ottimizziamo i circuiti per il backend scelto. Utilizziamo il pass manager preimpostato di Qiskit con optimization_level=3, che seleziona automaticamente un buon insieme di qubit fisici e instrada ogni circuito sulla topologia del dispositivo.

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

Eseguire circuiti più profondi su hardware reale richiede una mitigazione degli errori aggressiva. Abilitiamo il dynamical decoupling, il twirling di gate e misurazione, la mitigazione degli errori di misurazione, e l'estrapolazione a rumore zero (ZNE). Nota che i fattori di rumore ZNE che utilizziamo qui (1, 1.2, 1.4) sono più piccoli rispetto a uno scenario con circuito superficiale, poiché i costituenti MPF più profondi sono già vicini alla soglia di rumore e ampie amplificazioni di rumore li spingerebbero oltre il punto in cui l'estrapolazione ZNE è affidabile.

Inviamo tutti e quattro i circuiti (tre costituenti MPF a $k_j = [3, 4, 6]$ più la baseline $k = 10$) in un unico job 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)
    ]
)

Estraiamo i valori di aspettazione e le deviazioni standard per singolo circuito dal risultato del job, quindi li combiniamo con ciascun insieme di coefficienti MPF esattamente come nell'esempio su piccola scala: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$, con varianza propagata $\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()

Alcune osservazioni sui risultati hardware sopra riportati:

• Andare più in profondità non è gratuito sull'hardware. Le baseline a singolo circuito raccontano la storia direttamente: il circuito $k = 6$ è essenzialmente esatto ($-0.256$ rispetto al riferimento $-0.244$), eppure la baseline più profonda $k = 10$ è peggiore ($-0.061$, con uno scarto di $\sim 0.18$), non migliore. Una volta che l'errore di Trotter è già piccolo, aggiungere passi principalmente approfondisce il circuito e accumula più rumore di gate e decoerenza. Questo è precisamente il regime per cui le MPF sono costruite: raggiungere l'accuratezza di un circuito profondo utilizzando solo costituenti superficiali.

• Un MPF a norma piccola batte il singolo circuito profondo. L'MPF statico approssimato (limitato a $\|x\|_1 \approx 2$) si attesta a $-0.259$, entro $\sim 0.015$ dal riferimento e molto più vicino della baseline $k = 10$. L'MPF dinamico ($-0.127$) batte anch'esso comodamente quella baseline. Entrambi combinano solo i circuiti superficiali $k_j = [3, 4, 6]$, eppure recuperano una risposta che il singolo circuito profondo non è riuscito a ottenere.

• La norma dei coefficienti conta più dell'ottimalità matematica. La MPF statica esatta ha $\|x\|_1 = 4.66$ ed è il peggiore stimatore di tutti ($-0.567$, con uno scarto superiore a $0.3$): la grande norma dei coefficienti amplifica il rumore di gate residuo, la decoerenza e l'errore ZNE su ciascun $\langle A \rangle_{k_j}$ approssimativamente dello stesso fattore, sopraffacendo la cancellazione dell'errore di Trotter che ottiene. Limitare la norma (il solver statico approssimato, $\|x\|_1 \approx 2$) elimina questo effetto e fornisce la stima migliore — anche se i suoi coefficienti non cancellano più esattamente l'errore di Trotter principale.

• I singoli circuiti superficiali possono ancora essere competitivi. Il solo costituente $k = 6$ ($-0.256$) è di per sé essenzialmente esatto qui — in questa esecuzione è persino marginalmente più vicino della MPF statica approssimata. Il problema è che non sai in anticipo quale singolo $k$ si trova nel punto ottimale di "convergente ma non ancora limitato dal rumore", e la scelta apparentemente sicura di semplicemente andare più in profondità ($k = 10$) per garantire la convergenza di Trotter è esattamente quella che fallisce. La MPF fornisce una combinazione motivata di circuiti superficiali che non richiede di indovinare la profondità giusta.

L'insegnamento pratico è che sull'hardware, le MPF dovrebbero essere abbinate a una forte mitigazione degli errori su ciascun singolo $\langle A \rangle_{k_j}$, la norma $L_1$ dei coefficienti dovrebbe essere mantenuta modesta (usa il solver approssimato, o la MPF dinamica), e i passi di Trotter $k_j$ dovrebbero essere scelti in modo che $t/k_{\min} \lesssim 1$ — qui $k_{\min} = 3$ a $t = 3$ dà $t/k_{\min} = 1$, mantenendo i costituenti all'interno del regime convergente in cui il modello di errore principale su cui si basa la MPF statica è valido. Con queste scelte, le MPF a norma piccola qui eguagliano un singolo circuito convergente mentre la baseline ingenua "vai semplicemente più in profondità" non lo fa, recuperando il vantaggio profondità-contro-accuratezza mostrato nel Rif. [3]. Nota anche che le esecuzioni individuali sono rumorose — in un diverso invio dello stesso job (o un backend diverso), l'ordinamento esatto può cambiare; le tendenze robuste sono che le MPF a $\|x\|_1$ piccolo funzionano bene, la MPF statica esatta a $\|x\|_1$ grande è amplificata dal rumore hardware, e il singolo circuito eccessivamente profondo è limitato dal rumore.

Prossimi passi

Raccomandazioni

Se hai trovato interessante questo lavoro, potresti essere interessato al seguente materiale:

• Come scegliere i passi di Trotter per un MPF — indicazioni pratiche sulla selezione dei valori $k_j$ per evitare instabilità
• Come utilizzare il modello approssimato — regolazione del vincolo sulla norma $L_1$ e delle opzioni del solver per l'MPF statico approssimato
• Riferimento API qiskit-addon-mpf — documentazione completa per i moduli statici, dinamici e backend

References

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