Sample-based Krylov quantum diagonalization of a fermionic lattice model

Bruuchsschätzig: Nüün Sekunde uf emne Heron r2 Prozessor (HIIWYS: Das isch nume e Schätzig. Dini Laufzyt cha variiere.)

Hingergrund

Das Tutorial zeigt, wie mer d stichproobebasisrti Quantediagonalisirig (SQD) verwändet, zum d Grundzuestandsenergie vomene fermionische Gittermodell z schätze. Spezifisch studiere mir s eididimänsional Einzelstörschtellemodell vo Anderson (SIAM), wo mer bruucht, zum magnetischi Störschtelle i Metall z beschriibe.

Das Tutorial folgt emene ähnliche Arbetsablauf wie s verwandt Tutorial Stichproobebasisrti Quantediagonalisirig vonere Chemie-Hamiltonian. Aber en wichtige Ungerschiid lit dinn, wie d Quanteschaltkreis ufbaut wärde. S ander Tutorial verwändet en heuristischi Variationsaasatz, wo für Chemie-Hamiltonians mit potenziell Millione vo Wächselwirkigstärm attraktiv isch. Das Tutorial da gäge verwändet Schaltkreis, wo d Zytentwiglig dür dr Hamiltonian approximiere. Sötchi Schaltkreis chöi tüüf si, was dä Vorgang besser für Aawändige uf Gittermodäll macht. D Zuestandsvektore, wo vo dene Schaltkreis präpariert wärde, bilde d Basis fürne Krylov-Unggerruem, und als Resultat konvergiert dr Algorithmus bewiislech und effizient zum Grundzuestand under gäignete Aanahme.

Dr Aasatz i däm Tutorial cha als e Kombination vo de Tächnike i SQD und Krylov-Quantediagonalisirig (KQD) betrachtet wärde. Dr kombiniert Aasatz wird mänggisch als stichproobebasisrti Krylov-Quantediagonalisirig (SQKD) bezeichnit. Lueg Krylov-Quantediagonalisirig vo Gitter-Hamiltonians füres Tutorial zur KQD-Methode.

Das Tutorial basiert uf dr Arbet "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", uf wo mer für mee Details verwyst.

Einzelstörschtellemodell vo Anderson (SIAM)

Dr eididimänsional SIAM-Hamiltonian isch e Summe us drü Tärm:

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

wobii

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

Doo si $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ d fermionische Erzeugigs-/Vernichtigsoperatore für d $\mathbf{j}^{\textrm{ti}}$ Bad-Schtell mit Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ si Erzeugigs-/Vernichtigsoperatore füre Störschtällemodus, und $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ und $V$ si reelli Zahle, wo d Hüpf-, Vor-Ort- und Hybridisierigswächselwirkige beschriibe, und $\varepsilon$ isch e reelli Zahl, wo s chemisch Potenziaal spezifiziert.

Biacht, ass dr Hamiltonian e spezifischi Instanz vom allgemeine Wächselwirkigs-Elektrone-Hamiltonian isch,

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

wobii $H_1$ us Eichorpertärm bschtaht, wo quadratisch i de fermionische Erzeugigs- und Vernichtigsoperatore si, und $H_2$ us Zweichorpertärm bschtaht, wo quartisch si. Fürs SIAM gilt

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

und $H_1$ enthaut di reschtliche Tärm im Hamiltonian. Zum dr Hamiltonian programmatisch daarzuschteile, spychere mir d Matrix $h_{pq}$ und dr Tensor $h_{pqrs}$.

Orts- und Impulsbase

Wäge dr approximative Translationsymmetrii i $H_\textrm{bath}$ erwarte mir nöd, ass dr Grundzuestand i dr Ortsbasis (dr Orbitalbasis, i wo dr Hamiltonian obe spezifiziert isch) dünn bsetzt isch. D Leistig vo SQD isch nume garantiert, wänn dr Grundzuestand dünn bsetzt isch, das heisst, wänn er signifikants Gwicht nume uf nere chlyne Aazahl vo Rächnebasiszueständ het. Zum d Dünnbsetztig vom Grundzuestand z verbessere, füere mir d Simulazioon i dr Orbitalbasis düre, i wo $H_\textrm{bath}$ diagonal isch. Mir nänne die Basis d Impulsbasis. Wil $H_\textrm{bath}$ e quadratische fermionische Hamiltonian isch, cha är effizient dür e Orbitalrotazioon diagonalisiert wärde.

Approximativi Zytentwiglig dür dr Hamiltonian

Zum d Zytentwiglig dür dr Hamiltonian z approximiere, verwände mir e Trotter-Suzuki-Zerlegig vo zweiter Ornig,

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

Under dr Jordan-Wigner-Transformazioon entspricht d Zytentwiglig dür $H_2$ emene einzelne CPhase-Gate zwüsche de Spin-up- und Spin-down-Orbitale a dr Störschtell. Wil $H_1$ e quadratische fermionische Hamiltonian isch, entspricht d Zytentwiglig dür $H_1$ nere Orbitalrotazioon.

D Krylov-Basiszueständ $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, wobii $D$ d Dimänsioon vom Krylov-Unggerruem isch, wärde dür wiederholti Aawändig vomene einzelne Trotter-Schritt bildet, so dass

$$|\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 folgende SQD-basisierte Arbetsablauf ziehe mir Stichproobe us däm Satz vo Schaltkreis und verarbeite dr kombiniert Satz vo Bitfolge mit SQD noche. Dä Vorgang schtaht im Gägesatz zu däm im verwandte Tutorial Stichproobebasisrti Quantediagonalisirig vonere Chemie-Hamiltonian verwändete, wo Stichproobe us emene einzelne heuristische Variazioonsschtaltkreis zoge worde si.

Voorussetzig

Vor em Aafang vo däm Tutorial schteume sicher, ass du s Folgende installiert hesch:

• Qiskit SDK v1.0 oder höcher, mit Ungerschtützig für Visualisirig
• Qiskit Runtime v0.22 oder höcher (pip install qiskit-ibm-runtime)
• SQD Qiskit Addon v0.11 oder höcher (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Schritt 1: Problem ufne Quanteschtaltkreis abblde

Zueersch generiere mir dr SIAM-Hamiltonian i dr Ortsbasis. Dr Hamiltonian wird dür d Matrix $h_{pq}$ und dr Tensor $h_{pqrs}$ daargschtellt. Denn rotiere mir en i d Impulsbasis. I dr Ortsbasis platziere mir d Störschtell uf dr erschte Schtell. Aber wänn mir zur Impulsbasis rotiere, verschüebe mir d Störschtell uf e zentrali Schtell, zum Wächselwirkige mit andere Orbitale z erliichtere.

# 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ächschts generiere mir d Schaltkreis zur Erzeugig vo de Krylov-Basiszueständ. Für jedi Spin-Spezies isch dr Aafangszuestand $\ket{\psi_0}$ dür d Superposizioon vo allne möögliche Aaregige vo de drüü Elektrone, wo em Fermi-Niveau am nöchschte si, i d 4 nöchschte läre Modi gää, usgend vom Zuestand $|00\cdots 0011 \cdots 11\rangle$, und realisiert dür d Aawändig vo sibne XXPlusYYGates. D zytentwiglete Zueständ wärde dür sukzessivi Aawändige vomene Trotter-Schritt vo zweiter Ornig erzügt.

Füre mee detaillierti Beschryybig vo däm Modell und wie d Schaltkreis entwoorfe worde si, lueg "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 Quanteuusführig optimiere

Nochdem mir d Schaltkreis erschtellt hei, chöi mir si für e Ziel-Hardware optimiere. Mir wehle di am wenigschte uusglaschteti QPU mit mindischtens 127 Qubits. Lueg d Qiskit IBM® Runtime-Dokumentazioon für mee Informatioone.

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

Jetzt verwände mir Qiskit, zum d Schaltkreis fürs Ziel-Backend z transpiliere.

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: Uusführig mit Qiskit-Primitive

Nochdem mir d Schaltkreis für Hardware-Uusführig optimiert hei, si mir parat, si uf dr Ziel-Hardware uusz'füere und Stichproobe für d Grundzuestandsenergischätzig z sammle. Nochdem mir d Sampler-Primitivi verwändet hei, zum Bitfolge us jedem Schaltkreis z zie, kombiniere mir alli Resoltaat i emene einzelne Zähler-Wörterbuech und zeige di 20 am häufigschte zogene Bitfolge aa.

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: Nochverarbeitig und Rugggab vom Resoltaat im gwünscht klassische Format

Jetzt füere mir dr SQD-Algorithmus mit dr Funkzioon diagonalize_fermionic_hamiltonian uus. Lueg d API-Dokumentazioon für Erkläärige vo de Argumänt vo dere Funkzioon.

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

Di folgende Code-Zälle zeigt d Resoltaat aa. Di erschti Grafik zeigt di berächneti Energii als Funkzioon vo dr Aazahl vo de Konfiguratioonswiderherstelligsiterazioon, und di zweiti Grafik zeigt d durchschnittlichi Bsätzig vo jedem rüümliche Orbital noch dr letscht Iterazioon. Für d Referänzenergie verwände mir d Resoltaat vonere DMRG-Berächnig, wo separat düregfüert worde isch.

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

Verifizierig vo dr Energii

Di vo SQD zruggää Energii isch garantiert e oberi Gränze für di waari Grundzuestandsenergie. Dr Wärt vo dr Energii cha verifiziert wärde, wil SQD au d Koeffiziente vom Zuestandsvektor zrugggit, wo dr Grundzuestand approximiert. Du chasch d Energii us däm Zuestandsvektor mit sine 1- und 2-Partikel-reduzierte Dichtematritze berächne, wie i dr folgende Code-Zälle demonstriert wird.

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

Referänze

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