Bumuo ng mga Krylov time-evolution circuit (SKQD)

Important

Ang mga konsepto sa gabay na ito ay kasalukuyang available lang sa Python API. Magiging available ang katumbas na functionality sa C API sa isang susunod na release.

Ang Sample-based Krylov Quantum Diagonalization (SKQD) ay isang quantum-centric na variant ng subspace diagonalization. Sa halip na mag-sample ng iisang variational ansatz (tulad ng sa payak na sample-based quantum diagonalization (SQD)) o ng ensemble ng mga randomized na time-evolution circuit (tulad ng sa SqDRIFT), nagsa-sample ito ng Krylov basis: ang pamilya ng mga state

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

na nakukuha sa pag-evolve ng isang reference state $\lvert \psi_0 \rangle$ para sa pataas na multiple ng isang nakapirming time step $\Delta t$. Ang mga bitstring na na-sample mula sa $D$ circuit na ito ay bumubuo sa Krylov subspace kung saan idi-diagonalize nang klasikal ang Hamiltonian pagkatapos.

Ginagaya ng gabay na ito ang hakbang ng pagbuo ng circuit ng algorithm na SKQD para sa single-impurity Anderson model (SIAM), isang magnetic impurity na nakakabit sa isang non-interacting na bath, kasunod ng publikasyon ng SKQD.

Note

Gumagana lang ang hakbang ng sampling ng SKQD kung sparse ang ground state sa computational basis: dapat kayang buuin ng klasikal na diagonalization ito mula sa kakayanin na bilang ng mga na-sample na bitstring. Tulad ng makikita sa ibaba, hindi awtomatiko ang sparsity na iyon; katangian ito ng single-particle basis kung saan ipinahayag ang problema. Ang pagpili ng tamang basis ang pinakapuso ng konstruksiyong ito.

1. Buuin ang SIAM Hamiltonian bilang fermionic operator

Naglalagay ang single-impurity Anderson model ng isang interacting na impurity orbital (index $0$) na nakadikit sa isang non-interacting na bath na may $L$ na site na nakaayos bilang chain. Nahahati ang Hamiltonian nito sa isang one-body na bahagi na $H_1$ at isang two-body na bahagi na $H_2$,

$$H = H_1 + H_2,$$

kung saan kinokolekta ng one-body na bahagi ang on-site energy ng impurity (ang chemical potential na $\mu$), ang impurity–bath hybridization na $V$, at ang nearest-neighbor na bath hopping na $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),$$

at dala ng two-body na bahagi ang on-site na Coulomb repulsion na $U$ sa impurity lang,

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

Gumagana ang gabay na ito sa particle-hole-symmetric na punto na $\mu = -U/2$. Direktang nagma-map ang dalawang bahaging ito sa mga one-body integral na h1e (bath hopping, impurity hybridization, impurity on-site energy) at sa nag-iisang two-body integral na h2e (ang Coulomb term ng impurity). Buuin muna ang mga ito sa position basis, kung saan natural na tinutukoy ang model.

>>> 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. Lumipat sa momentum basis

Ang bahagi ng bath ng SIAM ay isang free-fermion chain, kaya dina-diagonalize ito ng paglipat sa momentum basis. Ito ang mapagpasyang hakbang para sa SKQD. Dahil translationally invariant ang bath, hindi sparse ang ground state sa position basis, pero sparse ito sa momentum basis, kung saan diagonal ang Hamiltonian ng bath. Ang single-particle orbital rotation na gumagawa ng paglipat na ito ay linear algebra na partikular sa problema. Dina-diagonalize nito ang bath hopping matrix habang hindi ginagalaw ang impurity orbital, saka inililipat ang impurity sa isang sentral na site:

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

Hindi tulad ng variational workflow, ina-rotate nito ang mga Hamiltonian integral papunta sa momentum basis, para ang buong circuit (Hamiltonian at state pareho) ay nasa basis kung saan sparse ang ground state at kaya may saysay ang mga na-sample na bitstring. Ang pag-rotate ng mga one- at two-body integral ay isang karaniwang basis change ng mga tensor ng electronic-structure:

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

Dahil hindi ginagalaw ng rotation ang impurity orbital (sa bath lang ito kumikilos), hindi nagbabago ang on-site na interaksiyon. Nananatiling iisang number-number term ang h2e_momentum, tulad ng h2e. Ang tanging pagbabago ay nailipat ang impurity sa isang sentral na index. Kitang-kita ito mismo sa mga integral tensor, na bawat isa ay may isang nonzero na entry:

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

Ang lokalidad na ito ang nagpapanatiling mura ang pag-synthesize ng interaksiyon sa bandang huli.

I-assemble ang mga one-body integral sa isang FermionOperator gamit ang constructor ng electronic-integral na from_1body_tril_spin_sym(), na umaasang naka-pack ang mga ito sa lower-triangular na ordering. Iisang number-number term ang two-body na bahagi, kaya sa halip na idaan ito sa makinarya ng two-body integral, buuin ito nang direkta gamit ang cre()/ann(): ang on-site na Coulomb operator na $U \, n_{p\uparrow} n_{p\downarrow}$ para sa impurity mode na $p$ (ang $\uparrow$ mode nito sa index na p, ang $\downarrow$ mode nito sa index na p + norb sa block-spin na layout). Dinodoble ng parehong constructor ang spin ng norb na spatial orbital papunta sa 2 * norb na fermionic mode.

Panatilihing magkahiwalay na operator ang one-body na bahagi na $H_1$ at ang two-body na bahagi na $H_2$ sa halip na pagsamahin agad ang mga ito, ayon sa mga term na $H_1$ at $H_2$ ng mga equation sa seksyon 1. Dahil minsan mo lang bubuuin ang bawat piraso (sa parehong basis), magagamit mo itong muli sa ibaba. Kapag pinagsama, ibibigay nila ang buong Hamiltonian para sa exact diagonalization sa alinmang basis. Nang magkahiwalay, ang $H_2$ (momentum basis) ang operator na ine-evolve ng Trotter step sa hakbang 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)
... )

Ang basis change ay isang unitary similarity transform, kaya hindi nito ginagalaw ang spectrum. Tiyaking tama ang operator sa pamamagitan ng paghahambing ng pinakamababang eigenvalue nito sa sector na $(\text{norb}, \text{nelec})$ laban sa exact diagonalization. Nagbibigay ang FermionOperator ng SciPy LinearOperator (na suportado ng native na full configuration interaction (FCI) matrix-vector kernel), kaya puwede itong ipasa nang diretso sa 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

Ang pakinabang ng basis change ay ang sparsity. Para gawin itong konkreto, bilangin kung ilang computational-basis determinant ang kailangan para makuha ang 99% ng bigat ng ground state, at ihambing ang momentum basis sa position basis. Hindi nakadepende sa basis ang energy, pero ang bilang ng mga determinant ay oo:

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

Magkapareho ang energy ng dalawang operator, pero nakakalat ang ground state sa position basis sa libu-libong determinant, habang pinagsisiksik ito ng momentum basis sa humigit-kumulang 20, kaya gumagana ang SKQD sampling sa momentum basis at hindi sa position basis.

3. Ihanda ang reference state

Ang SKQD reference na $\lvert \psi_0 \rangle$ ay ang superposition ng lahat ng excitation ng mga electron na pinakamalapit sa Fermi level papunta sa pinakamalapit na walang lamang momentum mode. Binubuo ito ng publikasyon ng SKQD gamit ang isang kamay na isinulat na network ng mga XXPlusYYGate. Mukha itong malakas ang correlation (mayroon itong daan-daang nonzero na computational-basis amplitude), pero ang network ng mga number-conserving na two-mode rotation ay isang single-particle (fermionic Gaussian) na operasyon, kaya ang state na nagagawa nito ay iisang Slater determinant. Ibig sabihin, maihahanda mo ito gamit ang isang PrepareSlaterDeterminant gate, kung ibibigay ang single-particle rotation na ipinatutupad ng network.

Dahil kumikilos ang bawat XXPlusYYGate bilang Givens rotation sa isang pares ng mode, maaaring buuin ang rotation na iyon nang direkta sa single-particle na pananaw (isang Givens matrix bawat gate) sa halip na isulat nang mano-mano ang mga two-qubit gate. Bumubuo ang mga gate ng brickwork ng mga nearest-neighbor na rotation na kumakalat mula sa Fermi level sa lumalawak na mga layer, na ikinakalat ang bawat electron na malapit sa Fermi sa mga walang lamang mode na nasa itaas lang nito:

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

Note

Ang rotation sa itaas ay ipinahayag sa momentum basis, na tugma sa Hamiltonian. Ang mga electron malapit sa Fermi level ay ang mga occupied na momentum mode na may pinakamataas na energy, at ine-excite sila ng network papunta sa mga walang lamang mode na may pinakamababang energy. Kaya dapat buuin ang reference pagkatapos ng basis change ng nakaraang hakbang, hindi bago nito.

Ginagawa ang paghahanda ng Slater determinant kada spin sector, dahil may sarili nitong hugis (mga electron, mga orbital) ang bawat sector. Isang PrepareSlaterDeterminant gate (isang occupation, na pinunan ang pinakamababang nocc na mode, kasama ang single-particle rotation) ang naghahanda ng isang sector; ang parehong gate ay inilalapat nang isang beses sa mga alpha mode na 0..norb at isang beses sa mga beta mode na norb..2*norb, dahil pareho ang reference occupation at rotation ng dalawang spin sector dito:

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

Pinapalitan ng deklaratibong gate na ito ang kamay na isinulat na network ng XXPlusYYGate; ang transpiler ang bumubuo ng Givens-rotation synthesis.

4. I-assemble ang mga Krylov circuit

Ang bawat Krylov circuit ay ang paghahanda ng reference state na sinusundan ng evolution sa ilalim ng $e^{-i k \, \Delta t \, H}$. Hatiin ang Hamiltonian sa one-body at two-body na bahagi nito, $H = H_1 + H_2$, at i-Trotterize ang evolution ayon sa paghahating iyon. Ito ang mapagpasyang pagpili para sa circuit depth. Ang one-body evolution na $e^{-i t H_1}$ ay isang free-fermion (fermionic Gaussian) na operasyon, kaya isa itong single-particle basis rotation (isang OrbitalRotation na may unitary na $e^{-i t \, \mathbf{h}_1}$) at hindi isang bagay na kailangang i-Pauli-Trotterize bawat term. Ang two-body evolution na $e^{-i t H_2}$ ay kumikilos lang sa iisang on-site na interaction term. Dahil sa ganitong pagbuo ng step, makakapag-emit ang transpiler ng mababaw na brickwork ng mga Givens rotation kasama ang iisang two-qubit na interaction gate, ang mababaw na estrukturang binubuo nang mano-mano ng publikasyon ng SKQD.

Mayroon na tayo ng two-body na bahagi. Ito ang operator na two_body_momentum na binuo sa hakbang 2 (ang iisang on-site na term). Hindi nakadepende sa spin ang one-body evolution na $e^{-i t H_1}$, kaya isa itong OrbitalRotation kada sector na may unitary na $e^{-i t \, \mathbf{h}_1}$ na inilalapat nang isang beses sa bawat spin sector, na kapareho ng parallel na estruktura ng paghahanda ng reference state.

Ang bawat Krylov circuit ay nag-e-evolve para sa kabuuang oras na $k \, \Delta t$. Natutupad ito sa pamamagitan ng second-order na Trotter product ng $k$ na hakbang na may laking $\Delta t$, kaya nananatiling nakapirmi ang error kada hakbang habang lumalaki ang Krylov dimension. Naglalagay ang bawat hakbang ng buong one-body rotation (isa kada spin sector) sa pagitan ng dalawang kalahating hakbang ng (mura, iisang term na) two-body evolution. Ang Krylov dimension na $D$ ay ang bilang ng mga ganoong circuit, na ang dim ay mula 0 (reference state lang) hanggang $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

Note

Wala pang mga measurement ang mga circuit na binuo dito, dahil ang parehong mga fermionic circuit ang nagpapatakbo rin ng exact simulation sa hakbang 6, at hindi unitary na operasyon ang measurement. Idinaragdag ito ng hakbang 5 gamit ang FermionicCircuit.measure_all(), na bumabasa ng mode i papunta sa classical bit i.

Mas mainam na sukatin ang fermionic na circuit kaysa sa na-transpile. Pinangangalanan ng fermionic na measurement ang mode na binabasa nito, kaya dala ito ng layout at routing kasama ang qubit na napuntahan ng mode na iyon, at naka-index ayon sa mode ang na-sample na bitstring. Kapag ang na-transpile na circuit ang sinukat, ikaw na ang kailangang mag-undo ng final_index_layout(), at tahimik itong magbabalik ng mga permuted na resulta kapag nakalimutan mo. Umaasa ito sa occupation-basis na Jordan-Wigner encoding na ginagamit ng preset, kung saan isang qubit ang may dala sa bawat mode.

5. I-transpile sa mga qubit circuit

Ang pag-transpile ng isang FermionicCircuit ay nagma-map sa mga fermionic gate nito papunta sa mga qubit gate. Gumagamit ito ng Jordan-Wigner na mapping na jordan_wigner(). Ang mga OrbitalRotation na gate ay sini-synthesize bilang brickwork ng mga Givens rotation, at ang iisang-term na Evolution ay nagiging iisang two-qubit na interaction gate.

Dahil may default na synthesis ang bawat fermionic gate dito, puwede mong gamitin nang direkta ang handa nang generate_preset_jw_pass_manager() sa halip na i-assemble nang mano-mano ang mga stage. Ipinapasa sa qubit stage ang anumang keyword argument; ipasa ang optimization_level=0 para sa tapat at hindi naka-optimize na larawan ng synthesized na depth. Ang pagbuo ng buong pamilya ng mga Krylov circuit ay isa na lang loop sa Krylov dimension. Itabi ang mga hindi pa na-transpile na FermionicCircuit kasama ng mga na-transpile na qubit circuit. Ang mga fermionic ang nagpapatakbo ng exact (statevector) na simulation sa hakbang 6, kaya dapat manatili silang walang measurement; ang measure_all(inplace=False) ay nagbabalik ng kopyang may measurement na ita-transpile at hindi ginagalaw ang orihinal, at ang kopyang iyon ang isasagawa ng isang 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)

Direktang makikita ang Krylov dimension sa two-qubit depth ng mga circuit. Ang paghahanda ng reference state ay nakapirming gastos, at bawat karagdagang Krylov power ay nagdaragdag ng isa pang second-order na Trotter step ng time-evolution operator.

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

Note

Pinapanatiling katamtaman ng paghahati ng Hamiltonian ang mga depth na ito. Ang one-body evolution ay isang tunay na orbital rotation, kaya sini-synthesize ito bilang XXPlusYYGate na brickwork sa magkakatabing qubit. Ang bawat kalahating hakbang ng on-site na interaksiyon ay iisang RZZGate na nag-uugnay sa dalawang spin mode ng impurity; ang murang term na inilalagay ng second-order na product sa magkabilang gilid ng rotation. Kung ibibigay ang buong Hamiltonian sa iisang Evolution gate, ite-Pauli-Trotterize nito ang one-body na bahagi bawat term, na magpapalaganap ng mahahabang Jordan-Wigner na Z-string sa (long-range) na mga impurity-bath coupling at magpapalobo ng depth nang isang order of magnitude.

Ang unang circuit ay naghahanda lang ng reference determinant, habang ang huli ang may pinakamaraming Trotter step (apat sa kasong ito):

>>> 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. Mag-sample ng mga bitstring mula sa mga Krylov circuit

Sa hardware, patatakbuhin mo ang mga na-transpile na circuit at susukatin. Ito ang noiseless (statevector) na bersyon ng hakbang na iyon, na nagpapakita ng mga bitstring na kokolektahin ng SKQD. I-simulate ang bawat hindi pa na-transpile na FermionicCircuit gamit ang ffsim.apply_unitary() at i-sample ang nagresultang statevector gamit ang ffsim.sample_state_vector().

Ang isang PrepareSlaterDeterminant gate ay validate-then-rotate sa simulation (sinusuri nito na ang papasok na state ay sumasakop sa reference determinant nito, saka ina-apply ang rotation). Kaya, ang input ay ang payak na occupation determinant: ang pinakamababang nocc na mode ay puno kada spin, walang rotation. Ang gate mismo ang nag-a-apply ng 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)

Kapag pinagsama ang mga sample mula sa lahat ng limang Krylov circuit, napakaliit ng support: ilang daang natatanging bitstring mula sa $\binom{8}{4}^2 = 4900$ na determinant sa sector na (4, 4). Ang konsentrasyong ito, direktang bunga ng momentum-basis sparsity na naitatag sa hakbang 2, ang nagpapadali sa susunod na klasikal na diagonalization.

Ang pag-plot ng mga count gamit ang qiskit.visualization.plot_histogram() ay nagpapakita ng sparsity. Ilang configuration lang ang nangingibabaw, pangunahan ng reference determinant.

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

7. I-diagonalize sa sampled subspace

Ang mga na-sample na bitstring ang input sa classical na bahagi ng SKQD. Ipo-project ang Hamiltonian sa subspace na sakop ng mga na-sample na configuration at doon ito dina-diagonalize. Ipinapasa ito sa package na qiskit-addon-sqd, na ang diagonalize_fermionic_hamiltonian() ang nagpapatakbo ng buong sample-based quantum diagonalization loop. Bumubuo ito ng subspace mula sa mga batch ng mga na-sample na configuration, dina-diagonalize ang Hamiltonian doon, at paulit-ulit na pinapaganda ang subspace gamit ang configuration recovery: binabaligtad ang mga occupation para ayusin ang mga configuration na lumalabag sa kilalang particle-number symmetry. Ang configuration recovery ay idinisenyo para ayusin ang pinsala ng hardware noise, na nagpapasira sa mga na-sample na bitstring at inililipat ang mga ito sa maling particle-number sector; noiseless ang sampling dito, kaya walang kailangang ayusin sa aspetong iyon.

Kinukuha ng solver ang mga na-sample na count bilang BitArray, na direktang binuo mula sa pinagsamang counts ng step 6. Ipasa nang tahasan ang num_bits: kung wala ito, hinuhulaan ang lapad mula sa pinakamalapad na na-sample na bitstring, na mas makitid kaysa sa 2 * norb kapag walang sample na sumasakop sa huling mode, at tatanggihan na ng solver ang array. Nagdi-diagonalize ang solver sa parehong momentum basis na ginagamit ng natitirang bahagi ng construction, kaya direkta nitong kinukuha ang mga integral na h1e_momentum at h2e_momentum mula sa step 2, nang walang repacking. Ang nabawing energy ay tumutugma sa eksaktong reference hanggang sa mas mababa pa sa isang milli-Hartree, mula lang sa ilang daang na-sample na configuration. Dahil ang SQD energy ay isang variational upper bound sa tunay na ground state, laging nasa o mas mataas ito sa reference_energy (nasa scope pa rin mula sa step 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)

I-visualize ang takbo: ang convergence ng energy sa bawat iteration, at ang average na occupancy ng bawat spatial orbital sa nabawing ground state. Ang energy kada iteration ay ang pinakamababa sa mga batch ng iteration na iyon, at pinagsasama ng occupancy ang orbital_occupancies sa dalawang spin sector. Kahit walang noise na itatama, bumababa ang energy sa bawat iteration. Nagdi-diagonalize ang bawat iteration sa loob ng num_batches batch na may samples_per_batch configuration lang na sinubsample mula sa pool, at ang paulit-ulit na refinement dahil sa bitstring carryover ay tuloy-tuloy na nagpapabuti kung aling mga configuration ang napupunta sa mga batch na iyon, kaya gumaganda ang nabawing energy:

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

Dito nagtatapos ang SKQD loop. Ang mga momentum-basis circuit ng steps 1–5 ay gumagawa ng sparse na sample (step 6), at ang sample-based diagonalization sa itaas ang nagbabalik sa sample na iyon bilang ground-state energy.

Mga susunod na hakbang

Pinatakbo ng gabay na ito ang buong SKQD pipeline sa isang noiseless statevector simulator. Sa hardware, ie-execute at susukatin ang na-transpile na circuits kapalit ng statevector sampling sa step 6, at ang mga resultang count ay ipapasok sa step-7 diagonalization nang walang pagbabago. Tingnan ang Qiskit documentation para sa tulong sa pagpapatakbo ng mga circuit, at ang SQD addon tutorials para sa impormasyon tungkol sa subspace-diagonalization post-processing, kasama ang pag-uugali nito sa mga noisy na sample, kung saan pinakamalaki ang trabaho ng configuration recovery.

Basahin ang ffsim relationship guide para maintindihan kung bakit gumagana nang native ang ffsim.apply_unitary() at ffsim.sample_state_vector() sa mga fermionic circuit ng package na ito, at kung bakit ang fixed-particle-number sector na ginamit sa step 6 ay ginagawang sample-efficient ang momentum-basis sparsity na itinatag sa step 2.