Generování obvodů časového vývoje Krylov (SKQD)

Důležité

Koncepty v tomto průvodci jsou v současnosti dostupné pouze v Python API. Ekvivalentní funkcionalita bude v C API zpřístupněna v budoucím vydání.

Sample-based Krylov Quantum Diagonalization (SKQD) je kvantově centrická varianta diagonalizace v podprostoru. Místo vzorkování jediného variačního ansatzu (jako v prosté sample-based quantum diagonalization (SQD)) nebo souboru randomizovaných obvodů časového vývoje (jako v SqDRIFT) vzorkuje Krylovovu bázi: rodinu stavů

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

získaných vývojem referenčního stavu $\lvert \psi_0 \rangle$ po rostoucí násobky pevného časového kroku $\Delta t$. Bitové řetězce vzorkované z těchto $D$ obvodů vyplňují Krylovův podprostor, v němž se hamiltonián následně klasicky diagonalizuje.

Tento průvodce reprodukuje krok konstrukce obvodu algoritmu SKQD pro Andersonův model jedné příměsi (SIAM), tedy magnetickou příměs vázanou na neinteragující lázeň, podle publikace SKQD.

Poznámka

Krok vzorkování SKQD funguje pouze tehdy, je-li základní stav v počítačové bázi řídký: klasická diagonalizace ho musí být schopna rekonstruovat z zvládnutelného počtu navzorkovaných bitových řetězců. Jak je ukázáno níže, tato řídkost není automatická, je vlastností jednočásticové báze, ve které je problém vyjádřen. Správná volba báze je jádrem této konstrukce.

1. Sestavení hamiltoniánu SIAM jako fermionového operátoru

Andersonův model jedné příměsi umisťuje jeden interagující orbital příměsi (index $0$) do kontaktu s neinteragující lázní o $L$ místech uspořádaných do řetězce. Jeho hamiltonián se rozpadá na jednočásticovou část $H_1$ a dvoučásticovou část $H_2$,

$$H = H_1 + H_2,$$

přičemž jednočásticová část zahrnuje místní energii příměsi (chemický potenciál $\mu$), hybridizaci příměsi s lázní $V$ a poskakování v lázni mezi nejbližšími sousedy $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),$$

a dvoučásticová část nese místní coulombovské odpuzování $U$ pouze na příměsi,

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

Tento průvodce pracuje v bodě s částicově-děrovou symetrií $\mu = -U/2$. Tyto dvě části se přímo mapují na jednočásticové integrály h1e (poskakování v lázni, hybridizace příměsi, místní energie příměsi) a jediný dvoučásticový integrál h2e (coulombovský člen příměsi). Nejprve je sestav v poziční bázi, ve které je model přirozeně definován.

>>> 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. Přechod do hybnostní báze

Lázeň v SIAM je řetězec volných fermionů, takže se diagonalizuje přechodem do hybnostní báze. To je pro SKQD rozhodující krok. Protože je lázeň translačně invariantní, základní stav není v poziční bázi řídký, ale v hybnostní bázi, kde je hamiltonián lázně diagonální, řídký je. Jednočásticová orbitální rotace, která tento přechod provádí, je lineární algebra specifická pro problém. Diagonalizuje matici poskakování v lázni, přičemž orbital příměsi nechává nedotčený, a poté přesune příměs na centrální místo:

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

Na rozdíl od variačního pracovního postupu se zde rotují integrály hamiltoniánu do hybnostní báze, takže celý obvod (hamiltonián i stav) žije v bázi, kde je základní stav řídký, a navzorkované bitové řetězce jsou proto informativní. Rotace jedno- a dvoučásticových integrálů je standardní změna báze tenzorů elektronové struktury:

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

Protože rotace nechává orbital příměsi pevný (působí pouze na lázeň), místní interakce se nemění. h2e_momentum zůstává jediným členem typu číslo-číslo, stejně jako h2e. Jedinou změnou je, že příměs byla přesunuta na centrální index. To je přímo vidět na integrálních tenzorech, z nichž každý nese jednu nenulovou položku:

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

Tato lokálnost zajišťuje, že interakci bude později levné syntetizovat.

Jednočásticové integrály sestav do FermionOperator pomocí konstruktoru elektronových integrálů from_1body_tril_spin_sym(), který očekává, že budou zabaleny v pořadí dolní trojúhelníkové matice. Dvoučásticová část je jediný člen typu číslo-číslo, takže ji místo vedení přes mechanismus dvoučásticových integrálů sestav přímo pomocí cre()/ann(): místní coulombovský operátor $U \, n_{p\uparrow} n_{p\downarrow}$ pro mód příměsi $p$ (jeho mód $\uparrow$ na indexu p, jeho mód $\downarrow$ na indexu p + norb v blokovém uspořádání spinů). Oba konstruktory zdvojnásobují norb prostorových orbitalů o spin na 2 * norb fermionových módů.

Jednočásticovou část $H_1$ a dvoučásticovou část $H_2$ ponech jako samostatné operátory a nesčítej je předčasně, v souladu se členy $H_1$ a $H_2$ z rovnic v části 1. Jednorázové sestavení každé části (v obou bázích) umožňuje jejich opakované použití níže. Po sečtení dávají úplný hamiltonián pro přesnou diagonalizaci v obou bázích. Samostatně je $H_2$ (hybnostní báze) operátor, který Trotterův krok vyvíjí v kroku 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)
... )

Změna báze je unitární podobnostní transformace, takže spektrum zůstává nedotčeno. Správnost operátoru ověř porovnáním jeho nejnižšího vlastního čísla v sektoru $(\text{norb}, \text{nelec})$ s přesnou diagonalizací. FermionOperator zpřístupňuje SciPy LinearOperator (podložený nativním jádrem násobení matice vektorem pro full configuration interaction (FCI)), takže ho lze předat přímo do 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

Přínosem změny báze je řídkost. Pro konkrétnost spočítej, kolik determinantů výpočetní báze je potřeba k zachycení 99 % váhy základního stavu, a porovnej hybnostní bázi s poziční. Energie na bázi nezávisí, ale počet determinantů ano:

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

Oba operátory sdílejí stejnou energii, ale základní stav v poziční bázi je rozprostřen přes tisíce determinantů, zatímco hybnostní báze ho soustředí do zhruba 20, a proto vzorkování SKQD funguje v hybnostní bázi a ne v poziční.

3. Příprava referenčního stavu

Referenční stav SKQD $\lvert \psi_0 \rangle$ je superpozice všech excitací elektronů nejblíže Fermiho hladině do nejbližších prázdných hybnostních módů. Publikace SKQD ho sestavuje ručně psanou sítí XXPlusYYGate. Vypadá silně korelovaně (má stovky nenulových amplitud ve výpočetní bázi), ale síť rotací dvou módů zachovávajících počet částic je jednočásticová (fermionová gaussovská) operace, takže stav, který vytváří, je jediný Slaterův determinant. To znamená, že ho můžeš připravit jediným hradlem PrepareSlaterDeterminant, pokud znáš jednočásticovou rotaci, kterou síť implementuje.

Protože každé XXPlusYYGate působí jako Givensova rotace na dvojici módů, lze tuto rotaci sestavit přímo v jednočásticovém obrazu (jedna Givensova matice na hradlo), místo aby se dvoukvbitová hradla vypisovala ručně. Hradla tvoří cihlový vzor rotací mezi nejbližšími sousedy, který se rozvětvuje od Fermiho hladiny v rozšiřujících se vrstvách a rozprostírá každý elektron poblíž Fermiho hladiny přes prázdné módy těsně nad ní:

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

Poznámka

Výše uvedená rotace je vyjádřena v hybnostní bázi, shodně s hamiltoniánem. Elektrony poblíž Fermiho hladiny jsou obsazené hybnostní módy s nejvyšší energií a síť je excituje do prázdných módů s nejnižší energií. Proto musí být reference sestavena po změně báze z předchozího kroku, ne před ní.

Příprava Slaterova determinantu se provádí pro každý spinový sektor zvlášť, protože každý sektor má svůj vlastní tvar (elektrony, orbitaly). Jediné hradlo PrepareSlaterDeterminant (obsazení, tedy zaplněných nejnižších nocc módů, spolu s jednočásticovou rotací) připraví jeden sektor; stejné hradlo se použije jednou na módy alfa 0..norb a jednou na módy beta norb..2*norb, protože oba spinové sektory zde sdílejí stejné referenční obsazení i rotaci:

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

Toto deklarativní hradlo nahrazuje ručně psanou síť XXPlusYYGate; syntézu Givensových rotací generuje transpiler.

4. Sestavení Krylovových obvodů

Každý Krylovův obvod je příprava referenčního stavu následovaná vývojem pod $e^{-i k \, \Delta t \, H}$. Hamiltonián rozděl na jednočásticovou a dvoučásticovou část, $H = H_1 + H_2$, a vývoj Trotterizuj podle tohoto rozdělení. To je rozhodující volba pro hloubku obvodu. Jednočásticový vývoj $e^{-i t H_1}$ je operace volných fermionů (fermionová gaussovská), takže jde o jednočásticovou rotaci báze (OrbitalRotation s unitární maticí $e^{-i t \, \mathbf{h}_1}$), a nikoli o něco, co by se muselo Pauliho-Trotterizovat člen po členu. Dvoučásticový vývoj $e^{-i t H_2}$ působí pouze na jediný člen místní interakce. Takto sestavený krok umožňuje transpileru vytvořit mělký cihlový vzor Givensových rotací plus jediné dvoukvbitové interakční hradlo, tedy mělkou strukturu, kterou publikace SKQD sestavuje ručně.

Dvoučásticovou část už máme. Je to operátor two_body_momentum sestavený v kroku 2 (jediný místní člen). Jednočásticový vývoj $e^{-i t H_1}$ nezávisí na spinu, takže jde o jedinou OrbitalRotation pro každý sektor s unitární maticí $e^{-i t \, \mathbf{h}_1}$, použitou jednou na každý spinový sektor, tedy se stejnou paralelní strukturou jako příprava referenčního stavu.

Každý Krylovův obvod se vyvíjí po celkovou dobu $k \, \Delta t$. To je realizováno Trotterovým součinem druhého řádu z $k$ kroků velikosti $\Delta t$, takže chyba na krok zůstává pevná, i když Krylovova dimenze roste. Každý krok vkládá úplnou jednočásticovou rotaci (jednu na spinový sektor) mezi dva poloviční kroky (levného, jednočlenného) dvoučásticového vývoje. Krylovova dimenze $D$ je počet takových obvodů, přičemž dim probíhá od 0 (pouze referenční stav) do $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

Poznámka

Zde sestavené obvody zatím neobsahují žádná měření, protože tytéž fermionové obvody pohánějí i přesnou simulaci v kroku 6 a měření není unitární operace. Krok 5 je přidává pomocí FermionicCircuit.measure_all(), které čte mód i do klasického bitu i.

Raději měř fermionový obvod než transpilovaný. Fermionové měření pojmenovává mód, který čte, takže layout a routing ho přenesou s qubitem, na kterém mód skončil, a navzorkovaný bitový řetězec je indexován podle módu. Měření transpilovaného obvodu naopak ponechává na tobě, abys sám vrátil final_index_layout(), a pokud na to zapomeneš, tiše vrací permutované výsledky. To spoléhá na Jordanovo-Wignerovo kódování v bázi obsazení, které preset používá, kde je každý mód nesen jedním qubitem.

5. Transpilace na qubitové obvody

Transpilace FermionicCircuit mapuje její fermionová hradla na qubitová hradla. Používá se Jordanovo-Wignerovo mapování jordan_wigner(). Hradla OrbitalRotation se syntetizují do cihlového vzoru Givensových rotací a jednočlenné Evolution se snižuje na jediné dvoukvbitové interakční hradlo.

Protože každé zdejší fermionové hradlo má výchozí syntézu, můžeš přímo použít hotový generate_preset_jw_pass_manager() místo ručního skládání jednotlivých fází. Všechny pojmenované argumenty se předávají qubitové fázi; předej optimization_level=0 pro věrný, neoptimalizovaný obraz syntetizované hloubky. Vygenerování celé rodiny Krylovových obvodů je pak smyčka přes Krylovovu dimenzi. Netranspilované FermionicCircuit si ponech vedle transpilovaných qubitových obvodů. Fermionové pohánějí přesnou (statevector) simulaci v kroku 6, takže musí zůstat bez měření; measure_all(inplace=False) vrací změřenou kopii k transpilaci a původní nechává nedotčený a právě tato kopie je to, co by backend spustil.

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

Krylovova dimenze se přímo projevuje v dvoukvbitové hloubce obvodů. Příprava referenčního stavu je pevná cena a každá další Krylovova mocnina přidává další Trotterův krok druhého řádu operátoru časového vývoje.

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

Poznámka

Rozdělení hamiltoniánu udržuje tyto hloubky skromné. Jednočásticový vývoj je skutečná orbitální rotace, takže se syntetizuje do cihlového vzoru XXPlusYYGate na sousedních qubitech. Každý poloviční krok místní interakce je jediné RZZGate, které spojuje dva spinové módy příměsi; je to levný člen, který součin druhého řádu vkládá kolem rotace. Předání celého hamiltoniánu jedinému hradlu Evolution by naopak Pauliho-Trotterizovalo jednočásticovou část člen po členu, šířilo by dlouhé Jordanovy-Wignerovy řetězce Z přes (dalekodosahové) vazby příměsi s lázní a zvýšilo by hloubku o řád.

První obvod připravuje pouze referenční determinant, zatímco poslední nese nejvíce Trotterových kroků (v tomto případě čtyři):

>>> 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. Vzorkování bitových řetězců z Krylovových obvodů

Na hardwaru bys transpilované obvody spustil a změřil. Toto je bezšumová (statevector) verze onoho kroku, která ukazuje bitové řetězce, jež by SKQD shromáždilo. Každý netranspilovaný FermionicCircuit simuluj pomocí ffsim.apply_unitary() a z výsledného stavového vektoru vzorkuj pomocí ffsim.sample_state_vector().

Hradlo PrepareSlaterDeterminant se v simulaci chová jako ověř-a-rotuj (zkontroluje, že vstupní stav obsazuje svůj referenční determinant, a poté aplikuje rotaci). Proto je vstupem prostý determinant obsazení: nejnižší nocc módů zaplněných pro každý spin, bez rotace. Hradlo samo aplikuje 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)

Při sloučení vzorků ze všech pěti Krylovových obvodů je nosič drobný: několik set odlišných bitových řetězců z $\binom{8}{4}^2 = 4900$ determinantů v sektoru (4, 4). Tato koncentrace, přímý důsledek řídkosti v hybnostní bázi zjištěné v kroku 2, činí následnou klasickou diagonalizaci zvládnutelnou.

Vykreslení počtů pomocí qiskit.visualization.plot_histogram() zviditelní řídkost. Dominuje několik konfigurací v čele s referenčním determinantem.

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

7. Diagonalizace ve vzorkovaném podprostoru

Vzorkované bitové řetězce jsou vstupem klasické poloviny SKQD. Hamiltonián se promítne do podprostoru, který vzorkované konfigurace generují, a v něm se diagonalizuje. Tento krok se předává balíčku qiskit-addon-sqd, jehož funkce diagonalize_fermionic_hamiltonian() spouští celou smyčku vzorkové kvantové diagonalizace. Sestaví podprostor z dávek vzorkovaných konfigurací, diagonalizuje v něm hamiltonián a podprostor iterativně zpřesňuje pomocí obnovy konfigurací: přepínáním obsazení opravuje konfigurace, které porušují známou symetrii zachování počtu částic. Obnova konfigurací je navržena tak, aby napravila škody způsobené šumem hardwaru, který poškozuje vzorkované bitové řetězce do nesprávného sektoru počtu částic; zdejší vzorkování je bezšumové, takže v tomto ohledu není co opravovat.

Řešič spotřebovává vzorkované počty jako BitArray, vytvořené přímo ze sloučených counts z kroku 6. Zadej num_bits explicitně: bez něj se šířka odvodí z nejširšího vzorkovaného bitového řetězce, který je užší než 2 * norb, kdykoli žádný vzorek neobsazuje poslední mód, a řešič pak pole odmítne. Řešič diagonalizuje ve stejné hybnostní bázi, jakou používá zbytek konstrukce, takže přímo přebírá integrály h1e_momentum a h2e_momentum z kroku 2, bez jakéhokoli přeskupování. Získaná energie odpovídá přesné reference s odchylkou hluboko pod miliHartree, přičemž stačí jen několik set vzorkovaných konfigurací. Protože je energie SQD variační horní mez skutečného základního stavu, vždy leží na hodnotě reference_energy (stále platné z kroku 2) nebo nad ní:

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

Vizualizuj průběh: konvergenci energie napříč iteracemi a průměrné obsazení každého prostorového orbitalu v získaném základním stavu. Energie v každé iteraci je nejnižší napříč dávkami dané iterace a obsazení sčítá orbital_occupancies přes dva spinové sektory. I bez šumu, který by se opravoval, energie klesá z iterace na iteraci. Každá iterace diagonalizuje v rámci num_batches dávek, každá jen s samples_per_batch konfiguracemi vybranými z množiny, a iterativní zpřesňování díky přenosu bitových řetězců neustále zlepšuje, které konfigurace se do těchto dávek dostanou, takže se získaná energie zlepšuje:

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

Tím se smyčka SKQD uzavírá. Obvody v hybnostní bázi z kroků 1–5 vytvoří řídký vzorek (krok 6) a vzorková diagonalizace výše tento vzorek převede zpět na energii základního stavu.

Další kroky

Tento návod provedl celý postup SKQD na bezšumovém simulátoru stavového vektoru. Na hardwaru by se transpilované circuits spustily a změřily namísto vzorkování stavového vektoru v kroku 6 a výsledné počty by vstupovaly do diagonalizace v kroku 7 beze změny. Pomoc se spouštěním obvodů najdeš v dokumentaci Qiskit a informace o následném zpracování diagonalizací podprostoru, včetně jeho chování na zašuměných vzorcích, kde obnova konfigurací odvede nejvíc práce, najdeš v tutoriálech doplňku SQD.

Přečti si návod k vztahu s ffsim, abys pochopil, proč ffsim.apply_unitary() a ffsim.sample_state_vector() fungují nativně na fermionových obvodech tohoto balíčku a proč sektor s pevným počtem částic použitý v kroku 6 činí řídkost v hybnostní bázi zjištěnou v kroku 2 efektivní z hlediska vzorků.