Krylov-Zeitentwicklungs-Circuits (SKQD) erzeugen

Wichtig

Die Konzepte in diesem Leitfaden sind derzeit nur in der Python-API verfügbar. Entsprechende Funktionalität wird in einer zukünftigen Version in der C-API bereitgestellt.

Sample-based Krylov Quantum Diagonalization (SKQD) ist eine quantenzentrierte Variante der Unterraumdiagonalisierung. Statt wie bei der einfachen Sample-based Quantum Diagonalization (SQD) einen einzelnen variationellen Ansatz zu samplen oder wie bei SqDRIFT ein Ensemble randomisierter Zeitentwicklungs-Circuits, wird hier eine Krylov-Basis gesampelt: die Familie von Zuständen

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

die durch Zeitentwicklung eines Referenzzustands $\lvert \psi_0 \rangle$ für wachsende Vielfache eines festen Zeitschritts $\Delta t$ entstehen. Bitstrings, die aus diesen $D$ Circuits gesampelt werden, besetzen den Krylov-Unterraum, in dem der Hamiltonoperator anschließend klassisch diagonalisiert wird.

Dieser Leitfaden reproduziert den Schritt der Circuit-Konstruktion des SKQD-Algorithmus für das Einzelstörstellen-Anderson-Modell (SIAM), eine magnetische Störstelle, die an ein nicht wechselwirkendes Bad gekoppelt ist, gemäß der SKQD-Publikation.

Hinweis

Der SKQD-Sampling-Schritt funktioniert nur, wenn der Grundzustand in der Rechenbasis dünnbesetzt ist: Die klassische Diagonalisierung muss ihn aus einer überschaubaren Anzahl gesampelter Bitstrings rekonstruieren können. Wie unten gezeigt wird, ist diese Dünnbesetztheit nicht selbstverständlich; sie ist eine Eigenschaft der Einteilchenbasis, in der das Problem ausgedrückt wird. Die richtige Basis zu wählen ist der Kern dieser Konstruktion.

1. Den SIAM-Hamiltonoperator als fermionischen Operator aufbauen

Das Einzelstörstellen-Anderson-Modell platziert ein wechselwirkendes Störstellenorbital (Index $0$) in Kontakt mit einem nicht wechselwirkenden Bad aus $L$ Plätzen, die als Kette angeordnet sind. Sein Hamiltonoperator zerfällt in einen Einteilchenanteil $H_1$ und einen Zweiteilchenanteil $H_2$,

$$H = H_1 + H_2,$$

wobei der Einteilchenanteil die On-Site-Energie der Störstelle (das chemische Potenzial $\mu$), die Hybridisierung $V$ zwischen Störstelle und Bad sowie das Hopping $t$ zwischen nächsten Nachbarn im Bad bündelt,

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

und der Zweiteilchenanteil die On-Site-Coulomb-Abstoßung $U$ allein auf der Störstelle trägt,

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

Dieser Leitfaden arbeitet am teilchen-loch-symmetrischen Punkt $\mu = -U/2$. Diese beiden Anteile entsprechen direkt den Einteilchenintegralen h1e (Bad-Hopping, Störstellen-Hybridisierung, On-Site-Energie der Störstelle) und dem einzelnen Zweiteilchenintegral h2e (dem Coulomb-Term der Störstelle). Stelle sie zunächst in der Ortsbasis zusammen, in der das Modell natürlicherweise definiert ist.

>>> 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. Zur Impulsbasis wechseln

Der Badanteil des SIAM ist eine Kette freier Fermionen und wird daher durch einen Wechsel zur Impulsbasis diagonalisiert. Das ist der entscheidende Schritt für SKQD. Da das Bad translationsinvariant ist, ist der Grundzustand in der Ortsbasis nicht dünnbesetzt, wohl aber in der Impulsbasis, in der der Bad-Hamiltonoperator diagonal ist. Die Einteilchen-Orbitalrotation, die diesen Wechsel vollzieht, ist problemspezifische lineare Algebra. Sie diagonalisiert die Hopping-Matrix des Bades, lässt das Störstellenorbital unberührt und verlegt die Störstelle dann an einen zentralen Platz:

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

Anders als bei einem variationellen Workflow rotiert dies die Integrale des Hamiltonoperators in die Impulsbasis, sodass der gesamte Circuit (Hamiltonoperator und Zustand gleichermaßen) in der Basis lebt, in der der Grundzustand dünnbesetzt ist und die gesampelten Bitstrings daher aussagekräftig sind. Die Rotation der Ein- und Zweiteilchenintegrale ist ein Standard-Basiswechsel der Tensoren der Elektronenstruktur:

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

Da die Rotation das Störstellenorbital unverändert lässt (sie wirkt nur auf das Bad), bleibt die On-Site-Wechselwirkung unverändert. h2e_momentum bleibt ein einzelner Zahl-Zahl-Term, genau wie h2e. Die einzige Änderung ist, dass die Störstelle an einen zentralen Index verlegt wurde. Das ist direkt in den Integraltensoren zu sehen, die jeweils einen von null verschiedenen Eintrag tragen:

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

Diese Lokalität sorgt dafür, dass sich die Wechselwirkung später günstig synthetisieren lässt.

Stelle die Einteilchenintegrale mit dem Konstruktor für elektronische Integrale from_1body_tril_spin_sym() zu einem FermionOperator zusammen; er erwartet sie in der Anordnung als untere Dreiecksmatrix. Der Zweiteilchenanteil ist ein einzelner Zahl-Zahl-Term. Statt ihn über die Zweiteilchenintegral-Maschinerie zu leiten, baust du ihn daher direkt mit cre()/ann() auf: der On-Site-Coulomb-Operator $U \, n_{p\uparrow} n_{p\downarrow}$ für die Störstellenmode $p$ (ihre $\uparrow$-Mode bei Index p, ihre $\downarrow$-Mode bei Index p + norb im Blockspin-Layout). Beide Konstruktoren verdoppeln die norb räumlichen Orbitale spinweise zu 2 * norb fermionischen Moden.

Behalte den Einteilchenanteil $H_1$ und den Zweiteilchenanteil $H_2$ als getrennte Operatoren, statt sie vorschnell zu summieren, passend zu den Termen $H_1$ und $H_2$ der Gleichungen aus Abschnitt 1. Wenn du jedes Stück einmal (in beiden Basen) aufbaust, kannst du sie unten wiederverwenden. Summiert ergeben sie den vollständigen Hamiltonoperator für die exakte Diagonalisierung in beiden Basen. Einzeln ist $H_2$ (Impulsbasis) der Operator, unter dem sich der Trotter-Schritt in Schritt 4 entwickelt.

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

Der Basiswechsel ist eine unitäre Ähnlichkeitstransformation und lässt das Spektrum daher unberührt. Bestätige, dass der Operator korrekt ist, indem du seinen niedrigsten Eigenwert im Sektor $(\text{norb}, \text{nelec})$ mit einer exakten Diagonalisierung vergleichst. Der FermionOperator stellt einen SciPy-LinearOperator bereit (gestützt auf einen nativen Matrix-Vektor-Kernel der vollständigen Konfigurationswechselwirkung (FCI)), sodass er direkt an scipy.sparse.linalg.eigsh() übergeben werden kann.

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

Der Gewinn des Basiswechsels ist die Dünnbesetztheit. Um das greifbar zu machen, zähle, wie viele Determinanten der Rechenbasis nötig sind, um 99 % des Gewichts des Grundzustands zu erfassen, und vergleiche die Impulsbasis mit der Ortsbasis. Die Energie ist basisunabhängig, die Anzahl der Determinanten aber nicht:

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

Die beiden Operatoren haben dieselbe Energie, doch der Grundzustand in der Ortsbasis verteilt sich auf Tausende von Determinanten, während ihn die Impulsbasis auf etwa 20 konzentriert. Deshalb funktioniert das SKQD-Sampling in der Impulsbasis und nicht in der Ortsbasis.

3. Den Referenzzustand vorbereiten

Die SKQD-Referenz $\lvert \psi_0 \rangle$ ist die Superposition aller Anregungen der Elektronen, die der Fermi-Energie am nächsten liegen, in die nächstgelegenen leeren Impulsmoden. Die SKQD-Publikation baut sie mit einem von Hand geschriebenen Netzwerk aus XXPlusYYGates auf. Sie wirkt stark korreliert (sie hat Hunderte von Amplituden in der Rechenbasis, die von null verschieden sind), aber ein Netzwerk aus teilchenzahlerhaltenden Zweimodenrotationen ist eine Einteilchenoperation (fermionisch gaußsch), sodass der damit erzeugte Zustand eine einzelne Slater-Determinante ist. Das heißt, du kannst ihn mit einem einzigen PrepareSlaterDeterminant-Gate vorbereiten, wenn du die Einteilchenrotation kennst, die das Netzwerk implementiert.

Da jedes XXPlusYYGate als Givens-Rotation auf einem Modenpaar wirkt, lässt sich diese Rotation direkt im Einteilchenbild aufbauen (eine Givens-Matrix pro Gate), statt die Zwei-Qubit-Gates von Hand auszugeben. Die Gates bilden ein Brickwork aus Rotationen zwischen nächsten Nachbarn, das sich in wachsenden Lagen von der Fermi-Energie aus auffächert und jedes Elektron nahe der Fermi-Energie auf die leeren Moden direkt darüber verteilt:

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

Hinweis

Die obige Rotation ist in der Impulsbasis ausgedrückt, passend zum Hamiltonoperator. Die Elektronen nahe der Fermi-Energie sind die höchstenergetischen besetzten Impulsmoden, und das Netzwerk regt sie in die niedrigstenergetischen leeren an. Deshalb muss die Referenz nach dem Basiswechsel des vorherigen Schritts aufgebaut werden, nicht davor.

Die Vorbereitung der Slater-Determinante erfolgt pro Spinsektor, da jeder Sektor seine eigene Form (Elektronen, Orbitale) hat. Ein einzelnes PrepareSlaterDeterminant-Gate (eine Besetzung, bei der die niedrigsten nocc Moden gefüllt sind, zusammen mit der Einteilchenrotation) bereitet einen Sektor vor; dasselbe Gate wird einmal auf die Alpha-Moden 0..norb und einmal auf die Beta-Moden norb..2*norb angewendet, da beide Spinsektoren hier dieselbe Referenzbesetzung und Rotation teilen:

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

Dieses deklarative Gate ersetzt das von Hand codierte XXPlusYYGate-Netzwerk; der Transpiler erzeugt die Synthese der Givens-Rotationen.

4. Die Krylov-Circuits zusammensetzen

Jeder Krylov-Circuit besteht aus der Vorbereitung des Referenzzustands, gefolgt von der Zeitentwicklung unter $e^{-i k \, \Delta t \, H}$. Teile den Hamiltonoperator in seinen Einteilchen- und Zweiteilchenanteil, $H = H_1 + H_2$, und trotterisiere die Zeitentwicklung über diese Aufteilung. Das ist die entscheidende Wahl für die Circuit-Tiefe. Die Einteilchen-Zeitentwicklung $e^{-i t H_1}$ ist eine Operation freier Fermionen (fermionisch gaußsch), also eine Einteilchen-Basisrotation (eine OrbitalRotation mit der Unitären $e^{-i t \, \mathbf{h}_1}$) und nicht etwas, das Term für Term Pauli-trotterisiert werden müsste. Die Zweiteilchen-Zeitentwicklung $e^{-i t H_2}$ wirkt nur auf den einzelnen On-Site-Wechselwirkungsterm. Wenn du den Schritt so aufbaust, kann der Transpiler ein flaches Brickwork aus Givens-Rotationen plus ein einzelnes Zwei-Qubit-Wechselwirkungs-Gate ausgeben, die flache Struktur, die die SKQD-Publikation von Hand aufbaut.

Den Zweiteilchenanteil haben wir bereits. Es ist der Operator two_body_momentum, der in Schritt 2 aufgebaut wurde (der einzelne On-Site-Term). Die Einteilchen-Zeitentwicklung $e^{-i t H_1}$ ist spinunabhängig, also eine einzelne OrbitalRotation pro Sektor mit der Unitären $e^{-i t \, \mathbf{h}_1}$, die einmal auf jeden Spinsektor angewendet wird, dieselbe parallele Struktur wie bei der Vorbereitung des Referenzzustands.

Jeder Krylov-Circuit entwickelt sich für eine Gesamtzeit $k \, \Delta t$. Das wird mit einem Trotter-Produkt zweiter Ordnung aus $k$ Schritten der Größe $\Delta t$ realisiert, sodass der Fehler pro Schritt konstant bleibt, wenn die Krylov-Dimension wächst. Jeder Schritt klemmt eine vollständige Einteilchenrotation (eine pro Spinsektor) zwischen zwei Halbschritte der (günstigen, einterm-)Zweiteilchen-Zeitentwicklung. Die Krylov-Dimension $D$ ist die Anzahl solcher Circuits, wobei dim von 0 (nur Referenzzustand) bis $D - 1$ reicht.

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

Hinweis

Die hier aufgebauten Circuits enthalten noch keine Messungen, weil dieselben fermionischen Circuits auch die exakte Simulation in Schritt 6 antreiben und eine Messung keine unitäre Operation ist. Schritt 5 fügt sie mit FermionicCircuit.measure_all() hinzu, das Mode i in das klassische Bit i ausliest.

Bevorzuge es, den fermionischen Circuit zu messen statt den transpilierten. Eine fermionische Messung benennt die Mode, die sie ausliest, sodass Layout und Routing sie mit dem Qubit mitführen, auf dem diese Mode gelandet ist, und der gesampelte Bitstring nach Moden indiziert ist. Wenn du stattdessen den transpilierten Circuit misst, musst du final_index_layout() selbst rückgängig machen, und du erhältst stillschweigend permutierte Ergebnisse, wenn du es vergisst. Das setzt die Jordan-Wigner-Kodierung in der Besetzungsbasis voraus, die das Preset verwendet und bei der jede Mode von einem Qubit getragen wird.

5. In Qubit-Circuits transpilieren

Das Transpilieren eines FermionicCircuit bildet seine fermionischen Gates auf Qubit-Gates ab. Dabei wird die Jordan-Wigner-Abbildung jordan_wigner() verwendet. Die OrbitalRotation-Gates werden zu einem Brickwork aus Givens-Rotationen synthetisiert, und die einterm-Evolution wird zu einem einzelnen Zwei-Qubit-Wechselwirkungs-Gate heruntergebrochen.

Da jedes fermionische Gate hier eine Standardsynthese hat, kannst du das fertige generate_preset_jw_pass_manager() direkt verwenden, statt die Stufen von Hand zusammenzusetzen. Alle Schlüsselwortargumente werden an die Qubit-Stufe weitergereicht; übergib optimization_level=0 für ein unverfälschtes, nicht optimiertes Bild der synthetisierten Tiefe. Die gesamte Familie der Krylov-Circuits zu erzeugen ist dann eine Schleife über die Krylov-Dimension. Behalte die nicht transpilierten FermionicCircuits neben den transpilierten Qubit-Circuits. Die fermionischen treiben die exakte (Statevector-)Simulation in Schritt 6 an und müssen daher messungsfrei bleiben; measure_all(inplace=False) gibt eine gemessene Kopie zurück, die transpiliert wird und das Original unberührt lässt, und diese Kopie würde ein Backend ausführen.

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

Die Krylov-Dimension zeigt sich direkt in der Zwei-Qubit-Tiefe der Circuits. Die Vorbereitung des Referenzzustands ist ein fester Aufwand, und jede zusätzliche Krylov-Potenz fügt einen weiteren Trotter-Schritt zweiter Ordnung des Zeitentwicklungsoperators hinzu.

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

Hinweis

Die Aufteilung des Hamiltonoperators hält diese Tiefen moderat. Die Einteilchen-Zeitentwicklung ist eine echte Orbitalrotation und wird daher zu einem XXPlusYYGate-Brickwork auf benachbarten Qubits synthetisiert. Jeder Halbschritt der On-Site-Wechselwirkung ist ein einzelnes RZZGate, das die beiden Spinmoden der Störstelle koppelt; der günstige Term, den das Produkt zweiter Ordnung um die Rotation klemmt. Würde man den gesamten Hamiltonoperator stattdessen an ein einziges Evolution-Gate übergeben, würde der Einteilchenanteil Term für Term Pauli-trotterisiert, wobei sich lange Jordan-Wigner-Z-Ketten über die (langreichweitigen) Kopplungen zwischen Störstelle und Bad ausbreiten und die Tiefe um eine Größenordnung aufblähen.

Der erste Circuit bereitet nur die Referenzdeterminante vor, während der letzte die meisten Trotter-Schritte trägt (hier vier):

>>> 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. Bitstrings aus den Krylov-Circuits sampeln

Auf Hardware würdest du die transpilierten Circuits ausführen und messen. Dies ist die rauschfreie (Statevector-)Version dieses Schritts, die die Bitstrings zeigt, die SKQD sammeln würde. Simuliere jeden nicht transpilierten FermionicCircuit mit ffsim.apply_unitary() und sample den resultierenden Statevector mit ffsim.sample_state_vector().

Ein PrepareSlaterDeterminant-Gate validiert in der Simulation zuerst und rotiert dann (es prüft, dass der eingehende Zustand seine Referenzdeterminante besetzt, und wendet dann die Rotation an). Daher ist die Eingabe die einfache Besetzungsdeterminante: die niedrigsten nocc Moden pro Spin gefüllt, ohne Rotation. Das Gate wendet reference_rotation selbst an.

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

Fasst man die Samples aus allen fünf Krylov-Circuits zusammen, ist der Träger winzig: einige hundert verschiedene Bitstrings aus den $\binom{8}{4}^2 = 4900$ Determinanten im Sektor (4, 4). Diese Konzentration, eine direkte Folge der in Schritt 2 festgestellten Dünnbesetztheit in der Impulsbasis, macht die anschließende klassische Diagonalisierung handhabbar.

Das Darstellen der Counts mit qiskit.visualization.plot_histogram() macht die Sparsity sichtbar. Wenige Konfigurationen dominieren, angeführt von der Referenzdeterminante.

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

7. In dem gesampelten Unterraum diagonalisieren

Die gesampelten Bitstrings sind die Eingabe für die klassische Hälfte von SKQD. Der Hamiltonoperator wird auf den Unterraum projiziert, den die gesampelten Konfigurationen aufspannen, und dort diagonalisiert. Dies wird an das Paket qiskit-addon-sqd übergeben, dessen diagonalize_fermionic_hamiltonian() die vollständige Schleife der Sample-based Quantum Diagonalization ausführt. Sie baut aus Batches der gesampelten Konfigurationen einen Unterraum auf, diagonalisiert den Hamiltonoperator darin und verfeinert den Unterraum iterativ mithilfe der Konfigurationswiederherstellung: Dabei werden Besetzungen umgekehrt, um Konfigurationen zu reparieren, die die bekannte Teilchenzahlsymmetrie verletzen. Die Konfigurationswiederherstellung soll die Schäden durch Hardware-Rauschen beheben, das gesampelte Bitstrings in den falschen Teilchenzahlsektor verfälscht; das Sampling ist hier rauschfrei, daher gibt es in dieser Hinsicht nichts zu reparieren.

Der Solver verarbeitet die gesampelten Counts als BitArray, das direkt aus den zusammengeführten counts von Schritt 6 aufgebaut wird. Übergib num_bits explizit: Ohne diese Angabe wird die Breite aus dem breitesten gesampelten Bitstring abgeleitet, der schmaler als 2 * norb ist, sobald keine Probe den letzten Modus besetzt, und der Solver lehnt das Array dann ab. Der Solver diagonalisiert in derselben Impulsbasis, die auch der Rest der Konstruktion verwendet, daher nimmt er die Integrale h1e_momentum und h2e_momentum aus Schritt 2 direkt, ohne Umpacken. Die wiederhergestellte Energie stimmt aus nur wenigen hundert gesampelten Konfigurationen auf deutlich weniger als ein Milli-Hartree mit der exakten Referenz überein. Da die SQD-Energie eine variationelle obere Schranke für den wahren Grundzustand ist, liegt sie immer bei oder über reference_energy (weiterhin gültig aus Schritt 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)

Visualisiere den Lauf: die Konvergenz der Energie über die Iterationen und die mittlere Besetzung jedes räumlichen Orbitals im wiederhergestellten Grundzustand. Die Energie pro Iteration ist die niedrigste über die Batches dieser Iteration, und die Besetzung summiert die orbital_occupancies über die beiden Spinsektoren. Auch ohne zu korrigierendes Rauschen sinkt die Energie von Iteration zu Iteration. Jede Iteration diagonalisiert innerhalb von num_batches Batches mit jeweils nur samples_per_batch aus dem Pool unterabgetasteten Konfigurationen, und die iterative Verfeinerung durch Bitstring-Übernahme verbessert stetig, welche Konfigurationen in diesen Batches landen, sodass sich die wiederhergestellte Energie verbessert:

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

Damit ist die SKQD-Schleife geschlossen. Die Impulsbasis-Circuits der Schritte 1–5 erzeugen ein dünnbesetztes Sample (Schritt 6), und die obige Sample-basierte Diagonalisierung macht aus diesem Sample wieder eine Grundzustandsenergie.

Nächste Schritte

Dieser Leitfaden hat die gesamte SKQD-Pipeline auf einem rauschfreien Statevector-Simulator ausgeführt. Auf Hardware würden die transpilierten circuits anstelle des Statevector-Samplings in Schritt 6 ausgeführt und gemessen, wobei die resultierenden Counts unverändert in die Diagonalisierung von Schritt 7 einfließen. Hilfe zum Ausführen von Circuits findest du in der Qiskit-Dokumentation, und Informationen zur Nachverarbeitung durch Unterraum-Diagonalisierung, einschließlich ihres Verhaltens bei verrauschten Samples, bei denen die Konfigurationswiederherstellung am meisten leistet, findest du in den SQD-Addon-Tutorials.

Lies den ffsim-Beziehungsleitfaden, um zu verstehen, warum ffsim.apply_unitary() und ffsim.sample_state_vector() nativ mit den fermionischen Circuits dieses Pakets funktionieren und warum der in Schritt 6 verwendete Sektor mit fester Teilchenzahl die in Schritt 2 festgestellte Sparsity der Impulsbasis sample-effizient macht.