Генерація схем часової еволюції Крилова (SKQD)

Важливо

Концепції в цьому посібнику наразі доступні лише в Python API. Еквівалентну функціональність буде додано до C API в одному з майбутніх випусків.

Вибіркова квантова діагоналізація Крилова (SKQD) — це квантово-центричний варіант діагоналізації підпростору. Замість того щоб семплувати один варіаційний анзац (як у звичайній вибірковій квантовій діагоналізації (SQD)) або ансамбль рандомізованих схем часової еволюції (як у SqDRIFT), вона семплує крилівський базис: сімейство станів

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

отриманих еволюцією опорного стану $\lvert \psi_0 \rangle$ протягом дедалі більших кратних фіксованого часового кроку $\Delta t$. Бітові рядки, отримані з цих $D$ схем, заповнюють підпростір Крилова, у якому гамільтоніан згодом діагоналізується класично.

Цей посібник відтворює крок побудови схем алгоритму SKQD для моделі Андерсона з однією домішкою (SIAM) — магнітної домішки, зв'язаної з неінтерагувальною ванною — слідуючи публікації SKQD.

Примітка

Крок семплування SKQD працює, лише якщо основний стан є розрідженим в обчислювальному базисі: класична діагоналізація має бути здатна відновити його з керованої кількості відібраних бітових рядків. Як показано нижче, така розрідженість не з'являється автоматично; це властивість одночастинкового базису, у якому виражено задачу. Правильний вибір базису — суть цієї побудови.

1. Побудова гамільтоніана SIAM як ферміонного оператора

Модель Андерсона з однією домішкою розміщує одну взаємодіючу орбіталь домішки (індекс $0$) у контакті з неінтерагувальною ванною з $L$ вузлів, розташованих у ланцюжок. Її гамільтоніан розкладається на однотільну частину $H_1$ і двотільну частину $H_2$,

$$H = H_1 + H_2,$$

де однотільна частина збирає енергію домішки на вузлі (хімічний потенціал $\mu$), гібридизацію домішки з ванною $V$ та стрибки у ванні між найближчими сусідами $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),$$

а двотільна частина містить внутрівузлове кулонівське відштовхування $U$ лише на домішці,

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

Цей посібник працює в частинково-діркосиметричній точці $\mu = -U/2$. Ці дві частини напряму відповідають однотільним інтегралам h1e (стрибки у ванні, гібридизація домішки, енергія домішки на вузлі) та єдиному двотільному інтегралу h2e (кулонівський доданок домішки). Спершу зберіть їх у позиційному базисі, де модель визначається природно.

>>> 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. Перехід до імпульсного базису

Частина SIAM, що описує ванну, — це ланцюжок вільних ферміонів, тож вона діагоналізується переходом до імпульсного базису. Це вирішальний крок для SKQD. Оскільки ванна трансляційно інваріантна, основний стан не є розрідженим у позиційному базисі, але є розрідженим в імпульсному базисі, де гамільтоніан ванни діагональний. Одночастинкове обертання орбіталей, що виконує цей перехід, — це лінійна алгебра, специфічна для задачі. Воно діагоналізує матрицю стрибків у ванні, не змінюючи орбіталь домішки, а потім переміщує домішку в центральний вузол:

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

На відміну від варіаційного робочого процесу, тут інтеграли гамільтоніана обертаються до імпульсного базису, тож уся схема (і гамільтоніан, і стан) живе в базисі, де основний стан розріджений, а відібрані бітові рядки тому інформативні. Обертання одно- та двотільних інтегралів — стандартна зміна базису тензорів електронної структури:

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

Оскільки обертання не змінює орбіталь домішки (воно діє лише на ванну), внутрівузлова взаємодія не змінюється. h2e_momentum залишається єдиним доданком число-число, так само як h2e. Єдина зміна — домішку переміщено до центрального індексу. Це видно безпосередньо в тензорах інтегралів, кожен з яких має один ненульовий елемент:

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

Ця локальність робить взаємодію дешевою для подальшого синтезу.

Зберіть однотільні інтеграли у FermionOperator за допомогою конструктора електронних інтегралів from_1body_tril_spin_sym(), який очікує їх упакованими в порядку нижнього трикутника. Двотільна частина — це єдиний доданок число-число, тому замість того щоб пропускати її через механізм двотільних інтегралів, побудуйте її напряму за допомогою cre()/ann(): внутрівузловий кулонівський оператор $U \, n_{p\uparrow} n_{p\downarrow}$ для моди домішки $p$ (її мода $\uparrow$ має індекс p, а мода $\downarrow$ — індекс p + norb у блочно-спіновому розташуванні). Обидва конструктори подвоюють за спіном norb просторових орбіталей у 2 * norb ферміонних мод.

Зберігайте однотільну частину $H_1$ і двотільну частину $H_2$ як окремі оператори, а не підсумовуйте їх завчасно, відповідно до доданків $H_1$ і $H_2$ з рівнянь розділу 1. Одноразова побудова кожної частини (в обох базисах) дає змогу використовувати їх повторно нижче. У сумі вони дають повний гамільтоніан для точної діагоналізації в будь-якому базисі. Окремо $H_2$ (імпульсний базис) — це оператор, який крок Троттера еволюціонує на кроці 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)
... )

Зміна базису — це унітарне перетворення подібності, тож вона не змінює спектр. Переконайся, що оператор правильний, порівнявши його найнижче власне значення в секторі $(\text{norb}, \text{nelec})$ з точною діагоналізацією. FermionOperator надає SciPy LinearOperator (на основі нативного ядра добутку матриці на вектор повної конфігураційної взаємодії (FCI)), тож його можна передати прямо в 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

Вигода від зміни базису — розрідженість. Щоб показати це наочно, порахуй, скільки детермінантів обчислювального базису потрібно, щоб охопити 99% ваги основного стану, і порівняй імпульсний базис із позиційним. Енергія не залежить від базису, а от кількість детермінантів залежить:

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

Ці два оператори мають однакову енергію, але основний стан у позиційному базисі розподілений по тисячах детермінантів, тоді як імпульсний базис концентрує його приблизно на 20, і саме тому семплування SKQD працює в імпульсному базисі, а не в позиційному.

3. Підготовка опорного стану

Опорний стан SKQD $\lvert \psi_0 \rangle$ — це суперпозиція всіх збуджень електронів, найближчих до рівня Фермі, у найближчі порожні імпульсні моди. Публікація SKQD будує його вручну написаною мережею XXPlusYYGate. Він виглядає сильно корельованим (має сотні ненульових амплітуд в обчислювальному базисі), але мережа обертань двох мод зі збереженням числа частинок — це одночастинкова (ферміонна гауссова) операція, тож стан, який вона створює, є єдиним детермінантом Слейтера. Це означає, що його можна підготувати одним гейтом PrepareSlaterDeterminant, якщо задано одночастинкове обертання, яке реалізує ця мережа.

Оскільки кожен XXPlusYYGate діє як обертання Гівенса на парі мод, таке обертання можна побудувати безпосередньо в одночастинковій картині (одна матриця Гівенса на гейт), замість того щоб вручну видавати двокубітні гейти. Гейти утворюють цегляну кладку обертань між найближчими сусідами, що розходиться від рівня Фермі розширюваними шарами, розподіляючи кожен приферміївський електрон по порожніх модах безпосередньо над ним:

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

Примітка

Обертання вище виражене в імпульсному базисі, відповідно до гамільтоніана. Електрони поблизу рівня Фермі — це зайняті імпульсні моди з найвищою енергією, і мережа збуджує їх у порожні моди з найнижчою енергією. Саме тому опорний стан потрібно будувати після зміни базису з попереднього кроку, а не до неї.

Підготовка детермінанта Слейтера виконується для кожного спінового сектора окремо, оскільки кожен сектор має власну форму (електрони, орбіталі). Один гейт PrepareSlaterDeterminant (заповнення, найнижчі nocc мод заповнені, разом з одночастинковим обертанням) готує один сектор; той самий гейт застосовується один раз до альфа-мод 0..norb і один раз до бета-мод norb..2*norb, оскільки обидва спінові сектори тут мають те саме опорне заповнення та обертання:

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

Цей декларативний гейт замінює вручну написану мережу XXPlusYYGate; транспайлер генерує синтез обертань Гівенса.

4. Збирання схем Крилова

Кожна схема Крилова — це підготовка опорного стану з подальшою еволюцією під дією $e^{-i k \, \Delta t \, H}$. Розбий гамільтоніан на однотільну та двотільну частини, $H = H_1 + H_2$, і виконай трotterизацію еволюції за цим розбиттям. Це вирішальний вибір для глибини схеми. Однотільна еволюція $e^{-i t H_1}$ — це операція вільних ферміонів (ферміонна гауссова), тож вона є одночастинковим обертанням базису (OrbitalRotation з унітарною матрицею $e^{-i t \, \mathbf{h}_1}$), а не тим, що потрібно трotterизувати по Паулі терм за термом. Двотільна еволюція $e^{-i t H_2}$ діє лише на єдиний внутрівузловий доданок взаємодії. Така побудова кроку дозволяє транспайлеру видати неглибоку цегляну кладку обертань Гівенса плюс один двокубітний гейт взаємодії — неглибоку структуру, яку публікація SKQD будує вручну.

Двотільну частину ми вже маємо. Це оператор two_body_momentum, побудований на кроці 2 (єдиний внутрівузловий доданок). Однотільна еволюція $e^{-i t H_1}$ не залежить від спіну, тож це одне OrbitalRotation на сектор з унітарною матрицею $e^{-i t \, \mathbf{h}_1}$, застосоване один раз до кожного спінового сектора, — та сама паралельна структура, що й у підготовці опорного стану.

Кожна схема Крилова еволюціонує протягом повного часу $k \, \Delta t$. Це реалізується добутком Троттера другого порядку з $k$ кроків розміру $\Delta t$, тож похибка на крок залишається фіксованою зі зростанням розмірності Крилова. Кожен крок затискає повне однотільне обертання (по одному на спіновий сектор) між двома половинами кроку (дешевої, однододанкової) двотільної еволюції. Розмірність Крилова $D$ — це кількість таких схем, де dim пробігає від 0 (лише опорний стан) до $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

Примітка

Збудовані тут схеми поки не містять вимірювань, оскільки ті самі ферміонні схеми також керують точною симуляцією на кроці 6, а вимірювання не є унітарною операцією. Крок 5 додає їх за допомогою FermionicCircuit.measure_all(), що зчитує моду i у класичний біт i.

Краще вимірювати ферміонну схему, а не транспільовану. Ферміонне вимірювання називає моду, яку зчитує, тож розміщення та маршрутизація переносять його разом із кубітом, на якому опинилася ця мода, а відібраний бітовий рядок індексується модою. Натомість вимірювання транспільованої схеми залишає тобі самому скасовувати final_index_layout(), і якщо забути, мовчки повертаються переставлені результати. Це спирається на кодування Йордана-Вігнера в базисі заповнення, яке використовує пресет, де кожну моду несе один кубіт.

5. Транспіляція в кубітні схеми

Транспіляція FermionicCircuit відображає її ферміонні гейти на кубітні гейти. Для цього використовується відображення Йордана-Вігнера jordan_wigner(). Гейти OrbitalRotation синтезуються в цегляну кладку обертань Гівенса, а однододанковий Evolution зводиться до одного двокубітного гейта взаємодії.

Оскільки кожен ферміонний гейт тут має синтез за замовчуванням, ти можеш напряму використати готовий generate_preset_jw_pass_manager(), а не збирати етапи вручну. Будь-які іменовані аргументи передаються на кубітний етап; передай optimization_level=0, щоб побачити вірну, неоптимізовану картину синтезованої глибини. Генерація всього сімейства схем Крилова тоді зводиться до циклу за розмірністю Крилова. Зберігай нетранспільовані FermionicCircuits поруч із транспільованими кубітними схемами. Ферміонні керують точною (statevector) симуляцією на кроці 6, тому вони мають залишатися без вимірювань; measure_all(inplace=False) повертає копію з вимірюваннями для транспіляції та не змінює оригінал, і саме цю копію виконував би бекенд.

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

Розмірність Крилова безпосередньо проявляється в двокубітній глибині схем. Підготовка опорного стану має фіксовану вартість, а кожен додатковий степінь Крилова додає ще один крок Троттера другого порядку оператора часової еволюції.

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

Примітка

Розбиття гамільтоніана утримує ці глибини помірними. Однотільна еволюція — це справжнє обертання орбіталей, тож вона синтезується в цегляну кладку XXPlusYYGate на сусідніх кубітах. Кожна половина кроку внутрівузлової взаємодії — це один RZZGate, що зв'язує дві спінові моди домішки; дешевий доданок, який добуток другого порядку затискає навколо обертання. Передавши весь гамільтоніан одному гейту Evolution, ми натомість виконали б трotterизацію однотільної частини по Паулі терм за термом, поширюючи довгі Z-рядки Йордана-Вігнера через (далекодійні) зв'язки домішки з ванною та роздуваючи глибину на порядок.

Перша схема готує лише опорний детермінант, тоді як остання містить найбільше кроків Троттера (у цьому випадку чотири):

>>> 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. Семплування бітових рядків зі схем Крилова

На реальному обладнанні ти б виконав транспільовані схеми та виміряв. Це безшумна (statevector) версія цього кроку, що показує бітові рядки, які збирав би SKQD. Симулюй кожну нетранспільовану FermionicCircuit за допомогою ffsim.apply_unitary() і семплуй отриманий вектор стану за допомогою ffsim.sample_state_vector().

Гейт PrepareSlaterDeterminant під час симуляції працює за принципом «перевір, потім оберни» (він перевіряє, що вхідний стан займає свій опорний детермінант, а потім застосовує обертання). Тому вхід — це простий детермінант заповнення: найнижчі nocc мод заповнені для кожного спіну, без обертання. Саме обертання 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)

Об'єднавши вибірки з усіх п'яти схем Крилова, отримуємо крихітний носій: кілька сотень різних бітових рядків із $\binom{8}{4}^2 = 4900$ детермінантів у секторі (4, 4). Ця концентрація, прямий наслідок розрідженості в імпульсному базисі, встановленої на кроці 2, робить подальшу класичну діагоналізацію здійсненною.

Побудова гістограми кількостей за допомогою qiskit.visualization.plot_histogram() робить розрідженість видимою. Кілька конфігурацій домінують, і першою серед них є опорний детермінант.

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

7. Діагоналізація у підпросторі зразків

Отримані бітові рядки є вхідними даними для класичної половини SKQD. Гамільтоніан проєктується на підпростір, натягнутий на вибрані конфігурації, і діагоналізується в ньому. Це передається пакету qiskit-addon-sqd, чия функція diagonalize_fermionic_hamiltonian() виконує повний цикл вибіркової квантової діагоналізації. Вона будує підпростір із пакетів вибраних конфігурацій, діагоналізує в ньому гамільтоніан та ітеративно уточнює підпростір за допомогою відновлення конфігурацій: перевертання заповнень, щоб виправити конфігурації, які порушують відому симетрію за числом частинок. Відновлення конфігурацій призначене для усунення шкоди від апаратного шуму, який псує вибрані бітові рядки, переводячи їх у хибний сектор числа частинок; семплування тут безшумне, тож на цьому фронті нічого виправляти.

Розв'язувач приймає вибрані кількості як BitArray, побудований безпосередньо з об'єднаних counts кроку 6. Явно передай num_bits: без нього ширина виводиться з найширшого вибраного бітового рядка, яка вужча за 2 * norb, щойно жоден зразок не займає останню моду, і тоді розв'язувач відхиляє масив. Розв'язувач діагоналізує в тому самому імпульсному базисі, який використовує решта побудови, тому він бере інтеграли h1e_momentum і h2e_momentum з кроку 2 безпосередньо, без перепакування. Відновлена енергія збігається з точним опорним значенням із точністю значно краще міліхартрі, лише з кількох сотень вибраних конфігурацій. Оскільки енергія SQD є варіаційною верхньою межею істинного основного стану, вона завжди лежить на рівні reference_energy або вище (який досі доступний із кроку 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)

Візуалізуй запуск: збіжність енергії за ітераціями та середнє заповнення кожної просторової орбіталі у відновленому основному стані. Енергія на ітерацію є найнижчою серед пакетів цієї ітерації, а заповнення підсумовує orbital_occupancies за двома спіновими секторами. Навіть без шуму, який потрібно виправляти, енергія зменшується від ітерації до ітерації. Кожна ітерація діагоналізує в межах num_batches пакетів лише по samples_per_batch конфігурацій, відібраних із пулу, а ітеративне уточнення завдяки перенесенню бітових рядків стабільно покращує те, які конфігурації потрапляють у ці пакети, тож відновлена енергія покращується:

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

Це завершує цикл SKQD. Кола в імпульсному базисі з кроків 1–5 дають розріджену вибірку (крок 6), а вибіркова діагоналізація вище перетворює цю вибірку назад в енергію основного стану.

Наступні кроки

Цей посібник виконав увесь конвеєр SKQD на безшумному симуляторі вектора стану. На обладнанні транспільовані circuits виконувалися б і вимірювалися замість семплування вектора стану на кроці 6, а отримані кількості без змін подавалися б у діагоналізацію кроку 7. Зверніться до документації Qiskit щодо запуску схем і до навчальних посібників з надбудови SQD щодо постобробки діагоналізації в підпросторі, зокрема її поведінки на зашумлених зразках, де відновлення конфігурацій працює найбільше.

Прочитай посібник про зв'язок з ffsim, щоб зрозуміти, чому ffsim.apply_unitary() і ffsim.sample_state_vector() працюють нативно з ферміонними схемами цього пакета та чому сектор із фіксованим числом частинок, використаний на кроці 6, робить розрідженість в імпульсному базисі, встановлену на кроці 2, ефективною для семплування.