Formule multiprodotto per ridurre l'errore di Trotter

Stima del tempo di esecuzione: quattro minuti su un processore Heron r2 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)

Risultati di apprendimento

Al termine di questo tutorial, avrai acquisito le seguenti conoscenze:

• In che modo le formule multiprodotto (MPF) riducono l’errore di Trotter nella simulazione hamiltoniana combinando i valori attesi provenienti da più circuiti poco profondi
• Quando le formule MPF sono più vantaggiose rispetto alle formule standard e quando non rappresentano lo strumento più adatto
• Come calcolare i coefficienti MPF statici e dinamici utilizzando il qiskit_addon_mpf pacchetto
• Come eseguire un flusso di lavoro MPF end-to-end su un hardware d IBM Quantum®, comprese la transpilazione, la mitigazione degli errori e la post-elaborazione

Prerequisiti

Consigliamo agli utenti di acquisire familiarità con i seguenti argomenti prima di seguire questo tutorial:

• Metodi di compilazione per circuiti di simulazione hamiltoniana — Introduzione ai circuiti di Trotter (formula del prodotto) in Qiskit.
• Le formule dei prodotti in Qiskit, in particolare le SuzukiTrotter classi di sintesi e LieTrotter .
• Qiskit primitives e l'interfaccia Estimator.

Sfondo

Cosa sono le formule multiprodotto?

Quando si simulano sistemi quantistici su un computer quantistico, un compito fondamentale consiste nell’approssimare l’operatore di evoluzione temporale $e^{-iHt}$ per un hamiltoniano $H$. L’approccio standard utilizza le formule del prodotto (PF), note anche come decomposizioni di Trotter-Suzuki. Questi scompongono $H = \sum_{a=1}^d F_a$ in termini i cui singoli operatori unitari $e^{-iF_a t}$ sono efficienti da implementare, per poi approssimare l'evoluzione completa come un prodotto ordinato di questi operatori unitari più semplici.

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

il che comporta un errore quadratico: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. Le formule simmetriche di ordine superiore $S_{2\chi}(t)$, dove $\chi$ indica l’ordine della formula del prodotto simmetrico (cfr. rif. [1] ), convergono più rapidamente secondo la formula $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, ma a costo di circuiti più complessi per ogni passo.

Per ridurre l'errore a un ordine fisso $\chi$, solitamente si suddivide il tempo totale di evoluzione $t$ in $k$ passi di Trotter più piccoli. Ogni fase approssima un $e^{-iHt/k}$ e mediante una formula di prodotto e le fasi vengono concatenate:

Per una formula simmetrica di ordine $2\chi$, l’errore residuo di Trotter varia proporzionalmente a $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Pertanto, aumentando $k$ si riduce rapidamente l’errore di Trotter, ma si aumenta anche linearmente la profondità del circuito e, su hardware soggetto a rumore, ciò comporta un maggiore accumulo di rumore di gate. Questa tensione tra l'errore di Trotter (che favorisce valori più grandi di $k$ ) e il rumore hardware (che favorisce valori più piccoli di $k$ ) è proprio ciò che le formule multiprodotto sono state progettate per risolvere. Si noti che gli MPF servono a combinare i risultati derivanti da diverse scelte di $k$ in un ordine fisso $\chi$ — non modificano l’ordine della formula del prodotto sottostante.

Le formule multiprodotto (MPF) [1] costruiscono una combinazione lineare ponderata dei valori attesi ottenuti da diversi circuiti di Trotter meno profondi, ciascuno dei quali utilizza un numero diverso di passi di Trotter $k_1, k_2, \ldots, k_r$ (una serie di $r$ step counts):

dove $\langle A \rangle_{k_j}(t)$ è il valore atteso 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 tale che i termini principali dell'errore di Trotter nella combinazione si annullino. Torneremo su questa espressione nel Passo 4, dove la calcoleremo esplicitamente per integrare i nostri risultati di Trotter. Il punto fondamentale dal punto di vista pratico è che il circuito più profondo nell’MPF richiede solo $k_{\max}$ passaggi, un numero di gran lunga inferiore rispetto al singolo $k$ che sarebbe necessario per raggiungere direttamente lo stesso errore effettivo di Trotter. I circuiti meno profondi rendono l'approccio MPF più adatto all'hardware soggetto a rumore.

Come vengono determinati i coefficienti?

Esistono due famiglie di coefficienti MPF:

I coefficienti statici sono indipendenti dall'hamiltoniano, dallo stato iniziale e dal tempo di evoluzione. Si ottengono risolvendo un sistema lineare $Ax = b$ che garantisce l'annullamento dei termini principali dell'errore di Trotter. Per una serie di passi di Trotter $\{k_j\}_{j=1}^r$ utilizzata con una formula del prodotto simmetrico di ordine $2\chi$, lo sviluppo dell'errore di Trotter in potenze inverse di $k_j$ porta a equazioni vincolanti della forma:

dove gli esponenti interi $\{\eta_n\}$ sono gli ordini dei termini successivi dell'errore di Trotter per la formula di prodotto scelta. Per un PF simmetrico di ordine $2\chi$, l’errore principale in $\left[S_{2\chi}(t/k)\right]^k$ varia proporzionalmente a $1/k^{2\chi}$, con correzioni successive in $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — quindi gli esponenti sono $\eta_n = 2\chi + 2n$. Per i PF non simmetrici, contribuiscono sia le potenze dispari che quelle pari e $\eta_n = 2\chi + n$. Si veda il rif. [1] per la derivazione completa. La prima equazione del sistema sopra riportato garantisce l’assenza di distorsioni (l’MPF riproduce il valore esatto dell’aspettativa nel limite $k_j \to \infty$ ), mentre le restanti equazioni $r-1$ annullano progressivamente i primi termini di errore di Trotter $r-1$. Quando la norma $L_1$ risultante $\|x\|_1$ è troppo elevata (il che amplifica il rumore di campionamento), è possibile risolvere invece un problema di ottimizzazione approssimativa che limiti $\|x\|_1$ minimizzando al contempo $\|Ax - b\|$.

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

dove $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ è la matrice di Gram delle sovrapposizioni tra gli stati evoluti secondo Trotter per un numero diverso 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 (approssimativo). In questo tutorial tali grandezze vengono calcolate in modo efficiente utilizzando metodi basati su reti di tensori, in particolare i backend TeNPy-based in qiskit_addon_mpf.

Quando utilizzare gli MPF

I fondi MPF offrono i maggiori vantaggi quando:

• La profondità del circuito rappresenta il collo di bottiglia. Se il rumore dell'hardware limita la profondità di esecuzione, utilizzare gli MPF per ottenere una maggiore precisione effettiva di Trotter con circuiti meno profondi.
• Servono valori attesi precisi, non una preparazione completa dello stato. Gli MPF operano a livello dei valori attesi: combinano numeri classici, non stati quantistici. Sono quindi ideali per la stima osservabile quando si utilizza la primitiva Estimator.
• Si combinano un numero modesto di passi di Trotter. In genere, combinando $r = 3$ – $5$, un numero diverso di passaggi $k_j$ è sufficiente per annullare diversi termini di errore di Trotter principali, mantenendo al contempo $\|x\|_1$ a livelli gestibili.

Quando i fondi pensionistici (MPF) potrebbero non essere d'aiuto

• Tempi di evoluzione molto brevi. Quando $t$ è sufficientemente piccolo da rendere già accurata una singola formula di Trotter di ordine basso, non è necessario sostenere il sovraccarico derivante dall'esecuzione di più circuiti.
• Attività di preparazione allo Stato. Gli MPF producono un valore atteso corretto, non uno stato quantistico corretto. Se è necessario lo stato effettivo in funzione del tempo (ad esempio, come input per un’altra subroutine quantistica), gli MPF non sono applicabili.
• Conteggi dei passi al trotto che violano il regime di convergenza. La derivazione del coefficiente statico espande ogni singolo $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ come una serie in $t/k_j$; tale espansione converge bene solo quando $t/k_{\min} \lesssim 1$. Se $k_{\min}$ viene scelto troppo piccolo per il dato $t$, il circuito più superficiale si trova ben al di fuori del regime perturbativo, i termini di errore di ordine superiore che l’MPF lascia non annullati diventano grandi e l’annullamento può richiedere coefficienti elevati. La norma " $L_1$ " $\|x\|_1$ costituisce un criterio diagnostico pratico: quando $\|x\|_1 \gg 1$, il sovraccarico di campionamento $\propto \|x\|_1^2$ potrebbe superare la riduzione dell'errore di Trotter. Per ulteriori dettagli, consulta la guida alla scelta dei gradini Trotter.

Argomenti trattati in questo tutorial

Questo tutorial illustra un flusso di lavoro MPF completo in due fasi. In primo luogo, un esempio di simulazione su piccola scala (catena di Heisenberg a 10 qubit) illustra come impostare il problema, calcolare i coefficienti MPF statici e dinamici e confrontare i valori attesi risultanti con quelli ottenuti tramite diagonalizzazione esatta. Successivamente, un esempio su hardware su larga scala (catena XXZ da 50 qubit) illustra come effettuare la transpilazione, eseguire il codice su un hardwar IBM Quantum e con mitigazione degli errori e post-elaborare i risultati utilizzando i coefficienti MPF. Durante tutto il processo, utilizziamo il qiskit_addon_mpf pacchetto insieme agli strumenti standard di Qiskit.

Requisiti

Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:

• Qiskit SDK v2.0 o versioni successive, con supporto alla visualizzazione
• Qiskit Runtime v0.22 o più tardi (pip install qiskit-ibm-runtime)
• Simulatore Qiskit Aer (pip install qiskit-aer)
• Componente aggiuntivo MPF Qiskit con il backend TeNPy (pip install "qiskit-addon-mpf[tenpy]")
• Utilità aggiuntive di Qiskit (pip install qiskit-addon-utils)
• SciPy (pip install scipy)

Configura

Di seguito riportiamo in un'unica cella tutte le importazioni di pacchetti utilizzate nel corso di questo tutorial. XXPlusYYGateDefiniamo inoltre un CollectAndCollapse passaggio del transpiler che fonde le rotazioni adiacenti rxx e ryy in un'unica rotazione. Questo passaggio viene applicato sia durante la costruzione del circuito nella Fase 1 (per mantenere basso il numero di porte) sia, indirettamente, quando estraiamo la struttura a strati per l’MPF dinamico nella Fase 4 (l’ TeNPy e prevede porte a due qubit, non coppie di rotazioni non fuse).

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 di simulatore su piccola scala

Fase 1: mappare gli input classici su un problema quantistico

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

dove $J$ è l'intensità di accoppiamento tra vicini più prossimi. Misuriamo il correlatore ZZ $Z_{L/2-1} Z_{L/2}$ su una coppia di qubit al centro della catena e utilizziamo i passi di Trotter $k_j = [1, 2, 4]$ con una formula di prodotto di 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)

Output:

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)

Output:

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 circuiti Trotter

Creiamo i circuiti implementando le evoluzioni temporali approssimative di Trotter per ciascun istante e per ciascun numero di passi di Trotter. Il CollectAndCollapse passaggio definito nella sezione “Setup” raggruppa le rotazioni XX e YY in singoli gate XX+YY, al fine di preparare una simulazione più efficiente della rete tensoriale in una fase successiva.

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

Output:

Fase 2: Ottimizzazione del problema per l'esecuzione su hardware quantistico

Per l'esempio su piccola scala, prendiamo come riferimento il simulatore Aer. Prima che i circuiti siano pronti per l'esecuzione, avvengono due trasformazioni:

1. Raccolta dei gate a livello di simulazione hamiltoniana. XXPlusYYGateNella cella “Setup” abbiamo creato un CollectAndCollapse passaggio che unisce le rotazioni adiacenti rxx e ryy in un’unica rotazione. Abbiamo già applicato questa fase quando abbiamo realizzato i circuiti di Trotter nel Passo 1 (la pm.run(...) chiamata). Ciò consente sia di ridurre il numero di gate a due qubit, sia di ottenere una struttura che si presta meglio alla simulazione tramite reti tensoriali per il successivo calcolo dei coefficienti dinamici.

2. Adeguamento all'ISA del simulatore. Di seguito eseguiamo il gestore di passaggi predefinito di Qiskit per optimization_level=3 adattare ciascun 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
]

Passaggio 3: eseguire utilizzando Qiskit primitives

Per l'esempio su piccola scala, eseguiamo i circuiti di Trotter ottimizzati con ISA attraverso la EstimatorV2 primitiva supportata da Aer. In questo modo otteniamo un valore di riferimento privo di rumore per ciascuna coppia $(k_j, t)$ : si tratta dei valori $\langle A \rangle_{k_j}(t)$ che l’MPF combinerà nella Fase 4. Esaminiamo i periodi di evoluzione in modo da poter tracciare in seguito la curva completa della serie temporale di ciascuna formula di prodotto e dell’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])

Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato

Il passaggio 4 è quello in cui viene effettivamente costruito l'MPF. Sebbene i coefficienti $x_j$ vengano calcolati in questa fase (e, per la variante dinamica, tale calcolo possa risultare molto oneroso), concettualmente essi costituiscono una formula classica per combinare le misurazioni quantistiche della Fase 3 in un unico valore atteso corretto; pertanto, consideriamo l'intero flusso di lavoro relativo ai coefficienti e alla combinazione come una fase di post-elaborazione.

Per valutare in che misura l’MPF riesca a riprodurre le dinamiche reali, calcoliamo innanzitutto i valori attesi esatti in funzione del tempo, elevando direttamente all’esponente l’Hamiltoniano. Ciò è fattibile solo perché $L = 10$; nell’esempio di hardware su larga scala riportato di seguito dovremo invece affidarci a stime basate su reti di tensori.

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

Gli MPF statici utilizzano coefficienti $x_j$ che sono indipendenti dal tempo di evoluzione, dall'hamiltoniano e dallo stato iniziale. Si imposta il sistema lineare $Ax = b$ descritto nella sezione "Contesto" e si calcolano i coefficienti. La matrice $A$ è determinata dal numero di passi di Trotter $k_j$, dall'ordine $\chi$ della formula del prodotto e dal fatto che la formula sia simmetrica (il che determina gli esponenti $\eta_n$ ).

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

La prima riga garantisce l’assenza di distorsioni ( $\sum_j x_j = 1$ ); la seconda e la terza riga annullano, rispettivamente, i termini di errore di Trotter di primo ordine $1/k^2$ e di ordine successivo $1/k^3$.

Configurare LSE

Utilizziamo setup_static_lse da qiskit_addon_mpf.static per costruire la matrice $A$ e il vettore del lato destro $b$ descritti sopra. La matrice $A$ dipende non solo da $k_j$, ma anche dalla formula del prodotto che scegliamo — in particolare dal suo ordine $\chi$ e dal fatto che sia simmetrica o meno. Il symmetric flag controlla lo schema dell'esponente $\eta_n$ (le formule simmetriche producono solo termini di errore di Trotter di potenza pari; cfr. rif. [1] ). Si noti che, come illustrato nel rif. [2], l'impostazione symmetric=True di non è strettamente necessaria nemmeno quando la funzione di prestazione (PF) sottostante è simmetrica: la LSE non simmetrica rimane valida (impone ulteriori vincoli non necessari).

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

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

Controllare la matrice $A$ e il vettore $b$ per verificare che corrispondano al sistema riportato sopra.

lse.A

Output:

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

lse.b

Output:

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

Una volta ottenuta l'equazione LSE, calcoliamo 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}"
)

Output:

The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]

Ottimizzare per l' $x$ e utilizzando un modello esatto

In alternativa al calcolo di $x = A^{-1}b$, è possibile utilizzare setup\_exact\_model per costruire un'istanza di cvxpy.Problem che utilizzi l'LSE come vincoli e la cui soluzione ottimale fornirà $x$.

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

Output:

[ 0.04761905 -0.57142857  1.52380952]

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

Output:

L1 norm of the exact coefficients: 2.1428571428556378

Ottimizzazione per l' $x$ e utilizzando un modello approssimativo

Potrebbe accadere che la norma $L_1$ per l'insieme scelto di valori $k_j$ sia ritenuta troppo elevata. Se è così e non è possibile scegliere un diverso insieme di valori per $k_j$, è possibile utilizzare una soluzione approssimativa che vincoli la norma $L_1$ a una soglia prescelta, minimizzando al contempo $\|Ax - b\|$. Consulta la guida su Come utilizzare il modello approssimativo.

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

Output:

[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

Coefficienti MPF dinamici

L'MPF statico annulla i termini di errore di Trotter in modo indipendente dall'hamiltoniano e dallo stato, pertanto non produce necessariamente l'errore di approssimazione più piccolo possibile per un dato hamiltoniano e uno stato iniziale. L’MPF dinamico (Rif. [2], [3] ) individua invece coefficienti dipendenti dal tempo $x_i(t)$ che minimizzano la distanza della norma di Frobenius $\|\rho(t) - \mu^D(t)\|_F^2$ in ogni istante $t$. Come illustrato nella sezione “Contesto”, ciò richiede la matrice di sovrapposizione $M_{ij}(t)$ tra gli stati evoluti secondo Trotter e la sovrapposizione $L_i(t)$ con lo stato esatto — entrambe le quali stimiamo utilizzando backend basati su reti tensoriali ( TeNPy ) in qiskit_addon_mpf.

Per configurare l'LSE dinamico occorrono tre elementi:

1. Una factory di evolver approssimativa che l'add-on eseguirà per ogni $k_j$ per generare $\rho_{k_j}(t)$ come MPS/MPO. Lo costruiamo a partire dalla struttura a strati del circuito di Trotter dell’ordine $2$ (uno strato per slice_by_depth), avvolto come un LayerwiseEvolver con parametri di troncamento TeNPy.
2. Una funzione di evoluzione esatta che produce un riferimento ad alta precisione $\rho(t)$. Utilizziamo un circuito di Suzuki-Trotter di quarto ordine a passo temporale piccolo (dt=0.1, order=4) come approssimazione dell'evoluzione esatta.
3. Una “fabbrica di identità” e un MPS di stato iniziale che fungono da seed per la simulazione “ TeNPy ”.

La cella sottostante crea la factory dell'evolver approssimativo.

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

dt=0.1Approssimiamo lo stato esatto in funzione del tempo con una formula di Suzuki-Trotter di quarto ordine, utilizzando un passo temporale piccolo. I parametri di troncamento dell' TeNPy e possono influire sulla precisione, pertanto è importante valutare una serie 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 un identity_factory che dia come risultato lo stato MPO iniziale e prepariamo lo stato iniziale di Néel come un MPS che corrisponda 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)

Una volta definite le fabbriche, calcoliamo ora i coefficienti dinamici in ciascun istante di evoluzione. Per ogni $t$, setup_dynamic_lse costruisce le matrici di sovrapposizione pertinenti tramite TeNPy, e setup_frobenius_problem restituisce un cvxpy.Problem che minimizza il costo della norma di Frobenius. Il risolutore restituisce i coefficienti $x_j(t)$ specifici per quel periodo; li raccogliamo in mpf_dynamic_coeffs_list. Se il risolutore non riesce a risolvere un dato $t$, si torna ai coefficienti pari a zero 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("")

Output:

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 attesi di Trotter con i coefficienti MPF

Ora calcoliamo il valore di " $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ " per ciascun insieme di coefficienti (static-exact, static-approximate e dynamic), propaghiamo gli errori standard per ciascun circuito e tracciamo le serie temporali risultanti 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()

Output:

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

• Errore di Trotter. Le formule dei singoli prodotti (indicatori grigi) si discostano sempre più dalla curva esatta con il passare del tempo. Il circuito $k=1$ presenta la deviazione maggiore ed è quello con la profondità minore, ma si trova già nel regime in cui $t/k \gtrsim 1$, quindi il termine di errore principale $1/k^{2}$ è elevato. Le combinazioni MPF (indicatori colorati) annullano molti di questi termini di errore di Trotter principali, quindi seguono la curva esatta in modo molto più fedele rispetto a qualsiasi singolo circuito di " $k_j$ ". Il divario residuo riflette i termini di Trotter di ordine superiore che l’MPF non annulla: un MPF statico di ordine $2$, $r=3$ elimina solo i primi due ordini di errore e, a valori elevati di $t/k_{\min}$, la coda non annullata finisce per prevalere — pertanto l’MPF non garantisce che i circuiti molto poco profondi rimangano accurati in momenti arbitrari.

• Errore di campionamento. Le barre di errore più ampie sulle curve MPF sono una conseguenza diretta della combinazione lineare: propagando gli errori standard indipendenti per ciascun circuito $\sigma_{k_j}$ si ottiene una varianza totale $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Pertanto, quanto maggiore è l’ $\|x\|_2$ (e, in pratica, l’ $\|x\|_1$, che è ciò che controlliamo), tanto più sono necessari per raggiungere una data incertezza target. Questo è il compromesso alla base dell’opzione “solutore approssimativo” in “Background”: limitiamo il valore di $\|x\|_1$ per mantenere questo sovraccarico a livelli gestibili. È fondamentale notare che, a differenza dell’errore di Trotter, l’errore di campionamento si riduce con l’aumentare dell’ $1/\sqrt{N_{\text{shots}}}$ e, quindi può sempre essere ridotto effettuando un numero maggiore di misurazioni.

Nell'esempio di hardware su larga scala riportato di seguito, il rumore dell'hardware costituisce una fonte di errore aggiuntiva su ogni $\langle A \rangle_{k_j}$, che viene a sua volta amplificata dai coefficienti MPF. In quella sezione vedremo come la mitigazione degli errori interagisce con gli MPF.

Esempio di hardware su larga scala

In questa sezione estendiamo il problema oltre i limiti di ciò che è possibile simulare con precisione. Riproduciamo alcuni dei risultati riportati nel riferimento [3], utilizzando una catena XXZ da 50 qubit al tempo $t = 3$. Seguiamo lo stesso flusso di lavoro in quattro fasi dell’esempio su piccola scala, puntando ora a hardware quantistico reale dotato di mitigazione degli errori. Come nel modello, ogni fase è contrassegnata direttamente nel codice, e una singola fase può estendersi su più celle quando vale la pena esaminarne i risultati intermedi.

La mappatura rispecchia l'esempio su piccola scala: definire un hamiltoniano, scegliere i parametri di Trotter, calcolare i coefficienti MPF (statici e dinamici) e costruire i circuiti. Le differenze principali sono:

• Un hamiltoniano XXZ su 50 siti con accoppiamenti casuali, tratto da $\mathcal{U}(0.5, 1.5)$ (Rif. [3] ).
• Una formula di Trotter simmetrica di secondo ordine con $k_j = [3, 4, 6]$ (quindi $\chi=1$, symmetric=True).
• Un unico tempo di evoluzione fisso $t = 3$. Con $k_{\min}=3$ si ottiene $t/k_{\min}=1$, mantenendo le componenti di bassa profondità all’interno del regime di convergenza di Trotter, dove è valido il modello dell’errore dominante su cui si basa l’MPF.
• Un'ulteriore serie di confronti a circuito singolo con i passi di Trotter dell' $k = 10$, utilizzati come riferimento. Abbiamo scelto l’ $k = 10$ perché la sua profondità di due qubit sull’hardware è superiore a quella del costituente MPF più profondo ( $k_{\max}=6$ ) più il sovraccarico derivante dall’esecuzione di più circuiti MPF — sufficientemente profonda da essere limitata dal rumore, ovvero il regime in cui ci si aspetta che la combinazione MPF superi le prestazioni del circuito singolo di riferimento. Si tratta di un confronto basato su un "circuito singolo profondo" rispetto alla combinazione MPF, non di un circuito che miri all'errore di Trotter effettivo dell'MPF (il che richiederebbe molti più passaggi).

Si noti che, sebbene ci troviamo ancora nella Fase 1 (mappatura e costruzione del circuito), in questa cella calcoliamo in anticipo sia i coefficienti dinamici che quelli statici. I coefficienti dinamici dipendono da $H$ e $t$, ma non dalle misurazioni quantistiche; pertanto, possono essere calcolati in qualsiasi momento prima della Fase 4. Lo facciamo ora per tenere tutte le impostazioni specifiche dell'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)

Output:

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. optimization_level=3Utilizziamo il gestore di passaggi predefinito di Qiskit, che seleziona automaticamente un insieme ottimale di qubit fisici e instrada ciascun 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
]

Output:

<IBMBackend('ibm_fez')>

L'esecuzione di circuiti più complessi su hardware reale richiede misure aggressive di mitigazione degli errori. Consentiamo il disaccoppiamento dinamico, la rotazione dei gate e delle misurazioni, la mitigazione degli errori di misurazione e l'estrapolazione a rumore zero (ZNE). Si noti che i fattori di rumore ZNE qui utilizzati (1, 1.2, 1.4) sono inferiori rispetto a quelli di uno scenario a circuito superficiale, poiché i costituenti MPF più profondi sono già vicini alla soglia di rumore e forti amplificazioni del rumore li spingerebbero oltre il punto in cui l’estrapolazione ZNE è affidabile.

Inviamo tutti e quattro i circuiti (tre componenti MPF all’indirizzo $k_j = [3, 4, 6]$ più il valore di riferimento $k = 10$ ) in un unico processo di 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 attesi e le deviazioni standard per ciascun circuito dai risultati del calcolo, quindi li combiniamo con ciascuna serie 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 la 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)

Output:

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

Output:

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

Output:

Alcune osservazioni sui risultati relativi all'hardware riportati sopra:

• L'approfondimento non è gratuito a livello hardware. I grafici di riferimento a circuito singolo parlano chiaro: il circuito $k = 6$ è sostanzialmente esatto ( $-0.256$ rispetto al riferimento $-0.244$ ), mentre il grafico di riferimento più approfondito $k = 10$ è peggiore ( $-0.061$, con uno scostamento di $\sim 0.18$ ), non migliore. Una volta che l’errore di Trotter è già ridotto, l’aggiunta di ulteriori passi non fa altro che aumentare la profondità del circuito e accumulare ulteriore rumore di gate e decoerenza. È proprio questo il contesto per cui sono stati concepiti gli MPF: raggiungere la precisione di un circuito profondo utilizzando solo componenti superficiali.

• Un MPF a norma ridotta supera il circuito singolo profondo. L'MPF approssimativo-statico (con un limite massimo di $\|x\|_1 \approx 2$ ) si attesta a $-0.259$, a soli $\sim 0.015$ dal valore di riferimento e molto più vicino rispetto al valore di base di $k = 10$. Anche il dinamico MPF ( $-0.127$ ) supera ampiamente tale valore di riferimento. Entrambe combinano solo i circuiti superficiali $k_j = [3, 4, 6]$, ma riescono comunque a ricavare una risposta che il singolo circuito profondo non era in grado di fornire.

• La norma del coefficiente è più importante dell'ottimalità matematica. L'MPF statico esatto presenta un $\|x\|_1 = 4.66$ e ed è il peggiore stimatore in assoluto ( $-0.567$, con uno scostamento superiore a $0.3$ ): l'elevata norma dei coefficienti amplifica il rumore residuo del gate, la decoerenza e l'errore ZNE su ciascun $\langle A \rangle_{k_j}$ all'incirca dello stesso fattore, vanificando la cancellazione dell'errore di Trotter che ne deriva. L'applicazione di un limite massimo alla norma (il risolutore approssimativamente statico, $\|x\|_1 \approx 2$ ) elimina questo sovraccarico e fornisce la stima migliore — anche se i suoi coefficienti non annullano più esattamente l'errore di Trotter principale.

• Anche i singoli circuiti poco profondi possono comunque essere competitivi. L'unico componente " $k = 6$ " ( $-0.256$ ) è di per sé sostanzialmente esatto in questo caso — in questa simulazione risulta addirittura leggermente più preciso dell'MPF "approximate-static". Il problema è che non si sa in anticipo quale singolo $k$ si trovi nel punto ottimale in cui “la convergenza è raggiunta ma non è ancora limitata dal rumore”, e la scelta apparentemente sicura di limitarsi ad aumentare la profondità ( $k = 10$ ) per garantire la convergenza di Trotter è proprio quella che fallisce. L'MPF offre una combinazione basata su principi di circuiti a bassa profondità che non richiede di indovinare la profondità corretta.

In pratica, ciò significa che, a livello hardware, gli MPF dovrebbero essere abbinati a una forte mitigazione degli errori su ogni singolo $\langle A \rangle_{k_j}$, la norma del coefficiente $L_1$ dovrebbe essere mantenuta modesta (utilizzare il risolutore approssimativo o l’MPF dinamico) e i passi di Trotter $k_j$ dovrebbero essere scelti in modo tale che $t/k_{\min} \lesssim 1$ — qui $k_{\min} = 3$ all’indirizzo $t = 3$ fornisce $t/k_{\min} = 1$, mantenendo i componenti all’interno del regime di convergenza in cui è valido il modello di errore principale su cui si basa l’MPF statico. Con queste scelte, gli MPF a norma piccola qui considerati eguagliano un singolo circuito convergente, mentre la linea di base “naïve” (“basta andare più in profondità”) non ci riesce, ripristinando il vantaggio in termini di profondità rispetto all’accuratezza illustrato nel rif. [3]. Si noti inoltre che le singole esecuzioni sono soggette a rumore: in un’altra esecuzione dello stesso lavoro (o su un backend diverso), l’ordine esatto può variare; le tendenze evidenti sono che gli MPF small- $\|x\|_1$ i ottengono buoni risultati, mentre l’MPF large- $\|x\|_1$ e exact-static è amplificato dal rumore dell’hardware e il circuito singolo eccessivamente profondo è limitato dal rumore.

Passi successivi

Riferimenti

[1] Vázquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Formule multiprodotto ben condizionate per la simulazione hamiltoniana ottimizzata per l'hardware. Quantum, 7, 1067 (2023)

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Limiti di errore di Trotter e formule dinamiche multiprodotto per la simulazione hamiltoniana. Physical Review Research, 6(3), 033309 (2024)

[3] Robertson, N. F., et al. Formule dinamiche multiprodotto potenziate tramite rete tensoriale. arXiv:2407.17405 (2024)