Sample-based Krylov quantum diagonalization of a fermionic lattice model

Nohzungsschätzung: Nüng Sekunde op enem Heron r2 Prozessor (HENWIES: Dat es nur en Schätzung. Ding Laufzick kann andisch sin.)

Hengerjrond

Hee dat Tutorial zeich, wi mer de stichprobebaseet Quantendiagonalisierung (SQD) nömme kann, öm de Jrondzostandsenergie vun enem fermionische Jittermodell ze schätze. Jesaht, studiere mer dat eendimensionale Einzelstörstelle-Anderson-Modell (SIAM), wat mer bruche för magnetische Störstelle en Metall ze beschrieve.

Dat Tutorial follich enem ähnliche Afloop wi dat verwandte Tutorial Stichprobebaseet Quantendiagonalisierung vun enem Chemie-Hamiltonian. Ävver ne wichtije Ongerscheed litt dren, wi de Quanteschaltkreise opjebaut sin. Dat ander Tutorial bruch ene heuristische Variationsansatz, dä för Chemie-Hamiltonians met potenziell Millione vun Wechselwirkungsterme aantrekklich es. Hee dat Tutorial hinjäje bruch Schaltkreise, di de Zickentwicklung dörch dä Hamiltonian approximeere. Sulch Schaltkreise künne deef sin, wat hee dä Vorjonn besser för Aanwendunge op Jittermodelle määt. De Zostandsvektore, di vun hee dä Schaltkreise präparreet wäde, bilde de Basis för ene Krylov-Ongerdraum, un als Resultat konverjeert dä Algorithmus bewieslech un effizient zom Jrondzostand onger jeejnete Aannahme.

Dä Ansatz en hee däm Tutorial kann als en Kombinazion vun de Technike vun SQD un Krylov-Quantendiagonalisierung (KQD) jesinn wäde. Dä kombineet Ansatz weed manchmol als stichprobebaseet Krylov-Quantendiagonalisierung (SQKD) bezeechnet. Luhr ens Krylov-Quantendiagonalisierung vun Jitter-Hamiltonians för e Tutorial övver de KQD-Metod.

Dat Tutorial baseet op dä Arbeid "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", op di mer för mieh Details verwies.

Einzelstörstelle-Anderson-Modell (SIAM)

Dä eendimensionale SIAM-Hamiltonian es en Summ us drei Terme:

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

wo

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

Hee, $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ sin de fermionische Erzeugongs-/Vernichtungsoperatore för de $\mathbf{j}^{\textrm{te}}$ Bad-Stell met Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ sin Erzeugongs-/Vernichtungsoperatore för dä Störstellemodus, un $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ un $V$ sin real Zahle, di de Hüpf-, Vor-Ort- un Hybridisierungswechselwirkunge beschrieve, un $\varepsilon$ es en real Zahl, di dat chemisch Potenziaal spezifizeert.

Maach Der bewoss, dat dä Hamiltonian en spezifisch Instanz vum allgemeine Wechselwirkungs-Elektrone-Hamiltonian es,

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

wo $H_1$ us Einkörpertermen bestonn deit, di quadratisch en de fermionische Erzeugongs- un Vernichtungsoperatore sin, un $H_2$ us Zweikörpertermen bestonn deit, di quartisch sin. För dat SIAM jilt

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

un $H_1$ enthält de reschteliche Terme em Hamiltonian. Öm dä Hamiltonian programmatesch darzestelle, speichere mer de Matrix $h_{pq}$ un dä Tensor $h_{pqrs}$.

Orts- un Impulsbasen

Wejen dä ungefähre Translazionssymmetrie en $H_\textrm{bath}$ erwade mer nit, dat dä Jrondzostand en dä Ortsbasis (dä Orbitalbasis, en dä dä Hamiltonian bovve spezifizeert es) dünn besatz es. De Leistung vun SQD es nur jarranteet, wann dä Jrondzostand dünn besatz es, dat heiß, wann hä sinnvoll Jewicht nur op en kleine Zahl vun Rechenbasis-Zöständ hät. Öm de Dünnbesatztheit vum Jrondzostand ze verbessere, maache mer de Simulazion en dä Orbitalbasis, en dä $H_\textrm{bath}$ diagonal es. Mer nenne hee die Basis de Impulsbasis. Weil $H_\textrm{bath}$ ne quadratische fermionische Hamiltonian es, kann hä effizient dörch en Orbitalrotazion diagonaliseert wäde.

Ungefähre Zickentwicklung dörch dä Hamiltonian

Öm de Zickentwicklung dörch dä Hamiltonian ze approximeere, bruche mer en Trotter-Suzuki-Zerlejung 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}.$$

Onger dä Jordan-Wigner-Transformazion entsprich de Zickentwicklung dörch $H_2$ enem einzelne CPhase-Gate zwesche de Spin-up- un Spin-down-Orbitale an dä Störstell. Weil $H_1$ ne quadratische fermionische Hamiltonian es, entsprich de Zickentwicklung dörch $H_1$ ener Orbitalrotazion.

De Krylov-Basiszöständ $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, wo $D$ de Dimension vum Krylov-Ongerdraum es, wäde dörch widderholt Aanwendung vun enem einzelne Trotter-Schrett jebildt, esu 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}.$$

Em folgende SQD-baseet Afloop wäde mer Stichprobe us hee däm Satz vun Schaltkreise träcke un dä kombineet Satz vun Bitfolge met SQD nohbearbeide. Hee dä Vorjonn stonn em Jejensatz zo däm em verwandte Tutorial Stichprobebaseet Quantendiagonalisierung vun enem Chemie-Hamiltonian jebruchte, wo Stichprobe us enem einzelne heuristische Variazionsschaltkreis jetrocke woode sin.

Vorussetzunge

Vörm Aanfang vun hee däm Tutorial sorg doför, dat De Foljendes installeert häs:

• Qiskit SDK v1.0 oder höher, met Ongerstö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)

Schrett 1: Problem op ene Quanteschaltkreis afbilde

Zoeesch erzeuche mer dä SIAM-Hamiltonian en dä Ortsbasis. Dä Hamiltonian weed dörch de Matrix $h_{pq}$ un dä Tensor $h_{pqrs}$ daarejestellt. Dann drieie mer en en de Impulsbasis. En dä Ortsbasis platzere mer de Störstell op dä eezte Stell. Ävver wann mer en de Impulsbasis drieie, verscheeve mer de Störstell en en zentrale Stell, öm Wechselwirkunge met andere Orbitale ze erleichtere.

# 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

Als Nöchstes erzeuche mer de Schaltkreise för de Erzeugung vun de Krylov-Basiszöständ. För jede Spin-Spezies es dä Aanfangszustand $\ket{\psi_0}$ dörch de Superposizion vun all mööjliche Aanrejunge vun de drei Elektrone, di däm Fermi-Niveau am nöchste sin, en de 4 nöchste ledije Mode jejovve, usjonn vum Zostand $|00\cdots 0011 \cdots 11\rangle$, un realiseet dörch de Aanwendung vun sibbe XXPlusYYGates. De zickentwickelte Zöständ wäde dörch sukzessive Aanwendunge vun enem Trotter-Schrett zweiter Ordnung erzeucht.

För en mieh detailleet Beschrievung vun hee däm Modell un wi de Schaltkreise entworf woode sin, luhr "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)

Schrett 2: Problem för Quantenusföhrung optimiere

Nohm mer de Schaltkreise erschtallt han, künne mer se för en Ziel-Hardware optimiere. Mer wähle de am winnichste usjelasdete QPU met mingeschdens 127 Qubits. Luhr de Qiskit IBM® Runtime-Dokumentazion för mieh Informazione.

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

Jetz bruche mer Qiskit, öm de Schaltkreise för dat Ziel-Backend ze transpileere.

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)

Schrett 3: Usföhrung met Qiskit-Primitive

Nohm mer de Schaltkreise för Hardware-Usföhrung optimeet han, sin mer paraat, se op dä Ziel-Hardware uszefööre un Stichprobe för de Jrondzostandsenergischätzung ze sammele. Nohm mer de Sampler-Primitive jebruch han, öm Bitfolge us jedem Schaltkreis ze träcke, kombinere mer all Resultat en enem einzelne Zähler-Wööterbooch un zeechne de 20 am häufichste jetrocke Bitfolge op.

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)

Schrett 4: Nohbearbeitung un Returgov vum Resultat em jewönschte klassische Format

Jetz fööre mer dä SQD-Algorithmus met dä Funkzion diagonalize_fermionic_hamiltonian us. Luur de API-Dokumentazion för Erklärunge vun de Argumente vun hee dä Funkzion.

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 foljende Code-Zell zeich de Resultat op. De eezte Jrafik zeich de berächnete Energie als Funkzion vun dä Zahl vun de Konfigurazioneswidderhärstellungsiterazione, un de zweite Jrafik zeich de durchschnettliche Besatzung vun jedem räumliche Orbital noh dä letzde Iterazion. För de Referenzenergije bruche mer de Resultat vun ener DMRG-Berächnung, di separat durchjeföhrt woode es.

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

Verifizeerung vun dä Energije

De vun SQD returgejovvene Energije es jarranteet en bovvere Jrenz för de woohre Jrondzostandsenergie. Dä Wäät vun dä Energije kann verifizeert wäde, weil SQD och de Koeffiziente vum Zostandsvektor returgitt, dä dä Jrondzostand approximeet. Do kanns de Energije us däm Zostandsvektor met singene 1- un 2-Partikel-reduzeet Dichtematrizze berächne, wi en dä foljende Code-Zell demonstreet weed.

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

Referenze

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