Krylov-Zeitentwicklungsschaltungen (SKQD) generieren

Die samplebasierte Krylov-Quanten-Diagonalisierung ( SKQD ) ist eine quantenzentrierte Variante der Unterraum-Diagonalisierung. Anstatt einen einzelnen Variationsansatz (wie bei der einfachen, auf Stichproben basierenden Quantendiagonalisierung ( SQD )) oder ein Ensemble aus zufälligen Zeitentwicklungsschaltungen (wie bei „ SqDRIFT “) zu untersuchen, untersucht es eine Krylov-Basis : die Familie der Zustände

erhalten durch die Entwicklung eines Referenzzustands $\lvert \psi_0 \rangle$ über steigende Vielfache eines festen Zeitschritts $\Delta t$. Aus diesen $D$ Schaltungen abgetastete Bitfolgen füllen den Krylov-Unterraum, in dem der Hamiltonoperator anschließend klassisch diagonalisiert wird.

Diese Anleitung gibt die Schritte zum Aufbau der Schaltung für den SKQD-Algorithmus für das Anderson-Modell mit einer einzelnen Verunreinigung (SIAM) wieder, bei dem es sich um eine magnetische Verunreinigung handelt, die an ein nicht wechselwirkendes Bad gekoppelt ist, und folgt dabei der SKQD-Veröffentlichung.

1. Der SIAM-Hamilton-Operator als fermionischer Operator

Das Anderson-Modell mit einer einzigen Verunreinigung sieht ein wechselwirkendes Verunreinigungsorbital (Index „ $0$ “) vor, das mit einem nicht wechselwirkenden Bad aus „ $L$ “-Stellen in Kontakt steht, die als Kette angeordnet sind. Sein Hamilton-Operator zerlegt sich in einen Ein-Körper-Anteil $H_1$ und einen Zwei-Körper-Anteil $H_2$,

wobei das Ein-Körper-Modell die Energie der Verunreinigung vor Ort (das chemische Potential $\mu$ ) erfasst, die Hybridisierung zwischen Verunreinigung und Bad $V$ und das Hopping zwischen den nächsten Nachbarn im Bad $t$,

und der Zwei-Körper-Teil, der die vor Ort auftretende Coulomb-Abstoßung $U$ allein auf die Verunreinigung bezieht,

Dieser Leitfaden basiert auf dem Partikel-Loch-symmetrischen Punkt $\mu = -U/2$. Diese beiden Teile lassen sich direkt auf die Ein-Körper-Integrale ( h1e Bath-Hopping, Hybridisierung der Verunreinigung, On-Site-Energie der Verunreinigung) und das einzige Zwei-Körper-Integral h2e (den Coulomb-Term der Verunreinigung) abbilden. Setzen Sie sie zunächst in der Positionsbasis zusammen, in der das Modell naturgemäß 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. Umstellung auf die Impulsbasis

Der Bad-Teil des SIAM ist eine Kette aus freien Fermionen, sodass er durch eine Umstellung auf die Impulsbasis diagonalisiert wird. Dies ist der entscheidende Schritt für SKQD. Da das Bad translationsinvariant ist, ist der Grundzustand in der Positionsbasis nicht spärlich, wohl aber in der Impulsbasis, in der der Bad-Hamilton-Operator diagonal ist. Die Orbitalrotation eines einzelnen Teilchens, die diese Veränderung bewirkt, ist problemspezifische lineare Algebra. Dabei wird die Bath-Hopping-Matrix diagonalisiert, während das Verunreinigungsorbital unverändert bleibt, und anschließend wird die Verunreinigung an eine zentrale Stelle verlagert:

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

Im Gegensatz zu einem variationalen Arbeitsablauf werden hier die Hamilton-Integrale in die Impulsbasis transformiert, sodass sich der gesamte Schaltkreis (sowohl der Hamilton-Operator als auch der Zustand) in der Basis befindet, in der der Grundzustand spärlich ist und die abgetasteten Bitfolgen daher aussagekräftig sind. Die Rotation der Ein- und Zweikörperintegrale stellt eine Standard-Basisänderung der Elektronenstruktur-Tensoren dar:

>>> 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 Verunreinigungsorbital unverändert lässt (sie wirkt sich nur auf das Bad aus), bleibt die On-Site-Wechselwirkung unverändert. h2e_momentum bleibt ein Term mit nur einer Zahl, genau wie h2e. Die einzige Änderung besteht darin, dass die Unreinheit in einen zentralen Index verschoben wurde. Dies lässt sich direkt an den Integraltensoren erkennen, die jeweils einen Eintrag ungleich Null aufweisen:

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

Dank dieser Lokalisierung lässt sich die Interaktion später kostengünstig synthetisieren.

Fassen Sie die Ein-Körper-Integrale FermionOperator mithilfe des Elektronenintegral-Konstruktors zu einem zusammen from_1body_tril_spin_sym(), der erwartet, dass diese in einer unterhalb dreieckigen Reihenfolge angeordnet sind. Der Zwei-Körper-Teil ist ein einzelner Zahlen-Zahlen-Term; daher sollte man ihn nicht über den Mechanismus der Zwei-Körper-Integrale leiten, sondern direkt mit cre()/ bildenann(): dem Coulomb-Operator vor Ort $U \, n_{p\uparrow} n_{p\downarrow}$ für den Verunreinigungsmodus $p$ (sein $\uparrow$ -Modus am Index p, sein $\downarrow$ -Modus am Index p + norb im Block-Spin-Layout). Beide Konstruktoren wandeln die räumlichen norb Orbitalen spin-doppelt in 2 * norb fermionische Moden um.

Behalten Sie den Ein-Körper-Teil $H_1$ und den Zwei-Körper-Teil $H_2$ als separate Operatoren bei, anstatt sie vorschnell zu summieren, entsprechend den Termen $H_1$ und $H_2$ der section-1 -Gleichungen. Wenn man jedes Teil einmal (in beiden Basen) baut, kann man es weiter unten wiederverwenden. Zusammengerechnet ergeben sie den vollständigen Hamilton-Operator für die exakte Diagonalisierung in beiden Basen. Im Einzelnen ist $H_2$ (auf der Impulsbasis) der Operator, in 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)
... )

Die Basisänderung ist eine unitäre Ähnlichkeitstransformation, sodass das Spektrum davon unberührt bleibt. Überprüfen Sie, ob der Operator korrekt ist, indem Sie seinen kleinsten Eigenwert im Sektor „ $(\text{norb}, \text{nelec})$ “ mit dem Ergebnis einer exakten Diagonalisierung vergleichen. Die Funktion stellt FermionOperator eine „ SciPy “ bereit ( LinearOperator die auf einem nativen Matrix-Vektor-Kernel für die vollständige Konfigurationsinteraktion (FCI) basiert), sodass sie direkt an übergeben werden kann 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

Der Vorteil der Basisänderung liegt in der Sparsamkeit. Um dies zu veranschaulichen, zählen Sie, wie viele Determinanten der Berechnungsbasis erforderlich sind, um 99 % des Grundzustandsgewichts zu erfassen, und vergleichen Sie die Impulsbasis mit der Positionsbasis. Die Energie ist basisunabhängig, die Anzahl der Determinanten jedoch 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 zwar dieselbe Energie, doch der Grundzustand in der Positionsbasis verteilt sich auf Tausende von Determinanten, während er sich in der Impulsbasis auf etwa 20 konzentriert; deshalb funktioniert das SKQD-Sampling in der Impulsbasis und nicht in der Positionsbasis.

3. Den Referenzzustand vorbereiten

Der SKQD-Referenz $\lvert \psi_0 \rangle$ us ist die Überlagerung aller Anregungen der Elektronen, die dem Fermi-Niveau am nächsten liegen, in die nächstgelegenen leeren Impulsmoden. Die SKQD -Publikation baut darauf ein handgeschriebenes Netzwerk aus XXPlusYYGates auf. Es scheint stark korreliert zu sein (es weist Hunderte von Amplituden der Berechnungsbasis ungleich Null auf), doch ein Netzwerk aus zahlenkonservierenden Zwei-Modus-Rotationen ist eine Ein-Teilchen-Operation (fermionische Gauß-Operation), sodass der dadurch erzeugte Zustand eine einzelne Slater-Determinante ist. Das bedeutet, dass Sie es mit einem einzigen Gate PrepareSlaterDeterminant vorbereiten können, da das Netzwerk die Rotation einzelner Teilchen implementiert.

Da jedes davon als Givens-Rotation auf ein Modenpaar XXPlusYYGate wirkt, lässt sich diese Rotation direkt im Einzelteilchenmodell aufbauen (eine Givens-Matrix pro Gatter), anstatt die Zwei-Qubit-Gatter manuell zu generieren. Die Gitter bilden ein Geflecht aus „Nearest-Neighbor“-Rotationen, das sich vom Fermi-Niveau aus in sich ausdehnenden Schichten auffächert und jedes Elektron in der Nähe des Fermi-Niveaus 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])

Die Berechnung der Slater-Determinante erfolgt pro Spinsektor, da jeder Sektor seine eigene Form (Elektronen, Orbitale) aufweist. Ein einzelnes Gate PrepareSlaterDeterminant (eine Besetzung, die gefülltsten niedrigsten nocc Moden sowie die Rotation eines einzelnen Teilchens) bereitet einen Sektor vor; dasselbe Gate wird einmal auf die Alpha-Moden und 0..norb einmal auf die Beta-Moden angewendet norb..2*norb, da beide Spin-Sektoren hier dieselbe Referenzbesetzung und -rotation aufweisen:

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

Dieses deklarative Gatter ersetzt das manuell programmierte XXPlusYYGate Netzwerk; der Transpiler generiert die Givens-Rotationssynthese.

4. Bauen Sie die Krylov-Schaltungen zusammen

Jede Krylov-Schaltung besteht aus der Vorbereitung des Referenzzustands, gefolgt von einer Entwicklung gemäß $e^{-i k \, \Delta t \, H}$. Man zerlegt den Hamiltonoperator in seinen Ein- und Zweikörperteil, $H = H_1 + H_2$, und wendet die Trotter-Transformation auf die Entwicklung über diese Zerlegung an. Dies ist die entscheidende Wahl für die Schaltungstiefe. Die Ein-Teilchen-Evolution $e^{-i t H_1}$ ist eine Operation mit freien Fermionen (fermionische Gauß-Verteilung), es handelt sich also um eine Basisrotation für ein einzelnes Teilchen (eine OrbitalRotation mit unitären $e^{-i t \, \mathbf{h}_1}$ en) und nicht um etwas, das Term für Term nach Pauli-Trotter transformiert werden muss. Die Zwei-Körper-Entwicklungs $e^{-i t H_2}$ e wirkt sich nur auf den einzelnen On-Site-Wechselwirkungsterm aus. Wenn man diesen Schritt auf diese Weise aufbaut, kann der Transpiler ein flaches Geflecht aus Givens-Rotationen sowie ein einzelnes 2-Qubit-Wechselwirkungsgatter erzeugen – genau jene flache Struktur, die in der SKQD- Veröffentlichung manuell aufgebaut wird.

Den Teil über das Zwei-Körper-Problem haben wir bereits. Es handelt sich um den in Schritt 2 gebildeten Operator two_body_momentum (den einzelnen Vor-Ort-Term). Die Ein-Teilchen-Entwicklungs $e^{-i t H_1}$ t spinunabhängig, es handelt sich also um eine einzige Operation pro Sektor, OrbitalRotation bei der die unitäre $e^{-i t \, \mathbf{h}_1}$ einmal auf jeden Spinsektor angewendet wird – dieselbe parallele Struktur wie bei der Vorbereitung des Referenzzustands.

Jede Krylov-Iteration läuft über eine Gesamtzeit von $k \, \Delta t$. Dies wird durch ein Trotter-Produkt zweiter Ordnung aus $k$ Schritten der Größe $\Delta t$ realisiert, sodass der Fehler pro Schritt konstant bleibt, während die Krylov-Dimension zunimmt. Jeder Schritt umfasst eine vollständige Ein-Körper-Rotation (eine pro Spin-Sektor) zwischen zwei Halbschritten der (einfachen, aus einem Term bestehenden) Zwei-Körper-Entwicklung. Die Krylov-Dimension $D$ ist die Anzahl solcher Schaltkreise, wobei im Bereich von 0 (nur Referenzzustand) bis $D - 1$ liegt dim .

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

5. In Qubit-Schaltungen transpilieren

Beim Transpiling werden die fermionischen Gatter auf Qubit-Gatter abgebildet FermionicCircuit . Hierbei wird die Jordan-Wigner-Abbildung jordan_wigner() verwendet. Die Gatter OrbitalRotation lassen sich zu einem Verbund aus Givens-Rotationen zusammenfassen, und der Ein-Term-Ausdruck Evolution lässt sich auf ein einziges Zwei-Qubit-Wechselwirkungsgatter reduzieren.

Da jedes Fermionen-Gatter hier über eine Standardsynthese verfügt, können Sie die vorgefertigten Module direkt generate_preset_jw_pass_manager() verwenden, anstatt die Stufen von Hand zusammenzusetzen. Alle Schlüsselwortargumente werden an die Qubit-Stufe weitergeleitet; verwenden Sie „pass“, um optimization_level=0 ein originalgetreues, nicht optimiertes Bild der synthetisierten Tiefe zu erhalten. Die Erzeugung der gesamten Familie von Krylov-Schaltungen erfolgt dann mittels einer Schleife über die Krylov-Dimension. Bewahren Sie die nicht transpilierten FermionicCircuits neben den transpilierten Qubit-Schaltkreisen auf. Die fermionischen Varianten dienen als Grundlage für die exakte (Zustandsvektor-)Simulation in Schritt 6, während die transpilierten Varianten (mit hinzugefügten Messungen) dem entsprechen, was ein Backend ausführen würde.

>>> 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)
...         circuit = pass_manager.run(fermionic_circuit)
...         circuit.measure_all()
...         circuits.append(circuit)
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

Die Krylov-Dimension spiegelt sich direkt in der Zwei-Qubit-Tiefe der Schaltungen wider. Die Vorbereitung des Referenzzustands ist mit festen Kosten verbunden, und jede zusätzliche Krylov-Potenz führt zu einem weiteren Trotter-Schritt zweiter Ordnung des Zeitentwicklungsoperators.

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

Die erste Schleife berechnet lediglich die Referenzdeterminante, während die letzte die meisten Trotter-Schritte (in diesem Fall vier) umfasst:

>>> 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. Beispiel-Bitstrings aus den Krylov-Schaltungen

Auf der Hardware würde man die transpilierten Schaltungen ausführen und Messungen durchführen. Dies ist die rauschfreie (Zustandsvektor-)Version dieses Schritts, die die von SKQD erfassten Bitfolgen zeigt. Simuliere jedes nicht transpilierte Element mit FermionicCircuit und ffsim.apply_unitary() erzeuge eine Stichprobe des resultierenden Zustandsvektors mit ffsim.sample_state_vector().

Ein Gate PrepareSlaterDeterminant arbeitet in der Simulation nach dem Prinzip „Validieren und dann rotieren“ (es prüft, ob der eingehende Zustand seine Referenzdeterminante einnimmt, und wendet dann die Rotation an). Daher ist die Eingabe die einfache Besetzungsdeterminante: die niedrigsten pro Spin besetzten Moden nocc , ohne Rotation. Das Tor öffnet sich reference_rotation von selbst.

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

Führt man die Stichproben aus allen fünf Krylov-Schaltungen zusammen, ist die Unterstützungsmenge winzig: nur wenige hundert verschiedene Bitstrings von den insgesamt $\binom{8}{4}^2 = 4900$ Determinanten in diesem (4, 4) Sektor. Diese Konzentration, eine direkte Folge der in Schritt 2 festgestellten Sparsität auf der Impulsbasis, macht die anschließende klassische Diagonalisierung überschaubar.

Wenn man die Zählwerte mit darstellt, qiskit.visualization.plot_histogram() wird die Sparsamkeit sichtbar. Es gibt einige wenige Konfigurationen, die 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. Diagonalisierung im abgetasteten Unterraum

Die abgetasteten Bitstrings bilden die Eingabe für den klassischen Teil von SKQD. Der Hamilton-Operator wird auf den Unterraum projiziert, den die abgetasteten Konfigurationen aufspannen, und dort diagonalisiert. Dies wird an das Paket „qiskit-addon-sqd“ übergeben, das die gesamte beispiellbasierte Quanten-Diagonalisierungsschleife diagonalize_fermionic_hamiltonian() ausführt. Es wird ein Unterraum aus Chargen der abgetasteten Konfigurationen gebildet, der Hamiltonoperator darin diagonalisiert und der Unterraum durch Konfigurationswiederherstellung iterativ verfeinert: Dabei werden Besetzungen umgekehrt, um Konfigurationen zu korrigieren, die gegen die bekannte Teilchenzahlsymmetrie verstoßen. Die Konfigurationswiederherstellung dient dazu, die durch Hardware-Rauschen verursachten Schäden rückgängig zu machen, die dazu führen, dass abgetastete Bitfolgen in den Sektor mit der falschen Teilchenanzahl verschoben werden; die Abtastung erfolgt hier rauschfrei, sodass in dieser Hinsicht nichts zu reparieren ist.

Der Solver verwendet die Stichprobenzahlen als BitArray, die direkt aus dem zusammengefassten counts aus Schritt 6 gebildet werden. Es diagonalisiert in derselben Impulsbasis, die auch für den Rest der Konstruktion verwendet wird, sodass die Integrale h2e_momentum h1e_momentum und aus Schritt 2 direkt übernommen werden können, ohne dass eine Umformulierung erforderlich ist. Die zurückgewonnene Energie stimmt bis auf weit unter ein Milli-Hartree genau mit dem Referenzwert überein, und das bei nur wenigen hundert untersuchten Konfigurationen. Da die SQD-Energie eine variationale Obergrenze für den tatsächlichen Grundzustand darstellt, liegt sie immer auf oder über reference_energy (was weiterhin aus Schritt 2 gilt):

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts))
>>>
>>> # 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)

Stellen Sie sich den Durchlauf vor: die Konvergenz der Energie über die Iterationen hinweg und die durchschnittliche Besetzung jedes Raumorbital im wiederhergestellten Grundzustand. Die Energie pro Iteration ist die niedrigste aller Batches dieser Iteration, und die Belegung summiert sich über orbital_occupancies die beiden Spinsektoren. Selbst wenn kein Rauschen zu korrigieren ist, nimmt die Energie von Iteration zu Iteration ab. Bei jeder Iteration erfolgt die Diagonalisierung innerhalb von num_batches Chargen mit jeweils nur samples_per_batch Konfigurationen, die aus dem Pool als Stichprobe entnommen werden, und durch die iterative Verfeinerung aufgrund des Bitfolge-Übertrags verbessert sich stetig, welche Konfigurationen in diesen Chargen landen, sodass sich die zurückgewonnene 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 der SKQD-Kreislauf abgeschlossen. Die auf dem Impuls basierenden Schaltkreise der Schritte 1–5 erzeugen eine spärliche Stichprobe (Schritt 6), und die oben beschriebene stichprobenbasierte Diagonalisierung wandelt diese Stichprobe wieder in eine Grundzustandsenergie um.

Nächste Schritte

In dieser Anleitung wurde die gesamte SKQD-Pipeline auf einem rauschfreien Zustandsvektor-Simulator ausgeführt. Auf der Hardware würde circuits der transpilierte Code anstelle der Zustandsvektor-Abtastung in Schritt 6 ausgeführt und gemessen werden, wobei die resultierenden Zählwerte unverändert in die Diagonalisierung nach step-7 einfließen würden. Hilfe zur Ausführung von Schaltkreisen finden Sie in der Qiskit-Dokumentation, und Informationen zur Nachbearbeitung durch Subraum-Diagonalisierung – einschließlich ihres Verhaltens bei verrauschten Stichproben, bei denen die Konfigurationswiederherstellung den größten Teil der Arbeit leistet – finden Sie in den Tutorials zum SQD-Add-on.

Lesen Sie den ffsim-Beziehungsleitfaden, um zu verstehen, warum und ffsim.apply_unitary() auf den fermionischen Schaltkreisen dieses Pakets nativ ffsim.sample_state_vector() funktionieren und warum der in Schritt 6 verwendete Sektor mit fester Teilchenzahl die in Schritt 2 festgelegte Sparsamkeit der Impulsbasis sample-effizient macht.