Sample-based Krylov quantum diagonalization of a fermionic lattice model

Bruukschätzung: Neggen Sekunnen op'n Heron r2 Prozessor (HENWIES: Dat is blots en Schätzung. Dien Looptiet kann anners ween.)

Achtergrund

Dit Tutorial wiest, wo man de stichpravenbaseerde Quantendiagonaliseerung (SQD) bruuken deit, üm de Grundtoestandsenergie vun en fermioonschen Gittermodell to schätten. Speziell studeert wi dat eendimensionale Einzelstöörstellenmodell vun Anderson (SIAM), dat bruukt warrt, üm magnetische Störstellen in Metallen to beschrieven.

Dit Tutorial folgt en ähnlichen Warklopen as dat verwandte Tutorial Stichpravenbaseerde Quantendiagonaliseerung vun en Chemie-Hamiltonian. Man en wichtigen Ünnerscheed liggt darin, wo de Quantenschaltkreisen upbaut sünd. Dat anner Tutorial bruukt en heuristischen Variatschonsansatz, de för Chemie-Hamiltonians mit potentiell Millionen vun Wesselwarkungsterms attraktiv is. Düt Tutorial bruukt liekers Schaltkreisen, de de Tietentwicklung dörch den Hamiltonian approximeert. Sölke Schaltkreisen künnt deep ween, wat düssen Vörgang beter för Anwennen op Gittermodellen maakt. De Toestandsvektoren, de vun düsse Schaltkreisen prepareerd warrt, billt de Basis för en Krylov-Ünnerruum, un as Resultat konvergeert de Algorithmus bewieslich un effizient to'n Grundtostand ünner passliche Annamen.

De Ansatz in düt Tutorial kann as en Kombinatschoon vun de Techniken in SQD un Krylov-Quantendiagonaliseerung (KQD) ankeken warrt. De kombineerte Ansatz warrt manchmaal as stichpravenbaseerde Krylov-Quantendiagonaliseerung (SQKD) bezeekt. Kiek Krylov-Quantendiagonaliseerung vun Gitter-Hamiltonians för en Tutorial över de KQD-Methood.

Düt Tutorial baseert op de Arbeit "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", na de för mehr Details wiest warrt kann.

Einzelstöörstellenmodell vun Anderson (SIAM)

De eendimensionale SIAM-Hamiltonian is en Summ ut dree Terms:

$$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 sünd $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ de fermioonschen Erzeugings-/Vernichtungsoperatoren för de $\mathbf{j}^{\textrm{te}}$ Bad-Stell mit Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ sünd Erzeugings-/Vernichtungsoperatoren för den Störstellenmodus, un $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ un $V$ sünd reelle Tallen, de de Hüpp-, Vör-Oort- un Hybridiseerungswesselwarkungen beschrievt, un $\varepsilon$ is en reelle Tall, de dat chemische Potential spezifizeert.

Beacht, dat de Hamiltonian en spezielle Instanz vun den algemenen Wesselwarkings-Elektronen-Hamiltonian is,

$$\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$ ut Eenkörperterms bestaht, de quadratisch in de fermioonschen Erzeugings- un Vernichtungsoperatoren sünd, un $H_2$ ut Tweekörperterms bestaht, de quartisch sünd. För dat SIAM gillt

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

un $H_1$ enthölt de resterlichen Terms in den Hamiltonian. Üm den Hamiltonian programmatisch darstostelln, spiekert wi de Matrix $h_{pq}$ un den Tensor $h_{pqrs}$.

Oorts- un Impulsbasen

Wegens de approximative Translatschonssymmetrie in $H_\textrm{bath}$ verwacht wi nich, dat de Grundtostand in de Oortsbasis (de Orbitalbasis, in de de Hamiltonian baven spezifizeert is) dünn besett is. De Leistung vun SQD is blots garanteert, wenn de Grundtostand dünn besett is, dat heit, wenn he bedüdend Gewicht blots op en lütte Antall vun Rekenbasistöständ hett. Üm de Dünnbesettheit vun den Grundtostand to verbeetern, föhrt wi de Simulatschoon in de Orbitalbasis dörch, in de $H_\textrm{bath}$ diagonal is. Wi nöömt düsse Basis de Impulsbasis. Wieldat $H_\textrm{bath}$ en quadratischen fermioonschen Hamiltonian is, kann he effizient dörch en Orbitalrotatschoon diagonaliseert warrn.

Approximative Tietentwicklung dörch den Hamiltonian

Üm de Tietentwicklung dörch den Hamiltonian to approximeern, bruukt wi en Trotter-Suzuki-Opdelung vun de twete 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}.$$

Ünner de Jordan-Wigner-Transformatschoon entspricht de Tietentwicklung dörch $H_2$ en enkel CPhase-Gate twischen de Spin-up- un Spin-down-Orbitalen an de Störstell. Wieldat $H_1$ en quadratischen fermioonschen Hamiltonian is, entspricht de Tietentwicklung dörch $H_1$ en Orbitalrotatschoon.

De Krylov-Basistöständ $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, wobei $D$ de Dimenschoon vun den Krylov-Ünnerruum is, warrt dörch wedderholt Anwennen vun en enkel Trotter-Schritt billt, so dat

$$|\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}.$$

In den folgenden SQD-baseerten Warklopen treckt wi Stichpraven ut düssen Satz vun Schaltkreisen un naaverarbeit den kombineerten Satz vun Bitfolgen mit SQD. Düsse Vörgang staht in Gegensatz to den in dat verwandte Tutorial Stichpravenbaseerde Quantendiagonaliseerung vun en Chemie-Hamiltonian bruukten, wo Stichpraven ut en enkel heuristischen Variatschonsschaltkreis trocken worrn sünd.

Vörutsettungen

Vör dat Anfangen vun düt Tutorial stell seker, dat du dat Folgende installeert hest:

• Qiskit SDK v1.0 oder höger, mit Ünner stütten för Visualiseerung
• Qiskit Runtime v0.22 oder höger (pip install qiskit-ibm-runtime)
• SQD Qiskit Addon v0.11 oder höger (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Schritt 1: Problem op en Quantenschaltkreis afbillen

Eerst toseers erzeught wi den SIAM-Hamiltonian in de Oortsbasis. De Hamiltonian warrt dörch de Matrix $h_{pq}$ un den Tensor $h_{pqrs}$ darstellt. Denn dreih wi em in de Impulsbasis. In de Oortsbasis platzeert wi de Störstell op de eerste Stell. Man wenn wi na de Impulsbasis dreih, verschuuvt wi de Störstell na en zentrale Stell, üm Wesselwarkungen mit annere Orbitalen to erlichtern.

# Added by doQumentation — installs packages not in the Binder environment
!pip install -q ffsim qiskit-addon-sqd

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

As Nächstes erzeught wi de Schaltkreisen för de Erzeugung vun de Krylov-Basistöständ. För jede Spin-Spezies is de Anfangstostand $\ket{\psi_0}$ dörch de Superpositsoon vun all möglichen Anregungen vun de dree Elektronen, de dat Fermi-Niveau an'n nächsten sünd, in de 4 nächsten leedigen Modi geven, utgohn vun den Tostand $|00\cdots 0011 \cdots 11\rangle$, un realiseert dörch dat Anwennen vun söven XXPlusYYGates. De tietentwickelten Töständ warrt dörch sukzessive Anwennen vun en Trotter-Schritt vun de twete Ordnung erzeught.

För en mehr detailleerte Beschrievung vun düt Modell un wo de Schaltkreisen entwaarpt worrn sünd, kiek "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 Quantenutföhrung optimeren

Nadem wi de Schaltkreisen opstellt hebbt, künnt wi se för en Ziel-Hardware optimeren. Wi wählt de am wenigsten utlastete QPU mit mindestens 127 Qubits. Kiek de Qiskit IBM® Runtime-Dokumentatschoon för mehr Informatschonen.

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

Nu bruukt wi Qiskit, üm de Schaltkreisen för dat Ziel-Backend to transpileren.

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: Utföhrung mit Qiskit-Primitiven

Nadem wi de Schaltkreisen för Hardware-Utföhrung optimeert hebbt, sünd wi parat, se op de Ziel-Hardware uttofören un Stichpraven för de Grundtostandsenergieschätzung to sammeln. Nadem wi de Sampler-Primitive bruukt hebbt, üm Bitfolgen ut jeden Schaltkreis to trecken, kombineert wi all Resultaten in en enkel Teller-Wöörbook un wiest de 20 an't häufigsten trocken Bitfolgen.

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: Naaverarbeitung un Returgaav vun dat Resultat in dat wünschte klassische Format

Nu föhrt wi den SQD-Algorithmus mit de Funkschoon diagonalize_fermionic_hamiltonian ut. Kiek de API-Dokumentatschoon för Verklaren vun de Argumente vun düsse Funkschoon.

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

De folgende Code-Zell wiest de Resultaten. De eerste Grafik wiest de bereken Energie as Funkschoon vun de Antall vun de Konfiguratschonswedderharstellungsiteratschonen, un de twete Grafik wiest de döörsnittliche Besettung vun jeden rüümlichen Orbital na de letzde Iteratschoon. För de Referenzenergie bruukt wi de Resultaten vun en DMRG-Bereknung, de separat döörföhrt worrn is.

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

Verifizeren vun de Energie

De vun SQD returgeven Energie is garanteert en bövere Grenz för de wohre Grundtostandsenergie. De Weert vun de Energie kann verifizeert warrn, wieldat SQD ook de Koeffizienten vun den Toestandsvektor returgifft, de den Grundtostand approximeert. Du kannst de Energie ut düssen Toestandsvektor mit siene 1- un 2-Partikel-reduzeerten Dichtematrizen bereken, so as in de folgende Code-Zell demonstreert warrt.

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)