Generowanie obwodów ewolucji czasowej Kryłowa (SKQD)

Ważne

Pojęcia z tego przewodnika są obecnie dostępne tylko w API Pythona. Równoważna funkcjonalność zostanie udostępniona w API C w przyszłym wydaniu.

Sample-based Krylov Quantum Diagonalization (SKQD) to kwantowo-centryczny wariant diagonalizacji podprzestrzeni. Zamiast próbkować pojedynczy wariacyjny ansatz (jak w zwykłej próbkowej diagonalizacji kwantowej (SQD)) lub zespół losowych obwodów ewolucji czasowej (jak w SqDRIFT), próbkuje bazę Kryłowa: rodzinę stanów

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

otrzymanych przez ewolucję stanu referencyjnego $\lvert \psi_0 \rangle$ przez rosnące wielokrotności stałego kroku czasowego $\Delta t$. Ciągi bitów próbkowane z tych $D$ obwodów zapełniają podprzestrzeń Kryłowa, w której hamiltonian jest następnie diagonalizowany klasycznie.

Ten przewodnik odtwarza krok konstrukcji obwodu algorytmu SKQD dla modelu Andersona z pojedynczą domieszką (SIAM), czyli domieszki magnetycznej sprzężonej z nieoddziałującą kąpielą, zgodnie z publikacją SKQD.

Uwaga

Krok próbkowania SKQD działa tylko wtedy, gdy stan podstawowy jest rzadki w bazie obliczeniowej: klasyczna diagonalizacja musi być w stanie odtworzyć go z rozsądnej liczby próbkowanych ciągów bitów. Jak pokazano poniżej, ta rzadkość nie jest automatyczna; jest właściwością bazy jednocząstkowej, w której wyrażono problem. Dobranie właściwej bazy jest sednem tej konstrukcji.

1. Zbuduj hamiltonian SIAM jako operator fermionowy

Model Andersona z pojedynczą domieszką umieszcza jeden oddziałujący orbital domieszki (indeks $0$) w kontakcie z nieoddziałującą kąpielą złożoną z $L$ węzłów ułożonych w łańcuch. Jego hamiltonian dzieli się na część jednociałową $H_1$ i część dwuciałową $H_2$,

$$H = H_1 + H_2,$$

przy czym część jednociałowa zbiera energię lokalną domieszki (potencjał chemiczny $\mu$), hybrydyzację domieszki z kąpielą $V$ oraz przeskoki kąpieli między najbliższymi sąsiadami $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 część dwuciałowa niesie lokalne odpychanie kulombowskie $U$ wyłącznie na domieszce,

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

Ten przewodnik pracuje w punkcie symetrycznym względem cząstek i dziur $\mu = -U/2$. Te dwie części odwzorowują się bezpośrednio na całki jednociałowe h1e (przeskoki kąpieli, hybrydyzacja domieszki, energia lokalna domieszki) oraz pojedynczą całkę dwuciałową h2e (człon kulombowski domieszki). Najpierw złóż je w bazie położeniowej, w której model jest definiowany w naturalny sposób.

>>> 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. Przejdź do bazy pędowej

Część kąpielowa SIAM to łańcuch swobodnych fermionów, więc jest diagonalizowana przez zmianę na bazę pędową. To decydujący krok dla SKQD. Ponieważ kąpiel jest niezmiennicza względem translacji, stan podstawowy nie jest rzadki w bazie położeniowej, ale jest rzadki w bazie pędowej, w której hamiltonian kąpieli jest diagonalny. Jednocząstkowa rotacja orbitali, która dokonuje tej zmiany, to algebra liniowa specyficzna dla problemu. Diagonalizuje ona macierz przeskoków kąpieli, pozostawiając orbital domieszki nietknięty, a następnie przenosi domieszkę do węzła centralnego:

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

W przeciwieństwie do przepływu wariacyjnego, tutaj obracamy całki hamiltonianu do bazy pędowej, tak aby cały obwód (zarówno hamiltonian, jak i stan) żył w bazie, w której stan podstawowy jest rzadki, a próbkowane ciągi bitów są dzięki temu informatywne. Obrót całek jedno- i dwuciałowych to standardowa zmiana bazy tensorów struktury elektronowej:

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

Ponieważ rotacja pozostawia orbital domieszki niezmieniony (działa tylko na kąpiel), oddziaływanie lokalne się nie zmienia. h2e_momentum pozostaje pojedynczym członem liczba-liczba, tak jak h2e. Jedyną zmianą jest to, że domieszka została przeniesiona do centralnego indeksu. Widać to bezpośrednio w tensorach całek, z których każdy zawiera jeden niezerowy element:

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

Ta lokalność sprawia, że oddziaływanie jest później tanie w syntezie.

Złóż całki jednociałowe w FermionOperator za pomocą konstruktora całek elektronowych from_1body_tril_spin_sym(), który oczekuje ich spakowania w uporządkowaniu dolnotrójkątnym. Część dwuciałowa to pojedynczy człon liczba-liczba, więc zamiast przepuszczać ją przez mechanizm całek dwuciałowych, zbuduj ją bezpośrednio za pomocą cre()/ann(): lokalny operator kulombowski $U \, n_{p\uparrow} n_{p\downarrow}$ dla modu domieszki $p$ (jej mod $\uparrow$ pod indeksem p, jej mod $\downarrow$ pod indeksem p + norb w układzie blokowym względem spinu). Oba konstruktory podwajają spinowo norb orbitali przestrzennych do 2 * norb modów fermionowych.

Zachowaj część jednociałową $H_1$ i część dwuciałową $H_2$ jako osobne operatory, zamiast sumować je od razu, zgodnie ze składnikami $H_1$ i $H_2$ z równań w sekcji 1. Zbudowanie każdej części raz (w obu bazach) pozwala użyć ich ponownie poniżej. Zsumowane dają pełny hamiltonian do dokładnej diagonalizacji w obu bazach. Osobno $H_2$ (baza pędowa) jest operatorem, który krok Trottera ewoluuje w 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)
... )

Zmiana bazy jest unitarną transformacją podobieństwa, więc nie zmienia widma. Potwierdź poprawność operatora, porównując jego najniższą wartość własną w sektorze $(\text{norb}, \text{nelec})$ z dokładną diagonalizacją. FermionOperator udostępnia SciPy LinearOperator (oparty na natywnym jądrze mnożenia macierzy przez wektor pełnego oddziaływania konfiguracji (FCI)), więc można go przekazać wprost 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

Zyskiem ze zmiany bazy jest rzadkość. Aby to uzmysłowić, policz, ile wyznaczników w bazie obliczeniowej potrzeba do uchwycenia 99% wagi stanu podstawowego, i porównaj bazę pędową z położeniową. Energia nie zależy od bazy, ale liczba wyznaczników zależy:

>>> 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 operatory mają tę samą energię, ale stan podstawowy w bazie położeniowej jest rozłożony na tysiące wyznaczników, podczas gdy baza pędowa koncentruje go na około 20, dlatego próbkowanie SKQD działa w bazie pędowej, a nie w położeniowej.

3. Przygotuj stan referencyjny

Referencją SKQD $\lvert \psi_0 \rangle$ jest superpozycja wszystkich wzbudzeń elektronów najbliższych poziomowi Fermiego do najbliższych pustych modów pędowych. Publikacja SKQD buduje ją za pomocą ręcznie napisanej sieci XXPlusYYGate. Wygląda na silnie skorelowaną (ma setki niezerowych amplitud w bazie obliczeniowej), ale sieć zachowujących liczbę cząstek rotacji dwóch modów to operacja jednocząstkowa (fermionowa gaussowska), więc wytwarzany przez nią stan jest pojedynczym wyznacznikiem Slatera. Oznacza to, że możesz przygotować go jedną bramką PrepareSlaterDeterminant, znając jednocząstkową rotację, którą ta sieć implementuje.

Ponieważ każda XXPlusYYGate działa jak rotacja Givensa na parze modów, tę rotację można zbudować bezpośrednio w obrazie jednocząstkowym (jedna macierz Givensa na bramkę), zamiast ręcznie emitować bramki dwukubitowe. Bramki tworzą układ cegiełkowy rotacji między najbliższymi sąsiadami, który rozchodzi się od poziomu Fermiego w rozszerzających się warstwach, rozpraszając każdy elektron w pobliżu poziomu Fermiego na puste mody tuż nad nim:

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

Uwaga

Powyższa rotacja jest wyrażona w bazie pędowej, zgodnie z hamiltonianem. Elektrony w pobliżu poziomu Fermiego to zajęte mody pędowe o najwyższej energii, a sieć wzbudza je do pustych modów o najniższej energii. Dlatego referencję trzeba zbudować po zmianie bazy z poprzedniego kroku, a nie przed nią.

Przygotowanie wyznacznika Slatera wykonuje się osobno dla każdego sektora spinowego, ponieważ każdy sektor ma własny kształt (elektrony, orbitale). Pojedyncza bramka PrepareSlaterDeterminant (obsadzenie, z wypełnionymi najniższymi modami nocc, wraz z rotacją jednocząstkową) przygotowuje jeden sektor; ta sama bramka jest stosowana raz na modach alfa 0..norb i raz na modach beta norb..2*norb, ponieważ oba sektory spinowe mają tu to samo obsadzenie referencyjne i tę samą rotację:

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

Ta deklaratywna bramka zastępuje ręcznie kodowaną sieć XXPlusYYGate; transpiler generuje syntezę rotacji Givensa.

4. Złóż obwody Kryłowa

Każdy obwód Kryłowa to przygotowanie stanu referencyjnego, po którym następuje ewolucja pod działaniem $e^{-i k \, \Delta t \, H}$. Podziel hamiltonian na część jednociałową i dwuciałową, $H = H_1 + H_2$, i dokonaj trotteryzacji ewolucji względem tego podziału. To decydujący wybór dla głębokości obwodu. Ewolucja jednociałowa $e^{-i t H_1}$ jest operacją swobodnych fermionów (fermionową gaussowską), więc jest jednocząstkową rotacją bazy (OrbitalRotation z unitarną $e^{-i t \, \mathbf{h}_1}$), a nie czymś, co trzeba trotteryzować Pauli po składniku. Ewolucja dwuciałowa $e^{-i t H_2}$ działa tylko na pojedynczy człon oddziaływania lokalnego. Zbudowanie kroku w ten sposób pozwala transpilerowi wyemitować płytki układ cegiełkowy rotacji Givensa oraz pojedynczą dwukubitową bramkę oddziaływania, czyli płytką strukturę, którą publikacja SKQD buduje ręcznie.

Część dwuciałową już mamy. To operator two_body_momentum zbudowany w kroku 2 (pojedynczy człon lokalny). Ewolucja jednociałowa $e^{-i t H_1}$ nie zależy od spinu, więc jest pojedynczą rotacją OrbitalRotation na sektor z unitarną $e^{-i t \, \mathbf{h}_1}$, zastosowaną raz do każdego sektora spinowego, w tej samej równoległej strukturze co przygotowanie stanu referencyjnego.

Każdy obwód Kryłowa ewoluuje przez łączny czas $k \, \Delta t$. Realizuje się to drugorzędowym iloczynem Trottera z $k$ kroków o rozmiarze $\Delta t$, tak aby błąd na krok pozostawał stały wraz ze wzrostem wymiaru Kryłowa. Każdy krok umieszcza pełną rotację jednociałową (po jednej na sektor spinowy) między dwoma półkrokami (taniej, jednoczłonowej) ewolucji dwuciałowej. Wymiar Kryłowa $D$ to liczba takich obwodów, przy czym dim przebiega od 0 (tylko stan referencyjny) 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

Uwaga

Zbudowane tu obwody nie zawierają jeszcze pomiarów, ponieważ te same obwody fermionowe napędzają też dokładną symulację w kroku 6, a pomiar nie jest operacją unitarną. Krok 5 dodaje je za pomocą FermionicCircuit.measure_all(), który odczytuje mod i do bitu klasycznego i.

Lepiej mierzyć obwód fermionowy niż transpilowany. Pomiar fermionowy wskazuje nazwę odczytywanego modu, więc układ i routing niosą go wraz z kubitem, na którym ten mod się znalazł, a próbkowany ciąg bitów jest indeksowany według modu. Pomiar obwodu transpilowanego zostawia Ci natomiast samodzielne cofnięcie final_index_layout() i po cichu zwraca spermutowane wyniki, jeśli o tym zapomnisz. Opiera się to na kodowaniu Jordana-Wignera w bazie obsadzeń, którego używa preset, gdzie każdy mod jest niesiony przez jeden kubit.

5. Transpiluj do obwodów kubitowych

Transpilacja FermionicCircuit odwzorowuje jej bramki fermionowe na bramki kubitowe. Wykorzystuje to odwzorowanie Jordana-Wignera jordan_wigner(). Bramki OrbitalRotation są syntezowane do układu cegiełkowego rotacji Givensa, a jednoczłonowa Evolution redukuje się do pojedynczej dwukubitowej bramki oddziaływania.

Ponieważ każda bramka fermionowa ma tu domyślną syntezę, możesz użyć gotowego generate_preset_jw_pass_manager() bezpośrednio, zamiast ręcznie składać etapy. Wszelkie argumenty nazwane są przekazywane do etapu kubitowego; podaj optimization_level=0, aby uzyskać wierny, niezoptymalizowany obraz zsyntezowanej głębokości. Wygenerowanie pełnej rodziny obwodów Kryłowa to wtedy pętla po wymiarze Kryłowa. Zachowaj nietranspilowane FermionicCircuit obok transpilowanych obwodów kubitowych. Te fermionowe napędzają dokładną symulację (wektora stanu) w kroku 6, więc muszą pozostać bez pomiarów; measure_all(inplace=False) zwraca zmierzoną kopię do transpilacji i zostawia oryginał nietknięty, a to właśnie ta kopia byłaby wykonana przez 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)

Wymiar Kryłowa widać bezpośrednio w dwukubitowej głębokości obwodów. Przygotowanie stanu referencyjnego jest kosztem stałym, a każda kolejna potęga Kryłowa dodaje kolejny drugorzędowy krok Trottera operatora ewolucji czasowej.

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

Uwaga

Podział hamiltonianu utrzymuje te głębokości na umiarkowanym poziomie. Ewolucja jednociałowa jest prawdziwą rotacją orbitali, więc syntezuje się do układu cegiełkowego XXPlusYYGate na sąsiednich kubitach. Każdy półkrok oddziaływania lokalnego to pojedyncza RZZGate sprzęgająca dwa mody spinowe domieszki; tani człon, który drugorzędowy iloczyn umieszcza wokół rotacji. Przekazanie całego hamiltonianu do jednej bramki Evolution spowodowałoby natomiast trotteryzację Pauliego części jednociałowej składnik po składniku, propagując długie łańcuchy Z Jordana-Wignera przez (dalekozasięgowe) sprzężenia domieszka-kąpiel i zwiększając głębokość o rząd wielkości.

Pierwszy obwód przygotowuje tylko wyznacznik referencyjny, a ostatni niesie najwięcej kroków Trottera (w tym przypadku cztery):

>>> 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. Próbkuj ciągi bitów z obwodów Kryłowa

Na sprzęcie uruchomiłbyś transpilowane obwody i dokonał pomiaru. To bezszumowa wersja (wektora stanu) tego kroku, pokazująca ciągi bitów, które zebrałby SKQD. Zasymuluj każdy nietranspilowany FermionicCircuit za pomocą ffsim.apply_unitary() i spróbkuj otrzymany wektor stanu za pomocą ffsim.sample_state_vector().

Bramka PrepareSlaterDeterminant w symulacji działa na zasadzie „sprawdź, potem obróć” (sprawdza, czy stan wejściowy zajmuje swój wyznacznik referencyjny, a następnie stosuje rotację). Dlatego wejściem jest zwykły wyznacznik obsadzeń: najniższe mody nocc wypełnione dla każdego spinu, bez rotacji. Bramka sama stosuje 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)

Po zebraniu próbek ze wszystkich pięciu obwodów Kryłowa nośnik jest malutki: kilkaset różnych ciągów bitów spośród $\binom{8}{4}^2 = 4900$ wyznaczników w sektorze (4, 4). Ta koncentracja, bezpośrednia konsekwencja rzadkości w bazie pędowej ustalonej w kroku 2, sprawia, że późniejsza klasyczna diagonalizacja jest wykonalna.

Wykres zliczeń za pomocą qiskit.visualization.plot_histogram() uwidacznia rzadkość rozkładu. Dominuje kilka konfiguracji, na czele z wyznacznikiem referencyjnym.

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

7. Diagonalizacja w próbkowanej podprzestrzeni

Spróbkowane ciągi bitów są wejściem do klasycznej części SKQD. Hamiltonian jest rzutowany na podprzestrzeń rozpinaną przez spróbkowane konfiguracje i tam diagonalizowany. Odbywa się to w pakiecie qiskit-addon-sqd, którego funkcja diagonalize_fermionic_hamiltonian() uruchamia pełną pętlę diagonalizacji kwantowej opartej na próbkowaniu. Buduje ona podprzestrzeń z partii spróbkowanych konfiguracji, diagonalizuje w niej hamiltonian i iteracyjnie udoskonala podprzestrzeń za pomocą odzyskiwania konfiguracji: zmienia obsadzenia, aby naprawić konfiguracje naruszające znaną symetrię liczby cząstek. Odzyskiwanie konfiguracji ma cofać skutki szumu sprzętowego, który psuje spróbkowane ciągi bitów, przenosząc je do niewłaściwego sektora liczby cząstek; tutaj próbkowanie jest bezszumowe, więc nie ma w tym zakresie nic do naprawiania.

Solver przyjmuje spróbkowane zliczenia jako BitArray, zbudowany bezpośrednio ze zbiorczych counts z kroku 6. Przekaż num_bits jawnie: bez tego szerokość jest wnioskowana z najszerszego spróbkowanego ciągu bitów, która jest mniejsza niż 2 * norb, gdy żadna próbka nie obsadza ostatniego modu, a solver odrzuca wtedy tablicę. Solver diagonalizuje w tej samej bazie pędowej, której używa reszta konstrukcji, więc przyjmuje całki h1e_momentum i h2e_momentum z kroku 2 bezpośrednio, bez ponownego pakowania. Odzyskana energia zgadza się z dokładną wartością referencyjną z dokładnością znacznie lepszą niż miliHartree, z zaledwie kilkuset spróbkowanych konfiguracji. Ponieważ energia SQD jest wariacyjnym ograniczeniem górnym energii prawdziwego stanu podstawowego, zawsze leży na poziomie reference_energy (wciąż dostępnej z kroku 2) lub powyżej niego:

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

Zwizualizuj przebieg: zbieżność energii w kolejnych iteracjach oraz średnie obsadzenie każdego orbitalu przestrzennego w odzyskanym stanie podstawowym. Energia w danej iteracji to najniższa wartość spośród partii tej iteracji, a obsadzenie sumuje orbital_occupancies po dwóch sektorach spinowych. Nawet bez szumu do korygowania energia maleje z iteracji na iterację. Każda iteracja diagonalizuje w num_batches partiach zawierających tylko samples_per_batch konfiguracji podpróbkowanych z puli, a iteracyjne udoskonalanie dzięki przenoszeniu ciągów bitów stale poprawia to, które konfiguracje trafiają do tych partii, więc odzyskana energia się poprawia:

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

To domyka pętlę SKQD. Obwody w bazie pędowej z kroków 1–5 dają rzadką próbkę (krok 6), a diagonalizacja oparta na próbkowaniu powyżej zamienia tę próbkę z powrotem w energię stanu podstawowego.

Kolejne kroki

Ten przewodnik uruchomił cały potok SKQD na bezszumowym symulatorze wektora stanu. Na sprzęcie transpilowane circuits byłyby wykonywane i mierzone zamiast próbkowania wektora stanu w kroku 6, a uzyskane zliczenia zasilałyby diagonalizację z kroku 7 bez zmian. Zajrzyj do dokumentacji Qiskit, aby dowiedzieć się, jak uruchamiać obwody, oraz do samouczków dodatku SQD, aby poznać szczegóły dotyczące przetwarzania końcowego z diagonalizacją podprzestrzeni, w tym jego zachowania na zaszumionych próbkach, gdzie odzyskiwanie konfiguracji wykonuje najwięcej pracy.

Przeczytaj przewodnik o relacji z ffsim, aby zrozumieć, dlaczego ffsim.apply_unitary() i ffsim.sample_state_vector() działają natywnie na fermionowych obwodach tego pakietu oraz dlaczego sektor o ustalonej liczbie cząstek użyty w kroku 6 sprawia, że rzadkość w bazie pędowej ustalona w kroku 2 daje efektywne próbkowanie.