Sample-based Krylov-Quantendiagonalisierung eines fermionischen Gittermodells

Schätzung des Ressourcenverbrauchs: Neun Sekunden auf einem Heron-r2-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine tatsächliche Laufzeit kann abweichen.)

Lernziele

Nach diesem Tutorial sollten Nutzer folgendes verstehen:

• Wie man das SQD Qiskit Addon verwendet, um die Grundzustandsenergie eines Gittermodells mithilfe von Bitstrings zu approximieren, die von einer Quantenverarbeitungseinheit (QPU) gesampelt wurden.
• Wie man ffsim verwendet, um Zeitentwicklungsschaltkreise für die fermionische Simulation zu konstruieren.
• Wie man Samples aus mehreren Circuits für die Nachverarbeitung mit dem sample-basierten Krylov-Diagonalisierungsalgorithmus (SKQD) kombiniert.

Voraussetzungen

Wir empfehlen, dass du mit folgenden Themen vertraut bist, bevor du dieses Tutorial durcharbeitest:

• Sample-based Quantendiagonalisierung eines Chemie-Hamiltonoperators
• Krylov-Quantendiagonalisierung von Gitter-Hamiltonoperatoren
• Qiskit Primitives

Hintergrund

Dieses Tutorial zeigt, wie man die sample-basierte Quantendiagonalisierung (SQD) verwendet, um die Grundzustandsenergie eines fermionischen Gittermodells zu schätzen. Konkret untersuchen wir das eindimensionale Single-Impurity-Anderson-Modell (SIAM), das zur Beschreibung magnetischer Verunreinigungen in Metallen verwendet wird.

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

Der in diesem Tutorial verwendete Ansatz kann als Kombination der Techniken aus SQD und Krylov-Quantendiagonalisierung (KQD) betrachtet werden. Der kombinierte Ansatz wird manchmal als sample-basierte Krylov-Quantendiagonalisierung (SQKD) bezeichnet. Unter Krylov-Quantendiagonalisierung von Gitter-Hamiltonoperatoren findest du 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.

Single-Impurity-Anderson-Modell (SIAM)

Der eindimensionale SIAM-Hamiltonoperator 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 den $\mathbf{j}^{\textrm{th}}$ Badplatz mit Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ sind Erzeugungs-/Vernichtungsoperatoren für den Verunreinigungsmodus, und $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ und $V$ sind reelle Zahlen, die die Hüpf-, On-Site- und Hybridisierungswechselwirkungen beschreiben, und $\varepsilon$ ist eine reelle Zahl, die das chemische Potential angibt.

Beachte, dass der Hamiltonoperator eine spezifische Instanz des generischen wechselwirkenden Elektronen-Hamiltonoperators 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öperpotentialen besteht, die quadratisch in den fermionischen Erzeugungs- und Vernichtungsoperatoren sind, und $H_2$ aus Zweiköperpotentialen, die quartisch sind. Für das SIAM gilt:

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

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

Orts- und Impulsbasis

Aufgrund der näherungsweisen Translationssymmetrie in $H_\textrm{bath}$ erwarten wir nicht, dass der Grundzustand in der Ortsbasis dünn besetzt ist (die Orbitalbasis, in der der Hamiltonoperator oben angegeben ist). Die Leistung von SQD ist nur dann garantiert, wenn der Grundzustand dünn besetzt ist, das heißt, er hat signifikantes Gewicht auf nur einer kleinen Anzahl von Rechenbasiszuständen. Um die Dünnbesetzung 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 Hamiltonoperator ist, kann er effizient durch eine Orbitalrotation diagonalisiert werden.

Näherungsweise Zeitentwicklung durch den Hamiltonoperator

Um die Zeitentwicklung durch den Hamiltonoperator 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 am Verunreinigungsplatz. Da $H_1$ ein quadratischer fermionischer Hamiltonoperator 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 aus dieser Menge von Schaltkreisen sampeln und die kombinierte Menge von Bitstrings mit SQD nachverarbeiten. Dieser Ansatz unterscheidet sich von dem im verwandten Tutorial Sample-based Quantendiagonalisierung eines Chemie-Hamiltonoperators, bei dem Samples aus einem einzelnen heuristischen Variationsschaltkreis gezogen wurden.

Anforderungen

Bevor du dieses Tutorial beginnst, stelle sicher, dass du Folgendes installiert hast:

• Qiskit SDK v1.0 oder höher, mit Visualisierungsunterstützung
• 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 v0.0.72 oder höher (pip install ffsim)

Kleinskaliges Simulator-Beispiel

Schritt 1: Problem auf einen Quantenschaltkreis abbilden

Zunächst generieren wir den SIAM-Hamiltonoperator in der Ortsbasis. Der Hamiltonoperator 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 Verunreinigung am ersten Platz. Wenn wir jedoch in die Impulsbasis wechseln, verschieben wir die Verunreinigung an einen zentralen Platz, um Wechselwirkungen mit anderen Orbitalen zu erleichtern.

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

import numpy as np
import pyscf.fci

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

# 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

# Use PySCF to compute the exact ground state energy
reference_energy, _ = pyscf.fci.direct_spin1.kernel(h1e, h2e, norb, nelec)

Als nächstes generieren wir die Schaltkreise zur Erzeugung der Krylov-Basiszustände. Für jede Spinspezies ist der Anfangszustand $\ket{\psi_0}$ als Superposition aller möglichen Anregungen der drei Elektronen, die der Fermi-Energie am nächsten sind, in die 4 nächstliegenden 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 aufeinanderfolgende Anwendungen eines Trotter-Schritts zweiter Ordnung erzeugt.

Für eine detailliertere Beschreibung dieses Modells und wie die Schaltkreise entworfen werden, 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."""
    assert norb >= 8
    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 die Quantenausführung optimieren

Als nächstes optimieren wir den Schaltkreis für eine Zielhardware. Zunächst erstellen wir ein generisches Backend mit einer bestimmten Anzahl von Qubits und einem Gate-Set, in das die Zeitentwicklungsschaltkreise auf natürliche Weise zerlegt werden.

from qiskit.providers.fake_provider import GenericBackendV2

backend = GenericBackendV2(
    2 * norb, basis_gates=["cp", "xx_plus_yy", "p", "x"]
)

Nun verwenden wir Qiskit, um die Schaltkreise zum 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: Mit Qiskit Primitives ausführen

Nach der Optimierung der Schaltkreise für die Hardwareausführung sind wir bereit, sie auf der Zielhardware auszuführen und Samples für die Grundzustandsenergieschätzung zu sammeln. Nachdem wir das Sampler Primitive verwendet haben, um Bitstrings aus jedem Schaltkreis zu sampeln, kombinieren wir alle Ergebnisse in einem einzigen Counts-Dictionary und stellen die 20 am häufigsten gesampelten Bitstrings dar.

from qiskit.visualization import plot_histogram
from qiskit.primitives import StatevectorSampler

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

from qiskit.primitives import BitArray

# Combine the shots 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: Nachverarbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Jetzt führen wir den SQD-Algorithmus mit der Funktion diagonalize_fermionic_hamiltonian aus. Erklärungen 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: -13.4222953188441
		Subspace dimension: 529
	Subsample 1
		Energy: -13.42237556285828
		Subspace dimension: 784
	Subsample 2
		Energy: -13.422045397387413
		Subspace dimension: 529
Iteration 2
	Subsample 0
		Energy: -13.422379583305478
		Subspace dimension: 900
	Subsample 1
		Energy: -13.422376197704326
		Subspace dimension: 841
	Subsample 2
		Energy: -13.422421162849295
		Subspace dimension: 1089
Iteration 3
	Subsample 0
		Energy: -13.422421164670345
		Subspace dimension: 1156
	Subsample 1
		Energy: -13.422421492737689
		Subspace dimension: 1156
	Subsample 2
		Energy: -13.422421205869572
		Subspace dimension: 1156
Iteration 4
	Subsample 0
		Energy: -13.422421494558726
		Subspace dimension: 1225
	Subsample 1
		Energy: -13.422421492737689
		Subspace dimension: 1156
	Subsample 2
		Energy: -13.422421492737689
		Subspace dimension: 1156

Die folgende Code-Zelle stellt die Ergebnisse dar. Der erste Plot zeigt die berechnete Energie als Funktion der Anzahl der Konfigurationswiederherstellungsiterationen, und der zweite Plot zeigt die durchschnittliche Besetzung jedes räumlichen Orbitals nach der letzten Iteration. Da es sich um ein so kleines Problem handelt, bringt uns die erste Iteration bereits sehr nahe an die exakte Energie heran (beachte die Skala der y-Achse).

import matplotlib.pyplot as plt

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=reference_energy,
    color="#BF5700",
    linestyle="--",
    label="reference 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 energy: {reference_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - reference_energy):.5f}")
plt.tight_layout()
plt.show()

Reference energy: -13.42249
SQD energy: -13.42242
Absolute error: 0.00007

Energie verifizieren

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 mithilfe seiner Ein- und Zweiteilchen-Dichtematrizen berechnen, wie in der folgenden Code-Zelle gezeigt.

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

Großskaliges Hardware-Beispiel

Jetzt führen wir ein größeres Beispiel auf einer echten QPU aus. Als Referenzenergie verwenden wir die Ergebnisse einer DMRG-Berechnung, die separat durchgeführt wurde.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import QiskitRuntimeService

# Model parameters
norb = 20
nelec = (norb // 2, norb // 2)
n_bath = norb - 1
hybridization = 1.0
hopping = 1.0
onsite = 10.0
chemical_potential = -0.5 * onsite

# Generate Hamiltonian and orbital rotation
h1e, h2e = siam_hamiltonian(
    norb=norb,
    hopping=hopping,
    onsite=onsite,
    hybridization=hybridization,
    chemical_potential=chemical_potential,
)
orbital_rotation = momentum_basis(norb)
h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())
impurity_index = n_bath // 2

# Set reference energy to DMRG value computed separately
reference_energy = -28.70659686

# Algorithm parameters
time_step = 0.2
krylov_dim = 8

# Construct circuits
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)
for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):
    circuit.append(instruction)
circuit.measure_all()
circuits = [circuit.copy()]
one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)
for i in range(krylov_dim - 1):
    circuit.remove_final_measurements()
    for instruction in trotter_step(
        qubits,
        time_step,
        one_body_evolution,
        h2e_momentum,
        impurity_index,
        norb,
    ):
        circuit.append(instruction)
    circuit.measure_all()
    circuits.append(circuit.copy())

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

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

# Sample from the circuits
sampler = Sampler(backend)
sampler.options.environment.job_tags = ["TUT_SKQD"]
job = sampler.run(isa_circuits, shots=500)

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

# Run configuration recovery and diagonalization
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,
)

# Plot results
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])
x1 = range(len(result_history))
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].plot(x1, min_es, label="energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].axhline(
    y=reference_energy,
    color="#BF5700",
    linestyle="--",
    label="reference 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()
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 energy: {reference_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - reference_energy):.5f}")
plt.tight_layout()
plt.show()

Using backend ibm_boston
Iteration 1
	Subsample 0
		Energy: -28.63965951544449
		Subspace dimension: 9801
	Subsample 1
		Energy: -28.625588929202006
		Subspace dimension: 9409
	Subsample 2
		Energy: -28.647371834135498
		Subspace dimension: 8281
Iteration 2
	Subsample 0
		Energy: -28.67213260849567
		Subspace dimension: 29584
	Subsample 1
		Energy: -28.670340686158816
		Subspace dimension: 27225
	Subsample 2
		Energy: -28.669976379525988
		Subspace dimension: 31329
Iteration 3
	Subsample 0
		Energy: -28.68622875601382
		Subspace dimension: 36100
	Subsample 1
		Energy: -28.698569623143126
		Subspace dimension: 34225
	Subsample 2
		Energy: -28.694848533971882
		Subspace dimension: 33856
Iteration 4
	Subsample 0
		Energy: -28.69883392844593
		Subspace dimension: 42025
	Subsample 1
		Energy: -28.701289495200996
		Subspace dimension: 38025
	Subsample 2
		Energy: -28.699319594978245
		Subspace dimension: 45369
Iteration 5
	Subsample 0
		Energy: -28.701936886834154
		Subspace dimension: 51076
	Subsample 1
		Energy: -28.702468711812013
		Subspace dimension: 53824
	Subsample 2
		Energy: -28.702298147575938
		Subspace dimension: 52900
Reference energy: -28.70660
SQD energy: -28.70247
Absolute error: 0.00413

Nächste Schritte

Empfehlungen

Wenn du diese Arbeit interessant findest, könnte dich folgendes Material interessieren:

• Sample-based Quantendiagonalisierung eines Chemie-Hamiltonoperators - ein verwandtes Tutorial mit einem heuristischen Variationsansatz anstelle von Trotter-Schaltkreisen
• Krylov-Quantendiagonalisierung von Gitter-Hamiltonoperatoren - ein Tutorial zur KQD-Methode
• SQD Addon API-Dokumentation - Referenz für die Funktion diagonalize_fermionic_hamiltonian
• Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization - das Paper, auf dem dieses Tutorial basiert