Ridurre l'errore di Trotter della dinamica hamiltoniana con le formule multi-prodotto

In questo notebook imparerai come usare una Multi-Product Formula (MPF) per ottenere un errore di Trotter inferiore sul nostro osservabile rispetto a quello causato dal circuito di Trotter più profondo che eseguiremo effettivamente. Lo farai seguendo i passi di un pattern Qiskit:

• Passo 1: Mappare al problema quantistico
  • Inizializzare l'Hamiltoniano del nostro problema
  • Usare un MPF per generare i circuiti di evoluzione temporale Trotterizzati
• Passo 2: Ottimizzare il problema
  • Qui traspiliamo i nostri circuiti per un GenericBackendV2
• Passo 3: Eseguire gli esperimenti
  • Usare uno StatevectorEstimator per semplicità in questo notebook
• Passo 4: Ricostruire i risultati
  • Calcolare il valore di aspettazione dell'MPF

Passo 1: Mappare al problema quantistico

1a: Impostare il nostro Hamiltoniano

Usiamo il modello di Ising su una catena di 10 siti:

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

dove $J$ è la forza di accoppiamento tra due siti e $h$ è il campo magnetico esterno. Il pacchetto qiskit_addon_utils fornisce alcune funzionalità riutilizzabili per vari scopi.

Il suo modulo qiskit_addon_utils.problem_generators fornisce funzioni per generare Hamiltoniani simili a Heisenberg su un dato grafo di connettività. Questo grafo può essere un rustworkx.PyGraph oppure una CouplingMap, il che lo rende facile da usare nei flussi di lavoro incentrati su Qiskit.

Di seguito creiamo una semplice catena di 10 qubit usando il metodo CouplingMap.from_line.

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

from qiskit.transpiler import CouplingMap

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

from rustworkx.visualization import graphviz_draw

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

Successivamente generiamo lo SparsePauliOp sulla connettività fornita con le costanti desiderate.

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

# Get a qubit operator describing the Ising field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(0.0, 0.0, 1.0),
    ext_magnetic_field=(0.4, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],
              coeffs=[1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j,
 1. +0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,
 0.4+0.j, 0.4+0.j, 0.4+0.j])

L'osservabile che misureremo è la magnetizzazione totale, che possiamo costruire semplicemente come mostrato di seguito:

from qiskit.quantum_info import SparsePauliOp

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list([("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L)
print(observable)

SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j,
 0.05+0.j, 0.05+0.j, 0.05+0.j])

1b: Formule Multi-Prodotto

Le MPF riducono l'errore di Trotter della dinamica hamiltoniana attraverso una combinazione pesata di diverse esecuzioni di circuiti.

Per rendere questo più concreto, definiamo un MPF come:

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

dove $x_j$ sono i nostri coefficienti di ponderazione, $\rho^{k_j}_j$ è la matrice densità corrispondente allo stato puro ottenuto evolvendo lo stato iniziale con la formula prodotto, $S^{k_j}$, che coinvolge $k_j$ passi di Trotter, e $j$ indicizza il numero di formule prodotto che compongono l'MPF.

Il punto chiave è che l'errore di Trotter residuo è inferiore all'errore di Trotter che si otterrebbe semplicemente usando il valore $k_j$ più grande!

Puoi vedere l'utilità dell'MPF da due prospettive:

1. Per un budget fisso di passi di Trotter che puoi eseguire, puoi ottenere risultati con un errore di Trotter complessivamente inferiore.
2. Per un numero di passi di Trotter che produce circuiti profondi, puoi usare l'MPF per trovare diversi circuiti di profondità minore da eseguire che risultano in un errore di Trotter simile.

Introduzione alle MPF statiche

Le MPF statiche sono quelle in cui i valori $x_j$ NON dipendono dal tempo di evoluzione, $t$.

Determinare i coefficienti dell'MPF statica per un dato insieme di valori $k_j$ equivale a risolvere un sistema lineare di equazioni: $Ax=b$, dove $x$ sono i nostri coefficienti di interesse, $A$ è una matrice che dipende da $k_j$ e dal tipo di PF che usiamo ($S$), e $b$ è un vettore di vincoli. Per brevità, non approfondiremo ulteriormente e ti rimandiamo invece alla documentazione di LSE.

Possiamo trovare una soluzione per $x$ analiticamente come $x = A^{-1}b$, vedi ad esempio Carrera Vazquez et al., 2023 o Zhuk et al., 2023. Tuttavia, questa soluzione esatta può essere "mal condizionata" risultando in norme L1 molto grandi dei nostri coefficienti, $x$, il che può portare a prestazioni scarse dell'MPF. In alternativa, si può ottenere una soluzione approssimata che minimizza la norma L1 di $x$ al fine di tentare di ottimizzare il comportamento dell'MPF.

Di seguito imparerai come fare tutto questo.

Scegliere $k_j$

La scelta di $k_j$ spetta all'utente finale. In linea di principio si possono scegliere qualsiasi valore, ma alcuni $k_j$ portano a un'amplificazione del rumore maggiore sui dispositivi reali rispetto ad altre scelte. Pertanto è importante che si cerchi di trovare valori "buoni" di $k_j$.

Qui sceglieremo semplicemente alcuni valori fissi per $k_j$. Il valore minimo è motivato dal tempo di evoluzione target di $t=8.0$ che normalmente ci dice di soddisfare $t/k_{\text{min}} \lt 1$, ma empiricamente sappiamo che impostarlo uguale a $1$ di solito funziona anch'esso. Se vuoi saperne di più su questo argomento e su come scegliere gli altri valori di $k_j$, consulta la guida corrispondente: How to choose the Trotter steps for an MPF.

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

Impostare l'LSE

Ora che abbiamo scelto i nostri $k_j$, dobbiamo prima costruire l'LSE, $Ax=b$ come spiegato sopra. La matrice $A$ dipende non solo da $k_j$ ma anche dalla nostra scelta di formula prodotto (PF) -- in particolare dal suo ordine. Inoltre, si può tener conto se la PF è simmetrica o meno (vedi Carrera Vazquez et al., 2023), impostando symmetric=True. Tuttavia ciò non è necessario come mostrato da Zhuk et al., 2023.

Qui utilizzeremo una formula di Suzuki-Trotter del secondo ordine che produce order=2 e imposteremo symmetric=True.

from qiskit_addon_mpf.static import setup_static_lse

lse = setup_static_lse(trotter_steps, order=2, symmetric=True)
print(lse)

LSE(A=array([[1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
       [1.56250000e-02, 6.94444444e-03, 2.77008310e-03],
       [2.44140625e-04, 4.82253086e-05, 7.67336039e-06]]), b=array([1., 0., 0.]))

Risolvere analiticamente per $x$

Come accennato in precedenza, possiamo trovare $x$ analiticamente:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

Ottimizzare per $x$ usando un modello esatto

In alternativa al calcolo di $x=A^{-1}b$, puoi anche usare setup_exact_problem per costruire un'istanza di cvxpy.Problem che usa l'LSE come vincoli e la cui soluzione ottimale produrrà $x$.

Nella prossima sezione sarà chiaro perché questa interfaccia esiste.

from qiskit_addon_mpf.costs import setup_exact_problem

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

[ 0.17239057 -1.19447005  2.02207947]

Come indicatore del fatto che un MPF costruito con questi coefficienti produrrà buoni risultati, possiamo usare la norma L1 (vedi anche Carrera Vazquez et al., 2023).

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

3.3889400921655914

Ottimizzare per $x$ usando un modello approssimato

Può capitare che la norma L1 per l'insieme scelto di valori $k_j$ sia ritenuta troppo alta. In tal caso, se non puoi scegliere un insieme diverso di valori $k_j$, puoi usare una soluzione approssimata all'LSE anziché una esatta.

Per farlo, usa semplicemente setup_sum_of_squares_problem per costruire un'istanza diversa di cvxpy.Problem che vincola la norma L1 a una soglia scelta minimizzando al contempo la differenza di $Ax$ e $b$.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

model_approx, coeffs_approx = setup_sum_of_squares_problem(lse, max_l1_norm=3.0)
model_approx.solve()
print(coeffs_approx.value)
print(np.linalg.norm(coeffs_approx.value, ord=1))

[-0.40454257  0.57553173  0.8290123 ]
1.8090865903790838

Nota che hai completa libertà su come risolvere questo problema di ottimizzazione, il che significa che puoi cambiare il risolutore di ottimizzazione, le sue soglie di convergenza e così via. Consulta la guida corrispondente su How to use the approximate model.

1c: Impostare i circuiti di Trotter

A questo punto abbiamo trovato i nostri coefficienti di espansione, $x$, e tutto ciò che resta da fare è generare i circuiti quantistici Trotterizzati. Ancora una volta, il modulo qiskit_addon_utils.problem_generators viene in soccorso proprio per questo:

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

circuits = []
for k in trotter_steps:
    circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=2, reps=k),
        time=time,
    )
    circuits.append(circ)

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

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

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

Passo 2: Ottimizzare il problema

Normalmente questo è il passo del pattern durante il quale ottimizzi i tuoi circuiti per l'esecuzione su hardware. Qui, poiché usiamo solo un simulatore senza rumore, ci limitiamo a traspilare il nostro circuito per un GenericBackendV2.

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

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

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

Passo 3: Eseguire gli esperimenti quantistici

Come spiegato all'inizio, salteremo il passo di ottimizzazione 2 perché calcoleremo semplicemente i valori di aspettazione del nostro osservabile target usando un simulatore senza rumore, ovvero lo StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

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

Passo 4: Ricostruire i risultati

Prima estraiamo i singoli valori di aspettazione ottenuti per ciascuno dei circuiti di Trotter:

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

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

Successivamente li ricombiniamo semplicemente con i nostri coefficienti MPF per ottenere i valori di aspettazione totali dell'MPF. Di seguito lo facciamo per ciascuno dei diversi modi con cui abbiamo calcolato $x$.

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

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

Infine, per questo piccolo problema possiamo calcolare il valore di riferimento esatto usando scipy.linalg.expm come segue:

from scipy.linalg import expm

exp_H = expm(-1j * time * hamiltonian.to_matrix())

initial_state = np.zeros(exp_H.shape[0])
initial_state[0] = 1.0

time_evolved_state = exp_H @ initial_state

exact_obs = time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
print(exact_obs.real)

0.40060242487899755

Possiamo vedere chiaramente che l'MPF ha ridotto l'errore di Trotter rispetto a quello ottenuto con la PF individuale più profonda con $k_j=19$. Tuttavia vediamo anche che il modello approssimato non è privo di difetti, poiché ha effettivamente prodotto un valore di aspettazione peggiore della soluzione esatta. Questo mostra l'importanza di usare criteri di convergenza stretti sul modello approssimato, come imparerai nella guida How to use the approximate model.