Generează circuite Krylov de evoluție în timp (SKQD)

Important

Conceptele din acest ghid sunt disponibile în prezent doar în API-ul Python. Funcționalitatea echivalentă va fi disponibilă în API-ul C într-o versiune viitoare.

Diagonalizarea cuantică Krylov bazată pe eșantionare (SKQD) este o variantă cuantic-centrică a diagonalizării în subspațiu. În loc să eșantioneze un singur ansatz variațional (ca în diagonalizarea cuantică bazată pe eșantionare simplă (SQD)) sau un ansamblu de circuite aleatorizate de evoluție în timp (ca în SqDRIFT), ea eșantionează o bază Krylov: familia de stări

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

obținute prin evoluția unei stări de referință $\lvert \psi_0 \rangle$ pentru multipli crescători ai unui pas de timp fix $\Delta t$. Șirurile de biți eșantionate din aceste $D$ circuite populează subspațiul Krylov în care Hamiltonianul este apoi diagonalizat clasic.

Acest ghid reproduce pasul de construcție a circuitelor din algoritmul SKQD pentru modelul Anderson cu o singură impuritate (SIAM), o impuritate magnetică cuplată la o baie neinteractivă, urmând publicația SKQD.

Notă

Pasul de eșantionare SKQD funcționează doar dacă starea fundamentală este rară în baza computațională: diagonalizarea clasică trebuie să o poată reconstrui dintr-un număr rezonabil de șiruri de biți eșantionate. După cum se arată mai jos, această raritate nu apare automat; este o proprietate a bazei de particulă unică în care este exprimată problema. Alegerea corectă a bazei este esența acestei construcții.

1. Construiește Hamiltonianul SIAM ca operator fermionic

Modelul Anderson cu o singură impuritate plasează un orbital de impuritate interactiv (indicele $0$) în contact cu o baie neinteractivă de $L$ situsuri aranjate într-un lanț. Hamiltonianul său se împarte într-o parte cu un corp $H_1$ și o parte cu doi corpi $H_2$,

$$H = H_1 + H_2,$$

unde partea cu un corp cuprinde energia locală a impurității (potențialul chimic $\mu$), hibridizarea impuritate–baie $V$ și saltul între vecinii cei mai apropiați din baie $t$,

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

iar partea cu doi corpi conține repulsia Coulomb locală $U$ doar pe impuritate,

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

Acest ghid lucrează în punctul simetric particulă-gol $\mu = -U/2$. Aceste două părți corespund direct integralelor cu un corp h1e (saltul în baie, hibridizarea impurității, energia locală a impurității) și singurei integrale cu doi corpi h2e (termenul Coulomb al impurității). Mai întâi, asamblează-le în baza pozițiilor, în care modelul este definit în mod natural.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. Trecerea la baza impulsului

Partea de baie a SIAM este un lanț de fermioni liberi, deci este diagonalizată printr-o trecere la baza impulsului. Acesta este pasul decisiv pentru SKQD. Deoarece baia este invariantă la translații, starea fundamentală nu este rară în baza pozițiilor, dar este rară în baza impulsului, unde Hamiltonianul băii este diagonal. Rotația orbitală de particulă unică ce realizează această trecere este algebră liniară specifică problemei. Ea diagonalizează matricea de salt a băii, lăsând neatins orbitalul impurității, apoi mută impuritatea într-un situs central:

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

Spre deosebire de un flux de lucru variațional, aici se rotesc integralele Hamiltonianului în baza impulsului, astfel încât întregul circuit (Hamiltonianul și starea deopotrivă) trăiește în baza în care starea fundamentală este rară, iar șirurile de biți eșantionate sunt, prin urmare, informative. Rotirea integralelor cu un corp și cu doi corpi este o schimbare de bază standard a tensorilor de structură electronică:

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

Deoarece rotația lasă orbitalul impurității fix (acționează doar asupra băii), interacțiunea locală rămâne neschimbată. h2e_momentum rămâne un singur termen număr-număr, exact ca h2e. Singura schimbare este că impuritatea a fost mutată la un indice central. Acest lucru se vede direct în tensorii de integrale, fiecare având o singură intrare nenulă:

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

Această localitate este ceea ce menține interacțiunea ieftin de sintetizat mai târziu.

Asamblează integralele cu un corp într-un FermionOperator cu ajutorul constructorului de integrale electronice from_1body_tril_spin_sym(), care le așteaptă împachetate în ordinea triunghiului inferior. Partea cu doi corpi este un singur termen număr-număr, așa că, în loc să o treci prin mecanismul integralelor cu doi corpi, construiește-o direct cu cre()/ann(): operatorul Coulomb local $U \, n_{p\uparrow} n_{p\downarrow}$ pentru modul de impuritate $p$ (modul său $\uparrow$ la indicele p, modul său $\downarrow$ la indicele p + norb în aranjamentul pe blocuri de spin). Ambii constructori dublează prin spin cei norb orbitali spațiali în 2 * norb moduri fermionice.

Păstrează partea cu un corp $H_1$ și partea cu doi corpi $H_2$ ca operatori separați, în loc să îi aduni imediat, în acord cu termenii $H_1$ și $H_2$ din ecuațiile secțiunii 1. Construirea fiecărei piese o singură dată (în ambele baze) le face reutilizabile mai jos. Adunați, ei dau Hamiltonianul complet pentru diagonalizarea exactă în oricare dintre baze. Individual, $H_2$ (baza impulsului) este operatorul pe care pasul Trotter îl evoluează la pasul 4.

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

Schimbarea de bază este o transformare de similitudine unitară, deci lasă spectrul neschimbat. Confirmă că operatorul este corect comparând cea mai mică valoare proprie a sa în sectorul $(\text{norb}, \text{nelec})$ cu o diagonalizare exactă. FermionOperator expune un LinearOperator din SciPy (susținut de un nucleu nativ de produs matrice-vector pentru interacția de configurații completă (FCI)), deci poate fi transmis direct către scipy.sparse.linalg.eigsh().

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

Câștigul schimbării de bază este raritatea. Pentru a o concretiza, numără câți determinanți din baza computațională sunt necesari pentru a capta 99% din ponderea stării fundamentale și compară baza impulsului cu baza pozițiilor. Energia nu depinde de bază, dar numărul de determinanți depinde:

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

Cei doi operatori au aceeași energie, dar starea fundamentală din baza pozițiilor este împrăștiată pe mii de determinanți, în timp ce baza impulsului o concentrează pe aproximativ 20, motiv pentru care eșantionarea SKQD funcționează în baza impulsului și nu în baza pozițiilor.

3. Pregătește starea de referință

Referința SKQD $\lvert \psi_0 \rangle$ este suprapunerea tuturor excitațiilor electronilor cei mai apropiați de nivelul Fermi către cele mai apropiate moduri de impuls goale. Publicația SKQD o construiește cu o rețea scrisă de mână de XXPlusYYGate. Pare puternic corelată (are sute de amplitudini nenule în baza computațională), dar o rețea de rotații cu doi moduri care conservă numărul de particule este o operație de particulă unică (gaussiană fermionică), deci starea pe care o produce este un singur determinant Slater. Asta înseamnă că o poți pregăti cu o singură poartă PrepareSlaterDeterminant, dată fiind rotația de particulă unică pe care o implementează rețeaua.

Deoarece fiecare XXPlusYYGate acționează ca o rotație Givens pe o pereche de moduri, acea rotație poate fi construită direct în imaginea de particulă unică (o matrice Givens pentru fiecare poartă), în loc să emiți manual porțile cu doi qubiți. Porțile formează o structură de tip zidărie a rotațiilor între vecini apropiați, care se desfășoară din nivelul Fermi în straturi tot mai largi, răspândind fiecare electron apropiat de nivelul Fermi pe modurile goale aflate imediat deasupra lui:

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

Notă

Rotația de mai sus este exprimată în baza impulsului, în acord cu Hamiltonianul. Electronii din apropierea nivelului Fermi sunt modurile de impuls ocupate cu cea mai mare energie, iar rețeaua îi excită în cele goale cu cea mai mică energie. Acesta este motivul pentru care referința trebuie construită după schimbarea de bază din pasul anterior, nu înainte.

Pregătirea determinantului Slater se face pentru fiecare sector de spin, deoarece fiecare sector are propria formă (electroni, orbitali). O singură poartă PrepareSlaterDeterminant (o ocupare, cele mai joase moduri nocc fiind umplute, împreună cu rotația de particulă unică) pregătește un sector; aceeași poartă se aplică o dată pe modurile alfa 0..norb și o dată pe modurile beta norb..2*norb, deoarece ambele sectoare de spin au aici aceeași ocupare de referință și aceeași rotație:

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

Această poartă declarativă înlocuiește rețeaua XXPlusYYGate scrisă manual; transpilatorul generează sinteza rotațiilor Givens.

4. Asamblează circuitele Krylov

Fiecare circuit Krylov este pregătirea stării de referință urmată de evoluția sub $e^{-i k \, \Delta t \, H}$. Împarte Hamiltonianul în părțile sale cu un corp și cu doi corpi, $H = H_1 + H_2$, și aplică descompunerea Trotter evoluției pe baza acestei împărțiri. Aceasta este alegerea decisivă pentru adâncimea circuitului. Evoluția cu un corp $e^{-i t H_1}$ este o operație de fermioni liberi (gaussiană fermionică), deci este o rotație de bază de particulă unică (un OrbitalRotation cu unitara $e^{-i t \, \mathbf{h}_1}$) și nu ceva ce trebuie descompus Trotter pe termeni Pauli, termen cu termen. Evoluția cu doi corpi $e^{-i t H_2}$ acționează doar asupra singurului termen de interacțiune locală. Construirea pasului în acest fel permite transpilatorului să emită o structură superficială de rotații Givens plus o singură poartă de interacțiune cu doi qubiți, structura superficială pe care publicația SKQD o construiește manual.

Avem deja partea cu doi corpi. Este operatorul two_body_momentum construit la pasul 2 (singurul termen local). Evoluția cu un corp $e^{-i t H_1}$ nu depinde de spin, deci este un singur OrbitalRotation pe sector, cu unitara $e^{-i t \, \mathbf{h}_1}$ aplicată o dată fiecărui sector de spin, aceeași structură paralelă ca la pregătirea stării de referință.

Fiecare circuit Krylov evoluează pentru un timp total $k \, \Delta t$. Aceasta se realizează cu un produs Trotter de ordinul doi din $k$ pași de mărime $\Delta t$, astfel încât eroarea pe pas rămâne fixă pe măsură ce dimensiunea Krylov crește. Fiecare pas include o rotație completă cu un corp (una pe sector de spin) între două jumătăți de pas ale evoluției cu doi corpi (ieftină, cu un singur termen). Dimensiunea Krylov $D$ este numărul acestor circuite, cu dim variind de la 0 (doar starea de referință) la $D - 1$.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

Notă

Circuitele construite aici nu au încă măsurători, deoarece aceleași circuite fermionice conduc și simularea exactă de la pasul 6, iar o măsurătoare nu este o operație unitară. Pasul 5 le adaugă cu FermionicCircuit.measure_all(), care citește modul i în bitul clasic i.

Preferă să măsori circuitul fermionic în locul celui transpilat. O măsurătoare fermionică numește modul pe care îl citește, deci amplasarea și rutarea o poartă împreună cu qubitul pe care a ajuns acel mod, iar șirul de biți eșantionat este indexat după mod. Dacă măsori în schimb circuitul transpilat, rămâne în sarcina ta să anulezi final_index_layout(), iar dacă uiți, rezultatele permutate sunt returnate în tăcere. Aceasta se bazează pe codificarea Jordan-Wigner în baza de ocupare folosită de preset, în care fiecare mod este purtat de un singur qubit.

5. Transpilează în circuite de qubiți

Transpilarea unui FermionicCircuit mapează porțile sale fermionice pe porți de qubiți. Aceasta folosește maparea Jordan-Wigner jordan_wigner(). Porțile OrbitalRotation se sintetizează într-o structură de tip zidărie de rotații Givens, iar Evolution cu un singur termen se reduce la o singură poartă de interacțiune cu doi qubiți.

Deoarece fiecare poartă fermionică de aici are o sinteză implicită, poți folosi direct generate_preset_jw_pass_manager(), gata făcut, în loc să asamblezi manual etapele. Orice argumente cu nume sunt transmise etapei de qubiți; transmite optimization_level=0 pentru o imagine fidelă, neoptimizată, a adâncimii sintetizate. Generarea întregii familii de circuite Krylov este apoi o buclă pe dimensiunea Krylov. Păstrează FermionicCircuit-urile netranspilate alături de circuitele de qubiți transpilate. Cele fermionice conduc simularea exactă (vector de stare) de la pasul 6, deci trebuie să rămână fără măsurători; measure_all(inplace=False) returnează o copie cu măsurători, pe care să o transpilezi, și lasă originalul neatins, iar acea copie este cea pe care ar rula-o un backend.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

Dimensiunea Krylov se vede direct în adâncimea circuitelor măsurată în porți cu doi qubiți. Pregătirea stării de referință este un cost fix, iar fiecare putere Krylov suplimentară adaugă încă un pas Trotter de ordinul doi al operatorului de evoluție în timp.

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

Notă

Împărțirea Hamiltonianului menține aceste adâncimi moderate. Evoluția cu un corp este o rotație orbitală autentică, deci se sintetizează într-o structură de tip zidărie de XXPlusYYGate pe qubiți adiacenți. Fiecare jumătate de pas a interacțiunii locale este o singură poartă RZZGate care cuplează cele două moduri de spin ale impurității; termenul ieftin pe care produsul de ordinul doi îl include în jurul rotației. Dacă am da întregul Hamiltonian unei singure porți Evolution, partea cu un corp ar fi descompusă Trotter pe termeni Pauli, termen cu termen, propagând lungi șiruri Z Jordan-Wigner peste cuplajele impuritate-baie (de rază lungă) și umflând adâncimea cu un ordin de mărime.

Primul circuit pregătește doar determinantul de referință, în timp ce ultimul poartă cei mai mulți pași Trotter (patru în acest caz):

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. Eșantionează șiruri de biți din circuitele Krylov

Pe hardware, ai rula circuitele transpilate și ai măsura. Aceasta este versiunea fără zgomot (vector de stare) a acelui pas, care arată șirurile de biți pe care le-ar colecta SKQD. Simulează fiecare FermionicCircuit netranspilat cu ffsim.apply_unitary() și eșantionează vectorul de stare rezultat cu ffsim.sample_state_vector().

O poartă PrepareSlaterDeterminant este de tip validează-apoi-rotește în simulare (verifică dacă starea de intrare ocupă determinantul ei de referință, apoi aplică rotația). Prin urmare, intrarea este determinantul de ocupare simplu: cele mai joase moduri nocc umplute pentru fiecare spin, fără rotație. Poarta aplică singură reference_rotation.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

Adunând eșantioanele din toate cele cinci circuite Krylov, suportul este foarte mic: câteva sute de șiruri de biți distincte din cei $\binom{8}{4}^2 = 4900$ de determinanți ai sectorului (4, 4). Această concentrare, o consecință directă a rarității din baza impulsului stabilite la pasul 2, face tratabilă diagonalizarea clasică ulterioară.

Reprezentarea grafică a numărătorilor cu qiskit.visualization.plot_histogram() face raritatea vizibilă. Câteva configurații domină, în frunte cu determinantul de referință.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. Diagonalizează în subspațiul eșantionat

Șirurile de biți eșantionate sunt intrarea pentru jumătatea clasică a SKQD. Hamiltonianul este proiectat pe subspațiul generat de configurațiile eșantionate și diagonalizat acolo. Acest pas este transmis pachetului qiskit-addon-sqd, a cărui funcție diagonalize_fermionic_hamiltonian() rulează întreaga buclă de diagonalizare cuantică bazată pe eșantionare. Ea construiește un subspațiu din loturi de configurații eșantionate, diagonalizează Hamiltonianul în el și rafinează iterativ subspațiul folosind recuperarea configurațiilor: inversarea ocupărilor pentru a repara configurațiile care încalcă simetria cunoscută a numărului de particule. Recuperarea configurațiilor este concepută pentru a anula daunele zgomotului hardware, care corupe șirurile de biți eșantionate în sectorul greșit al numărului de particule; eșantionarea de aici este fără zgomot, deci nu este nimic de reparat în această privință.

Solverul consumă numărătorile eșantionate sub forma unui BitArray, construit direct din counts agregate la pasul 6. Transmite num_bits în mod explicit: fără el, lățimea este dedusă din cel mai lat șir de biți eșantionat, care este mai îngust decât 2 * norb ori de câte ori nicio probă nu ocupă ultimul mod, iar solverul respinge apoi tabloul. Solverul diagonalizează în aceeași bază de impuls pe care o folosește restul construcției, așa că primește direct integralele h1e_momentum și h2e_momentum de la pasul 2, fără nicio reîmpachetare. Energia recuperată coincide cu referința exactă cu mult sub un mili-Hartree, pornind de la doar câteva sute de configurații eșantionate. Deoarece energia SQD este o margine superioară variațională pentru starea fundamentală reală, ea se află întotdeauna la sau deasupra lui reference_energy (încă în domeniu de la pasul 2):

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

Vizualizează rularea: convergența energiei de-a lungul iterațiilor și ocuparea medie a fiecărui orbital spațial în starea fundamentală recuperată. Energia pe iterație este cea mai mică dintre loturile acelei iterații, iar ocuparea însumează orbital_occupancies peste cele două sectoare de spin. Chiar și fără zgomot de corectat, energia scade de la o iterație la alta. Fiecare iterație diagonalizează în num_batches loturi de doar samples_per_batch configurații subeșantionate din grup, iar rafinarea iterativă datorată reportării șirurilor de biți îmbunătățește constant care configurații ajung în acele loturi, astfel că energia recuperată se îmbunătățește:

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

Aceasta încheie bucla SKQD. Circuitele din baza de impuls de la pașii 1–5 produc un eșantion rar (pasul 6), iar diagonalizarea bazată pe eșantion de mai sus transformă acel eșantion înapoi într-o energie a stării fundamentale.

Pașii următori

Acest ghid a rulat întregul flux SKQD pe un simulator de vector de stare fără zgomot. Pe hardware, circuits transpilate ar fi executate și măsurate în locul eșantionării vectorului de stare de la pasul 6, iar numărătorile rezultate ar alimenta diagonalizarea de la pasul 7 neschimbate. Consultă documentația Qiskit pentru ajutor la rularea circuitelor și tutorialele addon SQD pentru informații despre post-procesarea prin diagonalizare în subspațiu, inclusiv comportamentul ei pe eșantioane zgomotoase, unde recuperarea configurațiilor are cel mai mult de lucru.

Citește ghidul despre relația cu ffsim pentru a înțelege de ce ffsim.apply_unitary() și ffsim.sample_state_vector() funcționează nativ pe circuitele fermionice ale acestui pachet și de ce sectorul cu număr fix de particule folosit la pasul 6 face ca raritatea din baza de impuls stabilită la pasul 2 să fie eficientă la eșantionare.