Stichprobenbasierte Krylov-Quantendiagonalisierung eines fermionischen Gittermodells

Nutzungsschätzung: Neun Sekunden auf einem Heron r2 Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann variieren.)

Hintergrund

Dieses Tutorial zeigt, wie du die stichprobenbasierte Quantendiagonalisierung (SQD) verwendest, um die Grundzustandsenergie eines fermionischen Gittermodells zu schätzen. Konkret untersuchen wir das eindimensionale Einzelstörstellen-Anderson-Modell (SIAM), das zur Beschreibung magnetischer Störstellen in Metallen verwendet wird.

Dieses Tutorial folgt einem ähnlichen Arbeitsablauf wie das verwandte Tutorial Stichprobenbasierte Quantendiagonalisierung eines Chemie-Hamiltonians. Ein wesentlicher Unterschied liegt jedoch darin, wie die Quantenschaltkreise aufgebaut werden. Das andere Tutorial verwendet einen heuristischen Variationsansatz, der für Chemie-Hamiltonians mit potenziell Millionen von Wechselwirkungstermen attraktiv ist. Dieses Tutorial hingegen verwendet Schaltkreise, die die Zeitentwicklung durch den Hamiltonian approximieren. Solche Schaltkreise können tief sein, was diesen Ansatz besser für Anwendungen auf Gittermodelle macht. Die von diesen Schaltkreisen präparierten Zustandsvektoren bilden die Basis für einen Krylov-Unterraum, und als Ergebnis konvergiert der Algorithmus nachweislich und effizient zum Grundzustand unter geeigneten Annahmen.

Der in diesem Tutorial verwendete Ansatz kann als eine Kombination der in SQD und Krylov-Quantendiagonalisierung (KQD) verwendeten Techniken betrachtet werden. Der kombinierte Ansatz wird manchmal als stichprobenbasierte Krylov-Quantendiagonalisierung (SQKD) bezeichnet. Siehe Krylov-Quantendiagonalisierung von Gitter-Hamiltonians für ein Tutorial zur KQD-Methode.

Dieses Tutorial basiert auf der Arbeit "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", auf die für weitere Details verwiesen werden kann.

Einzelstörstellen-Anderson-Modell (SIAM)

Der eindimensionale SIAM-Hamiltonian ist eine Summe aus drei Termen:

$$H = H_{\textrm{imp}}+ H_\textrm{bath} + H_\textrm{hyb},$$

wobei

$$\begin{align*} H_\textrm{imp} &= \varepsilon \left( \hat{n}_{d\uparrow} + \hat{n}_{d\downarrow} \right) + U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}, \\ H_\textrm{bath} &= -t \sum_{\substack{\mathbf{j} = 0\\ \sigma\in \{\uparrow, \downarrow\}}}^{L-1} \left(\hat{c}^\dagger_{\mathbf{j}, \sigma}\hat{c}_{\mathbf{j}+1, \sigma} + \hat{c}^\dagger_{\mathbf{j}+1, \sigma}\hat{c}_{\mathbf{j}, \sigma} \right), \\ H_\textrm{hyb} &= V\sum_{\sigma \in \{\uparrow, \downarrow \}} \left(\hat{d}^\dagger_\sigma \hat{c}_{0, \sigma} + \hat{c}^\dagger_{0, \sigma} \hat{d}_{\sigma} \right). \end{align*}$$

Hier sind $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ die fermionischen Erzeugungs-/Vernichtungsoperatoren für die $\mathbf{j}^{\textrm{te}}$ Bad-Stelle mit Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ sind Erzeugungs-/Vernichtungsoperatoren für den Störstellenmodus, und $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ und $V$ sind reelle Zahlen, die die Hüpf-, Vor-Ort- und Hybridisierungswechselwirkungen beschreiben, und $\varepsilon$ ist eine reelle Zahl, die das chemische Potenzial spezifiziert.

Beachte, dass der Hamiltonian eine spezifische Instanz des generischen Wechselwirkungs-Elektronen-Hamiltonians ist,

$$\begin{align*} H &= \sum_{\substack{p, q \\ \sigma}} h_{pq} \hat{a}^\dagger_{p\sigma} \hat{a}_{q\sigma} + \sum_{\substack{p, q, r, s \\ \sigma \tau}} \frac{h_{pqrs}}{2} \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma} \\ &= H_1 + H_2, \end{align*}$$

wobei $H_1$ aus Einkörpertermen besteht, die quadratisch in den fermionischen Erzeugungs- und Vernichtungsoperatoren sind, und $H_2$ aus Zweikörpertermen besteht, die quartisch sind. Für das SIAM gilt

$$H_2 = U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}$$

und $H_1$ enthält die restlichen Terme im Hamiltonian. Um den Hamiltonian programmatisch darzustellen, speichern wir die Matrix $h_{pq}$ und den Tensor $h_{pqrs}$.

Orts- und Impulsbasen

Aufgrund der approximativen Translationssymmetrie in $H_\textrm{bath}$ erwarten wir nicht, dass der Grundzustand in der Ortsbasis (der Orbitalbasis, in der der Hamiltonian oben spezifiziert ist) dünn besetzt ist. Die Leistung von SQD ist nur garantiert, wenn der Grundzustand dünn besetzt ist, das heißt, wenn er signifikantes Gewicht auf nur einer kleinen Anzahl von Rechenbasis-Zuständen hat. Um die Dünnbesetztheit des Grundzustands zu verbessern, führen wir die Simulation in der Orbitalbasis durch, in der $H_\textrm{bath}$ diagonal ist. Wir nennen diese Basis die Impulsbasis. Da $H_\textrm{bath}$ ein quadratischer fermionischer Hamiltonian ist, kann er effizient durch eine Orbitalrotation diagonalisiert werden.

Approximative Zeitentwicklung durch den Hamiltonian

Um die Zeitentwicklung durch den Hamiltonian zu approximieren, verwenden wir eine Trotter-Suzuki-Zerlegung zweiter Ordnung,

$$e^{-i \Delta t H} \approx e^{-i\frac{\Delta t}{2} H_2} e^{-i\Delta t H_1} e^{-i\frac{\Delta t}{2} H_2}.$$

Unter der Jordan-Wigner-Transformation entspricht die Zeitentwicklung durch $H_2$ einem einzelnen CPhase-Gate zwischen den Spin-up- und Spin-down-Orbitalen an der Störstelle. Da $H_1$ ein quadratischer fermionischer Hamiltonian ist, entspricht die Zeitentwicklung durch $H_1$ einer Orbitalrotation.

Die Krylov-Basiszustände $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, wobei $D$ die Dimension des Krylov-Unterraums ist, werden durch wiederholte Anwendung eines einzelnen Trotter-Schritts gebildet, sodass

$$|\psi_k\rangle \approx \left[e^{-i\frac{\Delta t}{2} H_2} e^{-i\Delta t H_1} e^{-i\frac{\Delta t}{2} H_2} \right]^k\ket{\psi_0}.$$

Im folgenden SQD-basierten Arbeitsablauf werden wir Stichproben aus diesem Satz von Schaltkreisen ziehen und den kombinierten Satz von Bitfolgen mit SQD nachbearbeiten. Dieser Ansatz steht im Gegensatz zu dem im verwandten Tutorial Stichprobenbasierte Quantendiagonalisierung eines Chemie-Hamiltonians verwendeten, wo Stichproben aus einem einzelnen heuristischen Variationsschaltkreis gezogen wurden.

Anforderungen

Stelle vor Beginn dieses Tutorials sicher, dass du Folgendes installiert hast:

• Qiskit SDK v1.0 oder höher, mit Unterstützung für Visualisierung
• Qiskit Runtime v0.22 oder höher (pip install qiskit-ibm-runtime)
• SQD Qiskit Addon v0.11 oder höher (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Schritt 1: Problem auf einen Quantenschaltkreis abbilden

Zunächst erzeugen wir den SIAM-Hamiltonian in der Ortsbasis. Der Hamiltonian wird durch die Matrix $h_{pq}$ und den Tensor $h_{pqrs}$ dargestellt. Anschließend rotieren wir ihn in die Impulsbasis. In der Ortsbasis platzieren wir die Störstelle an der ersten Stelle. Wenn wir jedoch zur Impulsbasis rotieren, verschieben wir die Störstelle an eine zentrale Stelle, um Wechselwirkungen mit anderen Orbitalen zu erleichtern.

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy qiskit qiskit-addon-sqd qiskit-ibm-runtime scipy

import numpy as np

def siam_hamiltonian(
    norb: int,
    hopping: float,
    onsite: float,
    hybridization: float,
    chemical_potential: float,
) -> tuple[np.ndarray, np.ndarray]:
    """Hamiltonian for the single-impurity Anderson model."""
    # Place the impurity on the first site
    impurity_orb = 0

    # One body matrix elements in the "position" basis
    h1e = np.zeros((norb, norb))
    np.fill_diagonal(h1e[:, 1:], -hopping)
    np.fill_diagonal(h1e[1:, :], -hopping)
    h1e[impurity_orb, impurity_orb + 1] = -hybridization
    h1e[impurity_orb + 1, impurity_orb] = -hybridization
    h1e[impurity_orb, impurity_orb] = chemical_potential

    # Two body matrix elements in the "position" basis
    h2e = np.zeros((norb, norb, norb, norb))
    h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite

    return h1e, h2e

def momentum_basis(norb: int) -> np.ndarray:
    """Get the orbital rotation to change from the position to the momentum basis."""
    n_bath = norb - 1

    # Orbital rotation that diagonalizes the bath (non-interacting system)
    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)

    # Expand to include impurity
    orbital_rotation = np.zeros((norb, norb))
    # Impurity is on the first site
    orbital_rotation[0, 0] = 1
    orbital_rotation[1:, 1:] = vecs

    # Move the impurity to the center
    new_index = n_bath // 2
    perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]
    orbital_rotation = orbital_rotation[:, perm]

    return orbital_rotation

def rotated(
    h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    """Rotate the orbital basis of a Hamiltonian."""
    h1e_rotated = np.einsum(
        "ab,Aa,Bb->AB",
        h1e,
        orbital_rotation,
        orbital_rotation.conj(),
        optimize="greedy",
    )
    h2e_rotated = np.einsum(
        "abcd,Aa,Bb,Cc,Dd->ABCD",
        h2e,
        orbital_rotation,
        orbital_rotation.conj(),
        orbital_rotation,
        orbital_rotation.conj(),
        optimize="greedy",
    )
    return h1e_rotated, h2e_rotated

# Total number of spatial orbitals, including the bath sites and the impurity
# This should be an even number
norb = 20

# System is half-filled
nelec = (norb // 2, norb // 2)
# One orbital is the impurity, the rest are bath sites
n_bath = norb - 1

# Hamiltonian parameters
hybridization = 1.0
hopping = 1.0
onsite = 10.0
chemical_potential = -0.5 * onsite

# Generate Hamiltonian in position basis
h1e, h2e = siam_hamiltonian(
    norb=norb,
    hopping=hopping,
    onsite=onsite,
    hybridization=hybridization,
    chemical_potential=chemical_potential,
)

# Rotate to momentum basis
orbital_rotation = momentum_basis(norb)
h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())
# In the momentum basis, the impurity is placed in the center
impurity_index = n_bath // 2

Als Nächstes erzeugen wir die Schaltkreise zur Erzeugung der Krylov-Basiszustände. Für jede Spin-Spezies ist der Anfangszustand $\ket{\psi_0}$ durch die Superposition aller möglichen Anregungen der drei Elektronen, die dem Fermi-Niveau am nächsten sind, in die 4 nächsten leeren Moden gegeben, ausgehend vom Zustand $|00\cdots 0011 \cdots 11\rangle$, und wird durch die Anwendung von sieben XXPlusYYGates realisiert. Die zeitentwickelten Zustände werden durch sukzessive Anwendungen eines Trotter-Schritts zweiter Ordnung erzeugt.

Für eine detailliertere Beschreibung dieses Modells und wie die Schaltkreise entworfen wurden, siehe "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization".

from typing import Sequence

import ffsim
import scipy
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit import CircuitInstruction, Qubit
from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate

def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):
    """Prepare initial state."""
    x_gate = XGate()
    rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)
    for i in range(nocc):
        yield CircuitInstruction(x_gate, [qubits[i]])
        yield CircuitInstruction(x_gate, [qubits[norb + i]])
    for i in range(3):
        for j in range(nocc - i - 1, nocc + i, 2):
            yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])
            yield CircuitInstruction(
                rot, [qubits[norb + j], qubits[norb + j + 1]]
            )
    yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])
    yield CircuitInstruction(
        rot, [qubits[norb + j + 1], qubits[norb + j + 2]]
    )

def trotter_step(
    qubits: Sequence[Qubit],
    time_step: float,
    one_body_evolution: np.ndarray,
    h2e: np.ndarray,
    impurity_index: int,
    norb: int,
):
    """A Trotter step."""
    # Assume the two-body interaction is just the on-site interaction of the impurity
    onsite = h2e[
        impurity_index, impurity_index, impurity_index, impurity_index
    ]
    # Two-body evolution for half the time
    yield CircuitInstruction(
        CPhaseGate(-0.5 * time_step * onsite),
        [qubits[impurity_index], qubits[norb + impurity_index]],
    )
    # One-body evolution for the full time
    yield CircuitInstruction(
        ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits
    )
    # Two-body evolution for half the time
    yield CircuitInstruction(
        CPhaseGate(-0.5 * time_step * onsite),
        [qubits[impurity_index], qubits[norb + impurity_index]],
    )

# Time step
time_step = 0.2
# Number of Krylov basis states
krylov_dim = 8

# Initialize circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# Generate initial state
for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):
    circuit.append(instruction)
circuit.measure_all()

# Create list of circuits, starting with the initial state circuit
circuits = [circuit.copy()]

# Add time evolution circuits to the list
one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)
for i in range(krylov_dim - 1):
    # Remove measurements
    circuit.remove_final_measurements()
    # Append another Trotter step
    for instruction in trotter_step(
        qubits,
        time_step,
        one_body_evolution,
        h2e_momentum,
        impurity_index,
        norb,
    ):
        circuit.append(instruction)
    # Measure qubits
    circuit.measure_all()
    # Add a copy of the circuit to the list
    circuits.append(circuit.copy())

circuits[0].draw("mpl", scale=0.4, fold=-1)

circuits[-1].draw("mpl", scale=0.4, fold=-1)

Schritt 2: Problem für Quantenausführung optimieren

Nachdem wir die Schaltkreise erstellt haben, können wir sie für eine Ziel-Hardware optimieren. Wir wählen die am wenigsten ausgelastete QPU mit mindestens 127 Qubits. Weitere Informationen findest du in der Qiskit IBM® Runtime-Dokumentation.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=127
)
print(f"Using backend {backend.name}")

Using backend ibm_fez

Nun verwenden wir Qiskit, um die Schaltkreise für das Ziel-Backend zu transpilieren.

from qiskit.transpiler import generate_preset_pass_manager

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
isa_circuits = pass_manager.run(circuits)

Schritt 3: Ausführung mit Qiskit-Primitiven

Nachdem wir die Schaltkreise für die Hardware-Ausführung optimiert haben, sind wir bereit, sie auf der Ziel-Hardware auszuführen und Stichproben für die Grundzustandsenergieschätzung zu sammeln. Nachdem wir die Sampler-Primitive verwendet haben, um Bitfolgen aus jedem Schaltkreis zu ziehen, kombinieren wir alle Ergebnisse in einem einzigen Zähler-Wörterbuch und zeichnen die 20 am häufigsten gezogenen Bitfolgen auf.

from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import SamplerV2 as Sampler

# Sample from the circuits
sampler = Sampler(backend)
job = sampler.run(isa_circuits, shots=500)

from qiskit.primitives import BitArray

# Combine the counts from the individual Trotter circuits
bit_array = BitArray.concatenate_shots(
    [result.data.meas for result in job.result()]
)

plot_histogram(bit_array.get_counts(), number_to_keep=20)

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Nun führen wir den SQD-Algorithmus mit der Funktion diagonalize_fermionic_hamiltonian aus. Erläuterungen zu den Argumenten dieser Funktion findest du in der API-Dokumentation.

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
)

# List to capture intermediate results
result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

rng = np.random.default_rng(24)
result = diagonalize_fermionic_hamiltonian(
    h1e_momentum,
    h2e_momentum,
    bit_array,
    samples_per_batch=100,
    norb=norb,
    nelec=nelec,
    num_batches=3,
    max_iterations=5,
    symmetrize_spin=True,
    callback=callback,
    seed=rng,
)

Iteration 1
	Subsample 0
		Energy: -28.61321893815165
		Subspace dimension: 10609
	Subsample 1
		Energy: -28.628985564542244
		Subspace dimension: 13924
	Subsample 2
		Energy: -28.620151775558114
		Subspace dimension: 10404
Iteration 2
	Subsample 0
		Energy: -28.656893066053115
		Subspace dimension: 34225
	Subsample 1
		Energy: -28.65277622004119
		Subspace dimension: 38416
	Subsample 2
		Energy: -28.670856034959165
		Subspace dimension: 39601
Iteration 3
	Subsample 0
		Energy: -28.684787675404362
		Subspace dimension: 42436
	Subsample 1
		Energy: -28.676984757118426
		Subspace dimension: 50176
	Subsample 2
		Energy: -28.671581704249885
		Subspace dimension: 40804
Iteration 4
	Subsample 0
		Energy: -28.6859683054753
		Subspace dimension: 47961
	Subsample 1
		Energy: -28.69418206537316
		Subspace dimension: 51529
	Subsample 2
		Energy: -28.686083516445752
		Subspace dimension: 51529
Iteration 5
	Subsample 0
		Energy: -28.694665630711178
		Subspace dimension: 50625
	Subsample 1
		Energy: -28.69505984237118
		Subspace dimension: 47524
	Subsample 2
		Energy: -28.6942873883992
		Subspace dimension: 48841

Die folgende Code-Zelle zeichnet die Ergebnisse auf. Die erste Grafik zeigt die berechnete Energie als Funktion der Anzahl der Konfigurationswiederherstellungsiterationen, und die zweite Grafik zeigt die durchschnittliche Besetzung jedes räumlichen Orbitals nach der letzten Iteration. Für die Referenzenergie verwenden wir die Ergebnisse einer DMRG-Berechnung, die separat durchgeführt wurde.

import matplotlib.pyplot as plt

dmrg_energy = -28.70659686

min_es = [
    min(result, key=lambda res: res.energy).energy
    for result in result_history
]
min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])

# Data for energies plot
x1 = range(len(result_history))

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, min_es, label="energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].axhline(
    y=dmrg_energy, color="#BF5700", linestyle="--", label="DMRG energy"
)
axs[0].set_title("Approximated Ground State Energy vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Reference (DMRG) energy: {dmrg_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - dmrg_energy):.5f}")
plt.tight_layout()
plt.show()

Reference (DMRG) energy: -28.70660
SQD energy: -28.69506
Absolute error: 0.01154

Verifizierung der Energie

Die von SQD zurückgegebene Energie ist garantiert eine obere Schranke für die wahre Grundzustandsenergie. Der Wert der Energie kann verifiziert werden, da SQD auch die Koeffizienten des Zustandsvektors zurückgibt, der den Grundzustand approximiert. Du kannst die Energie aus dem Zustandsvektor unter Verwendung seiner 1- und 2-Teilchen-reduzierten Dichtematrizen berechnen, wie in der folgenden Code-Zelle demonstriert.

rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)
rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)

energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)

print(f"Recomputed energy: {energy:.5f}")

Recomputed energy: -28.69506

Referenzen

• Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization (SKQD paper)
• Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer (SQD paper)
• Diagonalization of large many-body Hamiltonians on a quantum processor (KQD paper)