Reducerea erorii Trotter în dinamica hamiltoniană cu formule multi-produs

În acest notebook, vei învăța cum să folosești o Formulă Multi-Produs (MPF) pentru a obține o eroare Trotter mai mică asupra observabilului nostru față de cea introdusă de cel mai adânc Circuit Trotter pe care îl vom executa efectiv. Vei face asta parcurgând pașii unui pattern Qiskit:

• Pasul 1: Maparea la problema cuantică
  • Inițializarea Hamiltonianului problemei noastre
  • Folosirea unui MPF pentru a genera Circuit-urile de evoluție temporală Trotterizate
• Pasul 2: Optimizarea problemei
  • Aici transpilăm Circuit-urile noastre pentru un GenericBackendV2
• Pasul 3: Executarea experimentelor
  • Folosim un StatevectorEstimator pentru simplitate în acest notebook
• Pasul 4: Reconstruirea rezultatelor
  • Calcularea valorii așteptate a MPF

Pasul 1: Maparea la problema cuantică

1a: Configurarea Hamiltonianului nostru

Folosim modelul Ising pe o linie de 10 site-uri:

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

unde $J$ este intensitatea cuplajului dintre două site-uri și $h$ este câmpul magnetic extern. Pachetul qiskit_addon_utils oferă unele funcționalități reutilizabile pentru diverse scopuri.

Modulul său qiskit_addon_utils.problem_generators oferă funcții pentru generarea de Hamiltonieni de tip Heisenberg pe un graf de conectivitate dat. Acest graf poate fi fie un rustworkx.PyGraph, fie un CouplingMap, ceea ce îl face ușor de utilizat în fluxuri de lucru centrate pe Qiskit.

În cele ce urmează, creăm o linie simplă de 10 Qubiți folosind metoda 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")

În continuare, generăm SparsePauliOp pe conectivitatea furnizată cu constantele dorite.

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

Observabilul pe care îl vom măsura este magnetizarea totală, pe care o putem construi simplu după cum se arată mai jos:

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

MPF-urile reduc eroarea Trotter în dinamica hamiltoniană printr-o combinație ponderată a mai multor execuții de Circuit.

Pentru a fi mai concreți, definim un MPF astfel:

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

unde $x_j$ sunt coeficienții noștri de ponderare, $\rho^{k_j}_j$ este matricea densitate corespunzătoare stării pure obținute prin evoluția stării inițiale cu formula produs, $S^{k_j}$, implicând $k_j$ pași Trotter, iar $j$ indexează numărul de formule produs care alcătuiesc MPF-ul.

Ideea esențială este că eroarea Trotter rămasă este mai mică decât eroarea Trotter pe care ai obține-o folosind pur și simplu cea mai mare valoare $k_j$!

Poți privi utilitatea MPF din două perspective:

1. Pentru un buget fix de pași Trotter pe care îi poți executa, poți obține rezultate cu o eroare Trotter mai mică în total.
2. Pentru un număr de pași Trotter care duce la Circuit-uri adânci, poți folosi MPF pentru a găsi mai multe Circuit-uri cu adâncime mai mică pe care să le rulezi și care să conducă la o eroare Trotter similară.

O introducere în MPF-urile statice

MPF-urile statice sunt cele în care valorile $x_j$ NU depind de timpul de evoluție, $t$.

Determinarea coeficienților MPF statici pentru un set dat de valori $k_j$ se reduce la rezolvarea unui sistem liniar de ecuații: $Ax=b$, unde $x$ sunt coeficienții de interes, $A$ este o matrice care depinde de $k_j$ și de tipul de PF folosit ($S$), iar $b$ este un vector de constrângeri. Pe scurt, nu vom intra în mai multe detalii aici și te îndreptăm în schimb spre documentația LSE.

Putem găsi o soluție pentru $x$ analitic ca $x = A^{-1}b$, vezi de ex. Carrera Vazquez et al., 2023 sau Zhuk et al., 2023. Totuși, această soluție exactă poate fi „prost condiționată", rezultând în norme L1 foarte mari ale coeficienților noștri $x$, ceea ce poate duce la o performanță slabă a MPF. În schimb, se poate obține și o soluție aproximativă care minimizează norma L1 a lui $x$, în încercarea de a optimiza comportamentul MPF.

În cele ce urmează, vei învăța cum să faci toate acestea.

Alegerea $k_j$

Alegerea $k_j$ îi revine utilizatorului final. În principiu, orice valori pot fi alese, dar unele $k_j$-uri vor duce la o amplificare mai mare a zgomotului pe dispozitivele reale față de alte alegeri. Prin urmare, este important să încerci să găsești valori $k_j$ „bune".

Aici vom alege pur și simplu niște valori fixe pentru $k_j$. Valoarea minimă este motivată de timpul țintă de evoluție $t=8.0$, care în mod normal ne spune să satisfacem $t/k_{\text{min}} \lt 1$, dar empiric știm că setarea egală cu $1$ funcționează de obicei și ea. Dacă vrei să afli mai multe despre asta și despre cum să alegi celelalte valori $k_j$, consultă ghidul corespunzător: How to choose the Trotter steps for an MPF.

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

Configurarea LSE

Acum că am ales $k_j$-urile noastre, trebuie mai întâi să construim LSE, $Ax=b$, după cum s-a explicat mai sus. Matricea $A$ nu depinde doar de $k_j$, ci și de alegerea formulei produs (PF) -- în special de ordinul acesteia. În plus, se poate lua în considerare dacă PF este simetrică sau nu (vezi Carrera Vazquez et al., 2023), setând symmetric=True. Totuși, acest lucru nu este obligatoriu, după cum arată Zhuk et al., 2023.

Aici vom folosi o formulă Suzuki-Trotter de ordinul doi, obținând order=2, și vom seta 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.]))

Rezolvarea analitică pentru $x$

Așa cum am menționat anterior, putem găsi $x$ analitic:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

Optimizarea pentru $x$ folosind un model exact

Ca alternativă la calculul $x=A^{-1}b$, poți folosi și setup_exact_problem pentru a construi o instanță cvxpy.Problem care folosește LSE drept constrângeri și a cărei soluție optimă va produce $x$.

În secțiunea următoare va fi clar de ce există această interfață.

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]

Ca indicator al faptului că un MPF construit cu acești coeficienți va produce rezultate bune, putem folosi norma L1 (vezi și Carrera Vazquez et al., 2023).

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

3.3889400921655914

Optimizarea pentru $x$ folosind un model aproximativ

Se poate întâmpla ca norma L1 pentru setul ales de valori $k_j$ să fie considerată prea mare. Dacă acesta este cazul și nu poți alege un alt set de valori $k_j$, poți folosi o soluție aproximativă la LSE în loc de una exactă.

Pentru aceasta, folosește pur și simplu setup_sum_of_squares_problem pentru a construi o altă instanță cvxpy.Problem care constrânge norma L1 la un prag ales, minimizând în același timp diferența dintre $Ax$ și $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

Reține că ai libertate deplină în privința modului de rezolvare a acestei probleme de optimizare, adică poți schimba solver-ul de optimizare, pragurile sale de convergență și altele. Consultă ghidul corespunzător despre How to use the approximate model.

1c: Configurarea circuitelor Trotter

În acest punct, am găsit coeficienții de expansiune, $x$, și tot ce rămâne de făcut este să generăm circuitele cuantice Trotterizate. Din nou, modulul qiskit_addon_utils.problem_generators vine în ajutor pentru a face exact asta:

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)

Pasul 2: Optimizarea problemei

În mod normal, acesta este pasul din tipar în care îți optimizezi circuitele pentru execuția pe hardware. Aici, deoarece folosim doar un simulator fără zgomot, transpilăm pur și simplu circuitul pentru 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]

Pasul 3: Executarea experimentelor cuantice

Așa cum am explicat la început, vom sări peste pasul 2 de optimizare deoarece vom calcula pur și simplu valorile așteptate ale observabilului țintă folosind un simulator fără zgomot, și anume StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

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

Pasul 4: Reconstituirea rezultatelor

Mai întâi, extragem valorile așteptate individuale obținute pentru fiecare dintre circuitele Trotter:

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

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

Apoi, le recombinăm pur și simplu cu coeficienții noștri MPF pentru a obține valorile așteptate totale ale MPF. Mai jos, facem acest lucru pentru fiecare dintre diferitele moduri prin care am calculat $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

În final, pentru această problemă mică putem calcula valoarea exactă de referință folosind scipy.linalg.expm după cum urmează:

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

Putem vedea clar că MPF a redus eroarea Trotter față de cea obținută cu cel mai adânc PF individual cu $k_j=19$. Totuși, observăm și că modelul aproximativ nu este perfect, deoarece a rezultat de fapt într-o valoare așteptată mai proastă decât soluția exactă. Aceasta arată importanța utilizării unor criterii de convergență stricte pentru modelul aproximativ, așa cum vei afla în ghidul Cum să folosești modelul aproximativ.