Buat circuit evolusi waktu Krylov (SKQD)

Penting

Konsep di panduan ini saat ini hanya tersedia di Python API. Fungsionalitas yang setara akan tersedia di C API pada rilis mendatang.

Sample-based Krylov Quantum Diagonalization (SKQD) adalah varian diagonalisasi subruang yang berpusat pada kuantum. Alih-alih men-sample satu ansatz variasional (seperti pada sample-based quantum diagonalization (SQD) biasa) atau ensemble circuit evolusi waktu acak (seperti pada SqDRIFT), metode ini men-sample basis Krylov: keluarga keadaan

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

yang diperoleh dengan mengevolusikan keadaan referensi $\lvert \psi_0 \rangle$ selama kelipatan yang makin besar dari langkah waktu tetap $\Delta t$. Bitstring yang di-sample dari $D$ circuit ini mengisi subruang Krylov tempat Hamiltonian kemudian didiagonalisasi secara klasik.

Panduan ini mereproduksi langkah konstruksi circuit dari algoritma SKQD untuk single-impurity Anderson model (SIAM), yaitu sebuah impuritas magnetik yang terhubung ke bath tanpa interaksi, mengikuti publikasi SKQD.

Catatan

Langkah sampling SKQD hanya bekerja kalau ground state-nya sparse di basis komputasional: diagonalisasi klasik harus bisa merekonstruksinya dari sejumlah bitstring hasil sampling yang masih terkelola. Seperti yang ditunjukkan di bawah, sparsity itu tidak otomatis ada; ia adalah properti dari basis partikel tunggal tempat masalahnya dinyatakan. Memilih basis yang tepat adalah kunci dari konstruksi ini.

1. Bangun Hamiltonian SIAM sebagai operator fermionik

Single-impurity Anderson model menempatkan satu orbital impuritas yang berinteraksi (indeks $0$) bersentuhan dengan bath tanpa interaksi berisi $L$ situs yang tersusun sebagai rantai. Hamiltoniannya terbagi menjadi bagian satu-badan $H_1$ dan bagian dua-badan $H_2$,

$$H = H_1 + H_2,$$

dengan bagian satu-badan mengumpulkan energi on-site impuritas (potensial kimia $\mu$), hibridisasi impuritas–bath $V$, dan hopping bath tetangga terdekat $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),$$

dan bagian dua-badan membawa tolakan Coulomb on-site $U$ pada impuritas saja,

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

Panduan ini bekerja di titik simetris partikel-lubang $\mu = -U/2$. Kedua bagian ini terpetakan langsung ke integral satu-badan h1e (hopping bath, hibridisasi impuritas, energi on-site impuritas) dan satu integral dua-badan h2e (suku Coulomb impuritas). Rakit keduanya dulu di basis posisi, tempat model ini secara alami didefinisikan.

>>> 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. Pindah ke basis momentum

Bagian bath dari SIAM adalah rantai fermion bebas, jadi ia didiagonalisasi dengan pindah ke basis momentum. Ini langkah penentu untuk SKQD. Karena bath-nya invarian terhadap translasi, ground state tidak sparse di basis posisi, tapi sparse di basis momentum, tempat Hamiltonian bath-nya diagonal. Rotasi orbital partikel tunggal yang melakukan perpindahan ini adalah aljabar linear khusus masalah. Rotasi ini mendiagonalisasi matriks hopping bath sambil membiarkan orbital impuritas tidak tersentuh, lalu memindahkan impuritas ke situs tengah:

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

Berbeda dengan workflow variasional, di sini integral Hamiltonian-lah yang dirotasi ke basis momentum, sehingga seluruh circuit (baik Hamiltonian maupun keadaannya) berada di basis tempat ground state sparse dan bitstring hasil sampling karena itu informatif. Merotasi integral satu- dan dua-badan adalah perubahan basis standar untuk tensor struktur elektronik:

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

Karena rotasi membiarkan orbital impuritas tetap (ia hanya bekerja pada bath), interaksi on-site tidak berubah. h2e_momentum tetap berupa satu suku number-number, persis seperti h2e. Satu-satunya perubahan adalah impuritas dipindahkan ke indeks tengah. Ini terlihat langsung di tensor integralnya, yang masing-masing membawa satu entri tak nol:

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

Lokalitas inilah yang membuat interaksi tetap murah untuk disintesis nanti.

Rakit integral satu-badan menjadi sebuah FermionOperator dengan konstruktor integral elektronik from_1body_tril_spin_sym(), yang mengharapkan integral dikemas dalam urutan segitiga bawah. Bagian dua-badan adalah satu suku number-number, jadi daripada melewatkannya lewat mesin integral dua-badan, bangun langsung dengan cre()/ann(): operator Coulomb on-site $U \, n_{p\uparrow} n_{p\downarrow}$ untuk mode impuritas $p$ (mode $\uparrow$-nya di indeks p, mode $\downarrow$-nya di indeks p + norb pada tata letak block-spin). Kedua konstruktor menggandakan spin norb orbital spasial menjadi 2 * norb mode fermionik.

Pertahankan bagian satu-badan $H_1$ dan bagian dua-badan $H_2$ sebagai operator terpisah, jangan langsung dijumlahkan, sesuai dengan suku $H_1$ dan $H_2$ pada persamaan di bagian 1. Membangun tiap bagian sekali saja (di kedua basis) membuatnya bisa dipakai ulang di bawah. Kalau dijumlahkan, keduanya membentuk Hamiltonian penuh untuk diagonalisasi eksak di basis mana pun. Secara terpisah, $H_2$ (basis momentum) adalah operator yang dievolusikan oleh langkah Trotter di langkah 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)
... )

Perubahan basis adalah transformasi kesamaan uniter, jadi spektrumnya tidak berubah. Pastikan operatornya benar dengan membandingkan eigenvalue terendahnya di sektor $(\text{norb}, \text{nelec})$ terhadap diagonalisasi eksak. FermionOperator menyediakan SciPy LinearOperator (didukung kernel perkalian matriks-vektor full configuration interaction (FCI) native), sehingga bisa langsung diberikan ke 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

Keuntungan dari perubahan basis adalah sparsity. Biar konkret, hitung berapa banyak determinan basis komputasional yang dibutuhkan untuk menangkap 99% bobot ground state, lalu bandingkan basis momentum dengan basis posisi. Energinya tidak bergantung pada basis, tapi jumlah determinannya tidak demikian:

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

Kedua operator punya energi yang sama, tapi ground state di basis posisi tersebar di ribuan determinan, sedangkan basis momentum memusatkannya ke sekitar 20, makanya sampling SKQD bekerja di basis momentum dan bukan di basis posisi.

3. Siapkan keadaan referensi

Referensi SKQD $\lvert \psi_0 \rangle$ adalah superposisi dari semua eksitasi elektron yang paling dekat dengan level Fermi ke mode momentum kosong terdekat. Publikasi SKQD membangunnya dengan jaringan XXPlusYYGate yang ditulis manual. Keadaan ini tampak berkorelasi kuat (punya ratusan amplitudo basis komputasional tak nol), tapi jaringan rotasi dua-mode yang menjaga jumlah partikel adalah operasi partikel tunggal (Gaussian fermionik), sehingga keadaan yang dihasilkannya adalah satu determinan Slater. Artinya, kamu bisa menyiapkannya dengan satu gate PrepareSlaterDeterminant, asalkan rotasi partikel tunggal yang diimplementasikan jaringan itu diketahui.

Karena setiap XXPlusYYGate bekerja sebagai rotasi Givens pada sepasang mode, rotasi itu bisa dibangun langsung di gambaran partikel tunggal (satu matriks Givens per gate) daripada mengeluarkan gate dua-qubit secara manual. Gate-gate itu membentuk brickwork rotasi tetangga terdekat yang menyebar dari level Fermi dalam lapisan-lapisan yang meluas, menyebarkan tiap elektron dekat-Fermi ke mode kosong tepat di atasnya:

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

Catatan

Rotasi di atas dinyatakan di basis momentum, sesuai dengan Hamiltonian. Elektron di dekat level Fermi adalah mode momentum terisi berenergi tertinggi, dan jaringan itu mengeksitasi mereka ke mode kosong berenergi terendah. Inilah sebabnya referensi harus dibangun setelah perubahan basis di langkah sebelumnya, bukan sebelumnya.

Persiapan determinan Slater dilakukan per sektor spin, karena tiap sektor punya bentuk (elektron, orbital) sendiri. Satu gate PrepareSlaterDeterminant (sebuah okupansi, yaitu nocc mode terendah terisi, bersama rotasi partikel tunggal) menyiapkan satu sektor; gate yang sama diterapkan sekali pada mode alpha 0..norb dan sekali pada mode beta norb..2*norb, karena kedua sektor spin di sini berbagi okupansi referensi dan rotasi yang sama:

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

Gate deklaratif ini menggantikan jaringan XXPlusYYGate yang di-coding manual; transpiler yang menghasilkan sintesis rotasi Givens.

4. Rakit circuit Krylov

Setiap circuit Krylov adalah persiapan keadaan referensi diikuti evolusi di bawah $e^{-i k \, \Delta t \, H}$. Pisahkan Hamiltonian menjadi bagian satu-badan dan dua-badan, $H = H_1 + H_2$, lalu lakukan Trotterisasi evolusi atas pemisahan itu. Ini pilihan penentu untuk kedalaman circuit. Evolusi satu-badan $e^{-i t H_1}$ adalah operasi fermion bebas (Gaussian fermionik), jadi ia berupa rotasi basis partikel tunggal (sebuah OrbitalRotation dengan uniter $e^{-i t \, \mathbf{h}_1}$) dan bukan sesuatu yang harus di-Trotterisasi Pauli suku demi suku. Evolusi dua-badan $e^{-i t H_2}$ hanya bekerja pada satu suku interaksi on-site. Dengan membangun langkah seperti ini, transpiler bisa menghasilkan brickwork rotasi Givens yang dangkal plus satu gate interaksi dua-qubit, struktur dangkal yang dibangun secara manual oleh publikasi SKQD.

Kita sudah punya bagian dua-badan. Itu adalah operator two_body_momentum yang dibangun di langkah 2 (satu suku on-site). Evolusi satu-badan $e^{-i t H_1}$ tidak bergantung pada spin, jadi ia berupa satu OrbitalRotation per sektor dengan uniter $e^{-i t \, \mathbf{h}_1}$ yang diterapkan sekali ke tiap sektor spin, dengan struktur paralel yang sama seperti persiapan keadaan referensi.

Setiap circuit Krylov berevolusi selama total waktu $k \, \Delta t$. Itu diwujudkan dengan produk Trotter orde kedua dari $k$ langkah berukuran $\Delta t$, sehingga galat per langkah tetap konstan seiring bertambahnya dimensi Krylov. Setiap langkah mengapit satu rotasi satu-badan penuh (satu per sektor spin) di antara dua setengah-langkah evolusi dua-badan (yang murah, satu suku). Dimensi Krylov $D$ adalah jumlah circuit semacam itu, dengan dim berkisar dari 0 (hanya keadaan referensi) sampai $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

Catatan

Circuit yang dibangun di sini belum membawa pengukuran, karena circuit fermionik yang sama juga dipakai untuk simulasi eksak di langkah 6, dan pengukuran bukan operasi uniter. Langkah 5 menambahkannya dengan FermionicCircuit.measure_all(), yang membaca mode i ke bit klasik i.

Lebih baik mengukur circuit fermionik daripada circuit hasil transpile. Pengukuran fermionik menyebut mode yang dibacanya, sehingga layout dan routing membawanya bersama qubit tempat mode itu berakhir, dan bitstring hasil sampling diindeks berdasarkan mode. Kalau kamu mengukur circuit hasil transpile, kamu harus membatalkan final_index_layout() sendiri, dan hasilnya diam-diam akan tertukar kalau kamu lupa. Ini mengandalkan encoding Jordan-Wigner basis okupansi yang dipakai preset, di mana tiap mode dibawa oleh satu qubit.

5. Transpile ke circuit qubit

Men-transpile sebuah FermionicCircuit memetakan gate fermioniknya ke gate qubit. Ini memakai pemetaan Jordan-Wigner jordan_wigner(). Gate OrbitalRotation disintesis menjadi brickwork rotasi Givens, dan Evolution satu suku diturunkan menjadi satu gate interaksi dua-qubit.

Karena setiap gate fermionik di sini punya sintesis default, kamu bisa langsung memakai generate_preset_jw_pass_manager() yang sudah jadi daripada merakit tahapannya sendiri. Keyword argument apa pun diteruskan ke tahap qubit; berikan optimization_level=0 untuk gambaran kedalaman hasil sintesis yang apa adanya dan tanpa optimasi. Menghasilkan seluruh keluarga circuit Krylov tinggal berupa loop atas dimensi Krylov. Simpan FermionicCircuit yang belum di-transpile berdampingan dengan circuit qubit hasil transpile. Yang fermionik menjalankan simulasi eksak (statevector) di langkah 6, jadi harus tetap bebas pengukuran; measure_all(inplace=False) mengembalikan salinan berpengukuran untuk di-transpile dan membiarkan aslinya tidak tersentuh, dan salinan itulah yang akan dieksekusi 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)

Dimensi Krylov terlihat langsung pada kedalaman dua-qubit circuit. Persiapan keadaan referensi adalah biaya tetap, dan setiap pangkat Krylov tambahan menambah satu langkah Trotter orde kedua lagi dari operator evolusi waktu.

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

Catatan

Memisahkan Hamiltonian menjaga kedalaman ini tetap sedang. Evolusi satu-badan adalah rotasi orbital sejati, jadi ia disintesis menjadi brickwork XXPlusYYGate pada qubit yang bersebelahan. Setiap setengah-langkah interaksi on-site adalah satu RZZGate yang menghubungkan dua mode spin impuritas; suku murah yang diapit produk orde kedua di sekitar rotasi. Kalau seluruh Hamiltonian diberikan ke satu gate Evolution, bagian satu-badan akan di-Trotterisasi Pauli suku demi suku, merambatkan Z-string Jordan-Wigner yang panjang melintasi kopling impuritas-bath (jarak jauh) dan menggelembungkan kedalaman hingga satu orde magnitudo.

Circuit pertama hanya menyiapkan determinan referensi, sedangkan yang terakhir membawa langkah Trotter terbanyak (empat dalam kasus ini):

>>> 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. Sample bitstring dari circuit Krylov

Di hardware, kamu akan menjalankan circuit hasil transpile lalu mengukurnya. Ini adalah versi tanpa noise (statevector) dari langkah itu, yang menunjukkan bitstring yang akan dikumpulkan SKQD. Simulasikan tiap FermionicCircuit yang belum di-transpile dengan ffsim.apply_unitary() lalu sample statevector hasilnya dengan ffsim.sample_state_vector().

Gate PrepareSlaterDeterminant bersifat validasi-lalu-rotasi saat simulasi (ia memeriksa bahwa keadaan masukan menempati determinan referensinya, lalu menerapkan rotasi). Karena itu, inputnya adalah determinan okupansi biasa: nocc mode terendah terisi per spin, tanpa rotasi. Gate itu menerapkan reference_rotation sendiri.

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

Dengan menggabungkan sampel dari kelima circuit Krylov, dukungannya sangat kecil: beberapa ratus bitstring berbeda dari $\binom{8}{4}^2 = 4900$ determinan di sektor (4, 4). Konsentrasi ini, konsekuensi langsung dari sparsity basis momentum yang ditetapkan di langkah 2, membuat diagonalisasi klasik berikutnya tertangani.

Memplot hasil hitungan dengan qiskit.visualization.plot_histogram() membuat sparsitasnya terlihat jelas. Beberapa konfigurasi mendominasi, dipimpin oleh determinan referensi.

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

7. Diagonalisasi dalam subruang hasil sampling

Bitstring hasil sampling menjadi input untuk bagian klasik dari SKQD. Hamiltonian diproyeksikan ke subruang yang direntang oleh konfigurasi hasil sampling, lalu didiagonalisasi di sana. Ini diteruskan ke paket qiskit-addon-sqd, yang fungsi diagonalize_fermionic_hamiltonian()-nya menjalankan seluruh loop sample-based quantum diagonalization. Fungsi ini membangun subruang dari batch konfigurasi hasil sampling, mendiagonalisasi Hamiltonian di dalamnya, dan menyempurnakan subruang secara iteratif dengan configuration recovery: membalik okupansi untuk memperbaiki konfigurasi yang melanggar simetri jumlah partikel yang diketahui. Configuration recovery dirancang untuk membalikkan kerusakan akibat noise hardware, yang merusak bitstring hasil sampling sehingga masuk ke sektor jumlah partikel yang salah; sampling di sini bebas noise, jadi tidak ada yang perlu diperbaiki dari sisi itu.

Solver memakai hasil hitungan sampling sebagai BitArray, yang dibangun langsung dari counts gabungan di langkah 6. Berikan num_bits secara eksplisit: tanpa itu, lebarnya disimpulkan dari bitstring hasil sampling terlebar, yang lebih sempit dari 2 * norb kapan pun tidak ada sampel yang menempati mode terakhir, dan solver lalu menolak array tersebut. Solver mendiagonalisasi dalam basis momentum yang sama dengan yang dipakai di seluruh konstruksi lainnya, jadi solver langsung memakai integral h1e_momentum dan h2e_momentum dari langkah 2, tanpa pengemasan ulang. Energi hasil pemulihan cocok dengan referensi eksak hingga jauh di bawah satu milli-Hartree, dari hanya beberapa ratus konfigurasi hasil sampling. Karena energi SQD adalah batas atas variasional untuk ground state sebenarnya, nilainya selalu sama dengan atau di atas reference_energy (yang masih tersedia dari langkah 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)

Visualisasikan hasil run: konvergensi energi di sepanjang iterasi, dan okupansi rata-rata tiap orbital spasial pada ground state hasil pemulihan. Energi per iterasi adalah yang terendah dari semua batch pada iterasi tersebut, dan okupansi menjumlahkan orbital_occupancies atas kedua sektor spin. Bahkan tanpa noise yang perlu dikoreksi, energi menurun dari iterasi ke iterasi. Setiap iterasi mendiagonalisasi dalam num_batches batch yang masing-masing hanya berisi samples_per_batch konfigurasi yang disubsampel dari kumpulan, dan penyempurnaan iteratif berkat bitstring carryover terus memperbaiki konfigurasi mana yang masuk ke batch tersebut, sehingga energi hasil pemulihan membaik:

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

Ini menutup loop SKQD. Circuit basis momentum dari langkah 1–5 menghasilkan sampel yang sparse (langkah 6), dan diagonalisasi berbasis sampel di atas mengubah sampel itu kembali menjadi energi ground state.

Langkah berikutnya

Panduan ini menjalankan seluruh pipeline SKQD pada simulator statevector tanpa noise. Pada hardware, circuits yang sudah ditranspilasi akan dieksekusi dan diukur sebagai pengganti sampling statevector di langkah 6, dengan hasil hitungan yang diperoleh masuk ke diagonalisasi langkah 7 tanpa perubahan. Lihat dokumentasi Qiskit untuk bantuan menjalankan circuit, dan tutorial addon SQD untuk informasi tentang post-processing diagonalisasi subruang, termasuk perilakunya pada sampel yang noisy, di mana configuration recovery bekerja paling banyak.

Baca panduan hubungan ffsim untuk memahami mengapa ffsim.apply_unitary() dan ffsim.sample_state_vector() bekerja secara native pada circuit fermionik paket ini, dan mengapa sektor jumlah partikel tetap yang dipakai di langkah 6 membuat sparsitas basis momentum yang ditetapkan di langkah 2 efisien dalam sampling.