Generare circuiti di evoluzione temporale di Krylov (SKQD)

Importante

I concetti di questa guida sono attualmente disponibili solo nell'API Python. Una funzionalità equivalente sarà resa disponibile nell'API C in una versione futura.

La diagonalizzazione quantistica di Krylov basata su campionamento (SKQD) è una variante quantum-centrica della diagonalizzazione di sottospazio. Invece di campionare un singolo ansatz variazionale (come nella semplice diagonalizzazione quantistica basata su campionamento (SQD)) o un insieme di circuiti di evoluzione temporale randomizzati (come in SqDRIFT), campiona una base di Krylov: la famiglia di stati

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

ottenuti facendo evolvere uno stato di riferimento $\lvert \psi_0 \rangle$ per multipli crescenti di un passo temporale fisso $\Delta t$. Le stringhe di bit campionate da questi $D$ circuiti popolano il sottospazio di Krylov in cui l'Hamiltoniano viene poi diagonalizzato in modo classico.

Questa guida riproduce il passaggio di costruzione dei circuiti dell'algoritmo SKQD per il modello di Anderson a singola impurezza (SIAM), un'impurezza magnetica accoppiata a un bagno non interagente, seguendo la pubblicazione di SKQD.

Nota

Il passaggio di campionamento di SKQD funziona solo se lo stato fondamentale è sparso nella base computazionale: la diagonalizzazione classica deve poterlo ricostruire da un numero gestibile di stringhe di bit campionate. Come mostrato di seguito, questa sparsità non è automatica; è una proprietà della base a particella singola in cui il problema è espresso. Scegliere correttamente la base è il punto cruciale di questa costruzione.

1. Costruire l'Hamiltoniano SIAM come operatore fermionico

Il modello di Anderson a singola impurezza pone un orbitale di impurezza interagente (indice $0$) a contatto con un bagno non interagente di $L$ siti disposti a catena. Il suo Hamiltoniano si divide in una parte a un corpo $H_1$ e una parte a due corpi $H_2$,

$$H = H_1 + H_2,$$

dove la parte a un corpo raccoglie l'energia on-site dell'impurezza (il potenziale chimico $\mu$), l'ibridazione impurezza-bagno $V$ e l'hopping tra primi vicini del bagno $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),$$

e la parte a due corpi contiene la repulsione coulombiana on-site $U$ sulla sola impurezza,

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

Questa guida lavora nel punto a simmetria particella-lacuna $\mu = -U/2$. Queste due parti si mappano direttamente sugli integrali a un corpo h1e (hopping del bagno, ibridazione dell'impurezza, energia on-site dell'impurezza) e sul singolo integrale a due corpi h2e (il termine coulombiano dell'impurezza). Assemblali prima nella base delle posizioni, in cui il modello è definito in modo naturale.

>>> 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. Passare alla base degli impulsi

La parte di bagno del SIAM è una catena di fermioni liberi, quindi si diagonalizza con un cambio alla base degli impulsi. Questo è il passaggio decisivo per SKQD. Poiché il bagno è invariante per traslazione, lo stato fondamentale non è sparso nella base delle posizioni, ma lo è nella base degli impulsi, dove l'Hamiltoniano del bagno è diagonale. La rotazione orbitale a particella singola che effettua questo cambio è algebra lineare specifica del problema. Diagonalizza la matrice di hopping del bagno lasciando intatto l'orbitale dell'impurezza, e poi sposta l'impurezza in un sito centrale:

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

A differenza di un flusso di lavoro variazionale, qui si ruotano gli integrali dell'Hamiltoniano nella base degli impulsi, in modo che l'intero circuito (Hamiltoniano e stato insieme) viva nella base in cui lo stato fondamentale è sparso e le stringhe di bit campionate sono quindi informative. La rotazione degli integrali a uno e a due corpi è un cambio di base standard dei tensori di struttura elettronica:

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

Poiché la rotazione lascia fisso l'orbitale dell'impurezza (agisce solo sul bagno), l'interazione on-site resta invariata. h2e_momentum rimane un singolo termine numero-numero, proprio come h2e. L'unica differenza è che l'impurezza è stata spostata su un indice centrale. Questo è visibile direttamente nei tensori degli integrali, ciascuno dei quali ha una sola voce non nulla:

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

Questa località è ciò che mantiene poco costosa la sintesi dell'interazione più avanti.

Assembla gli integrali a un corpo in un FermionOperator con il costruttore per integrali elettronici from_1body_tril_spin_sym(), che li richiede impacchettati in ordine triangolare inferiore. La parte a due corpi è un singolo termine numero-numero, quindi, invece di passare attraverso il meccanismo degli integrali a due corpi, costruiscila direttamente con cre()/ann(): l'operatore coulombiano on-site $U \, n_{p\uparrow} n_{p\downarrow}$ per il modo di impurezza $p$ (il suo modo $\uparrow$ all'indice p, il suo modo $\downarrow$ all'indice p + norb nella disposizione a blocchi di spin). Entrambi i costruttori raddoppiano per lo spin i norb orbitali spaziali in 2 * norb modi fermionici.

Mantieni la parte a un corpo $H_1$ e quella a due corpi $H_2$ come operatori separati, anziché sommarli subito, in accordo con i termini $H_1$ e $H_2$ delle equazioni della sezione 1. Costruire ciascuna parte una sola volta (in entrambe le basi) le rende riutilizzabili più avanti. Sommate, danno l'Hamiltoniano completo per la diagonalizzazione esatta in entrambe le basi. Singolarmente, $H_2$ (base degli impulsi) è l'operatore che il passo di Trotter fa evolvere nel passaggio 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)
... )

Il cambio di base è una trasformazione di similitudine unitaria, quindi lascia invariato lo spettro. Conferma che l'operatore sia corretto confrontando il suo autovalore più basso nel settore $(\text{norb}, \text{nelec})$ con una diagonalizzazione esatta. FermionOperator espone un LinearOperator di SciPy (supportato da un kernel nativo di prodotto matrice-vettore di interazione di configurazione completa (FCI)), quindi può essere passato direttamente a 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

Il vantaggio del cambio di base è la sparsità. Per renderla concreta, conta quanti determinanti della base computazionale servono per catturare il 99% del peso dello stato fondamentale e confronta la base degli impulsi con quella delle posizioni. L'energia non dipende dalla base, ma il numero di determinanti sì:

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

I due operatori hanno la stessa energia, ma lo stato fondamentale nella base delle posizioni è distribuito su migliaia di determinanti, mentre la base degli impulsi lo concentra su circa 20, ed è per questo che il campionamento SKQD funziona nella base degli impulsi e non in quella delle posizioni.

3. Preparare lo stato di riferimento

Il riferimento SKQD $\lvert \psi_0 \rangle$ è la sovrapposizione di tutte le eccitazioni degli elettroni più vicini al livello di Fermi nei modi di impulso vuoti più prossimi. La pubblicazione di SKQD lo costruisce con una rete scritta a mano di XXPlusYYGate. Sembra fortemente correlato (ha centinaia di ampiezze non nulle nella base computazionale), ma una rete di rotazioni a due modi che conservano il numero è un'operazione a particella singola (gaussiana fermionica), quindi lo stato che produce è un singolo determinante di Slater. Questo significa che puoi prepararlo con un solo Gate PrepareSlaterDeterminant, data la rotazione a particella singola che la rete implementa.

Poiché ogni XXPlusYYGate agisce come una rotazione di Givens su una coppia di modi, quella rotazione può essere costruita direttamente nella rappresentazione a particella singola (una matrice di Givens per Gate) invece di emettere a mano i Gate a due qubit. I Gate formano una struttura a mattoni di rotazioni tra primi vicini che si espande dal livello di Fermi a strati successivi, distribuendo ogni elettrone vicino al livello di Fermi sui modi vuoti immediatamente sopra di esso:

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

Nota

La rotazione di cui sopra è espressa nella base degli impulsi, coerentemente con l'Hamiltoniano. Gli elettroni vicini al livello di Fermi sono i modi di impulso occupati a più alta energia, e la rete li eccita nei modi vuoti a più bassa energia. Per questo il riferimento deve essere costruito dopo il cambio di base del passaggio precedente, non prima.

La preparazione del determinante di Slater viene eseguita per ciascun settore di spin, poiché ogni settore ha una propria forma (elettroni, orbitali). Un singolo Gate PrepareSlaterDeterminant (un'occupazione, con i nocc modi più bassi riempiti, insieme alla rotazione a particella singola) prepara un settore; lo stesso Gate viene applicato una volta ai modi alpha 0..norb e una volta ai modi beta norb..2*norb, poiché qui entrambi i settori di spin condividono la stessa occupazione di riferimento e la stessa rotazione:

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

Questo Gate dichiarativo sostituisce la rete di XXPlusYYGate scritta a mano; il Transpiler genera la sintesi delle rotazioni di Givens.

4. Assemblare i circuiti di Krylov

Ogni circuito di Krylov è la preparazione dello stato di riferimento seguita dall'evoluzione sotto $e^{-i k \, \Delta t \, H}$. Dividi l'Hamiltoniano nelle sue parti a un corpo e a due corpi, $H = H_1 + H_2$, e applica la trotterizzazione dell'evoluzione su questa suddivisione. Questa è la scelta decisiva per la profondità del circuito. L'evoluzione a un corpo $e^{-i t H_1}$ è un'operazione di fermioni liberi (gaussiana fermionica), quindi è una rotazione di base a particella singola (un OrbitalRotation con unitaria $e^{-i t \, \mathbf{h}_1}$) anziché qualcosa che deve essere trotterizzato in Pauli termine per termine. L'evoluzione a due corpi $e^{-i t H_2}$ agisce solo sull'unico termine di interazione on-site. Costruire il passo in questo modo permette al Transpiler di emettere una struttura a mattoni poco profonda di rotazioni di Givens più un singolo Gate di interazione a due qubit, la struttura poco profonda che la pubblicazione di SKQD costruisce a mano.

Abbiamo già la parte a due corpi. È l'operatore two_body_momentum costruito nel passaggio 2 (il singolo termine on-site). L'evoluzione a un corpo $e^{-i t H_1}$ è indipendente dallo spin, quindi è un singolo OrbitalRotation per settore con unitaria $e^{-i t \, \mathbf{h}_1}$, applicato una volta a ciascun settore di spin, la stessa struttura parallela della preparazione dello stato di riferimento.

Ogni circuito di Krylov evolve per un tempo totale $k \, \Delta t$. Ciò si realizza con un prodotto di Trotter del secondo ordine di $k$ passi di ampiezza $\Delta t$, in modo che l'errore per passo resti fisso al crescere della dimensione di Krylov. Ogni passo inserisce una rotazione completa a un corpo (una per settore di spin) tra due mezzi passi dell'evoluzione a due corpi (poco costosa, a singolo termine). La dimensione di Krylov $D$ è il numero di tali circuiti, con dim che va da 0 (solo stato di riferimento) a $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

Nota

I circuiti costruiti qui non contengono ancora misure, perché gli stessi circuiti fermionici guidano anche la simulazione esatta nel passaggio 6, e una misura non è un'operazione unitaria. Il passaggio 5 le aggiunge con FermionicCircuit.measure_all(), che legge il modo i nel bit classico i.

Preferisci misurare il circuito fermionico piuttosto che quello transpilato. Una misura fermionica indica il modo che legge, quindi layout e instradamento la portano con sé insieme al qubit su cui è finito quel modo, e la stringa di bit campionata è indicizzata per modo. Misurando invece il circuito transpilato, devi annullare tu stesso final_index_layout(), e se te ne dimentichi restituisce silenziosamente risultati permutati. Questo si basa sulla codifica di Jordan-Wigner in base di occupazione usata dal preset, in cui ogni modo è portato da un qubit.

5. Transpilare in circuiti di qubit

La transpilazione di un FermionicCircuit mappa i suoi Gate fermionici in Gate su qubit. Questo usa la mappatura di Jordan-Wigner jordan_wigner(). I Gate OrbitalRotation vengono sintetizzati in una struttura a mattoni di rotazioni di Givens, e il singolo Evolution a termine singolo si riduce a un singolo Gate di interazione a due qubit.

Poiché ogni Gate fermionico qui ha una sintesi predefinita, puoi usare direttamente il già pronto generate_preset_jw_pass_manager() invece di assemblare a mano le fasi. Eventuali argomenti per parola chiave vengono inoltrati alla fase dei qubit; passa optimization_level=0 per un'immagine fedele e non ottimizzata della profondità sintetizzata. Generare l'intera famiglia di circuiti di Krylov è quindi un ciclo sulla dimensione di Krylov. Conserva i FermionicCircuit non transpilati accanto ai circuiti di qubit transpilati. Quelli fermionici guidano la simulazione esatta (a vettore di stato) nel passaggio 6, quindi devono restare privi di misure; measure_all(inplace=False) restituisce una copia con misure da transpilare e lascia intatto l'originale, e quella copia è ciò che un Backend eseguirebbe.

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

La dimensione di Krylov emerge direttamente nella profondità a due qubit dei circuiti. La preparazione dello stato di riferimento è un costo fisso, e ogni ulteriore potenza di Krylov aggiunge un altro passo di Trotter del secondo ordine dell'operatore di evoluzione temporale.

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

Nota

Dividere l'Hamiltoniano mantiene contenute queste profondità. L'evoluzione a un corpo è una vera rotazione orbitale, quindi si sintetizza in una struttura a mattoni di XXPlusYYGate su qubit adiacenti. Ogni mezzo passo dell'interazione on-site è un singolo RZZGate che accoppia i due modi di spin dell'impurezza; il termine poco costoso che il prodotto del secondo ordine inserisce attorno alla rotazione. Passare l'intero Hamiltoniano a un unico Gate Evolution comporterebbe invece la trotterizzazione di Pauli della parte a un corpo termine per termine, propagando lunghe stringhe Z di Jordan-Wigner attraverso gli accoppiamenti impurezza-bagno (a lungo raggio) e aumentando la profondità di un ordine di grandezza.

Il primo circuito prepara solo il determinante di riferimento, mentre l'ultimo contiene il maggior numero di passi di Trotter (quattro in questo caso):

>>> 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. Campionare le stringhe di bit dai circuiti di Krylov

Su hardware, eseguiresti i circuiti transpilati e misureresti. Questa è la versione senza rumore (a vettore di stato) di quel passaggio, che mostra le stringhe di bit che SKQD raccoglierebbe. Simula ogni FermionicCircuit non transpilato con ffsim.apply_unitary() e campiona il vettore di stato risultante con ffsim.sample_state_vector().

Un Gate PrepareSlaterDeterminant in simulazione esegue prima la validazione e poi la rotazione (verifica che lo stato in ingresso occupi il suo determinante di riferimento, poi applica la rotazione). Pertanto, l'input è il semplice determinante di occupazione: i nocc modi più bassi riempiti per spin, senza rotazione. Il Gate applica da sé 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)

Raggruppando i campioni di tutti e cinque i circuiti di Krylov, il supporto è minuscolo: poche centinaia di stringhe di bit distinte su $\binom{8}{4}^2 = 4900$ determinanti nel settore (4, 4). Questa concentrazione, conseguenza diretta della sparsità nella base degli impulsi stabilita nel passaggio 2, rende trattabile la successiva diagonalizzazione classica.

Tracciare i conteggi con qiskit.visualization.plot_histogram() rende visibile la sparsità. Poche configurazioni dominano, guidate dal determinante di riferimento.

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

7. Diagonalizza nel sottospazio campionato

Le bitstring campionate sono l'input della metà classica di SKQD. L'Hamiltoniana viene proiettata sul sottospazio generato dalle configurazioni campionate e diagonalizzata al suo interno. Questo viene affidato al pacchetto qiskit-addon-sqd, la cui funzione diagonalize_fermionic_hamiltonian() esegue l'intero ciclo di diagonalizzazione quantistica basata su campioni. Costruisce un sottospazio a partire da batch delle configurazioni campionate, diagonalizza al suo interno l'Hamiltoniana e raffina iterativamente il sottospazio usando il recupero delle configurazioni: inverte le occupazioni per riparare le configurazioni che violano la simmetria nota del numero di particelle. Il recupero delle configurazioni è pensato per annullare il danno del rumore dell'hardware, che corrompe le bitstring campionate portandole nel settore sbagliato del numero di particelle; qui il campionamento è privo di rumore, quindi su questo fronte non c'è nulla da riparare.

Il solver consuma i conteggi campionati come BitArray, costruito direttamente a partire dai counts aggregati del passo 6. Passa num_bits in modo esplicito: senza di esso, la larghezza viene dedotta dalla bitstring campionata più larga, che è più stretta di 2 * norb ogni volta che nessun campione occupa l'ultimo modo, e il solver rifiuta quindi l'array. Il solver diagonalizza nella stessa base degli impulsi usata dal resto della costruzione, quindi prende direttamente gli integrali h1e_momentum e h2e_momentum del passo 2, senza alcun riordinamento. L'energia recuperata coincide con il riferimento esatto entro molto meno di un milli-Hartree, partendo da poche centinaia di configurazioni campionate. Poiché l'energia SQD è un limite superiore variazionale dello stato fondamentale reale, si trova sempre al livello di reference_energy o al di sopra (ancora in ambito dal passo 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)

Visualizza l'esecuzione: la convergenza dell'energia nel corso delle iterazioni e l'occupazione media di ciascun orbitale spaziale nello stato fondamentale recuperato. L'energia per iterazione è la più bassa tra i batch di quell'iterazione, e l'occupazione somma le orbital_occupancies sui due settori di spin. Anche senza rumore da correggere, l'energia diminuisce da un'iterazione all'altra. Ogni iterazione diagonalizza all'interno di num_batches batch di sole samples_per_batch configurazioni sottocampionate dal pool, e il raffinamento iterativo dovuto al riporto delle bitstring migliora progressivamente quali configurazioni finiscono in quei batch, così l'energia recuperata migliora:

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

Questo chiude il ciclo SKQD. I circuiti in base degli impulsi dei passi 1–5 producono un campione sparso (passo 6), e la diagonalizzazione basata su campioni qui sopra trasforma quel campione di nuovo in un'energia dello stato fondamentale.

Prossimi passi

Questa guida ha eseguito l'intera pipeline SKQD su un simulatore di vettore di stato privo di rumore. Su hardware, i circuits transpilati verrebbero eseguiti e misurati al posto del campionamento del vettore di stato nel passo 6, e i conteggi risultanti alimenterebbero invariati la diagonalizzazione del passo 7. Consulta la documentazione di Qiskit per un aiuto nell'esecuzione dei circuiti e i tutorial dell'addon SQD per informazioni sul post-processing di diagonalizzazione del sottospazio, incluso il suo comportamento su campioni rumorosi, dove il recupero delle configurazioni svolge la maggior parte del lavoro.

Leggi la guida sulla relazione con ffsim per capire perché ffsim.apply_unitary() e ffsim.sample_state_vector() funzionano nativamente sui circuiti fermionici di questo pacchetto, e perché il settore a numero di particelle fisso usato nel passo 6 rende efficiente in termini di campioni la sparsità in base degli impulsi stabilita nel passo 2.