Diagonalizacja kwantowa Krylova oparta na próbkowaniu dla fermionowego modelu sieciowego

Szacowane użycie: dziewięć sekund na procesorze Heron r2 (UWAGA: To jest tylko szacunek. Rzeczywisty czas wykonania może się różnić.)

Wprowadzenie

Ten samouczek pokazuje, jak użyć kwantowej diagonalizacji opartej na próbkowaniu (SQD) do oszacowania energii stanu podstawowego fermionowego modelu sieciowego. Konkretnie badamy jednowymiarowy jednodomieszkowymodel Andersona (SIAM), który służy do opisu magnetycznych domieszek osadzonych w metalach.

Ten samouczek podąża podobnym schematem pracy do powiązanego samouczka Kwantowa diagonalizacja oparta na próbkowaniu dla hamiltonianu chemicznego. Kluczowa różnica polega jednak na sposobie budowania obwodów kwantowych. Tamten samouczek używa heurystycznego wariacyjnego ansatzu, który sprawdza się dla hamiltonianów chemicznych z potencjalnie milionami wyrazów oddziaływań. Natomiast ten samouczek używa obwodów przybliżających ewolucję czasową przez hamiltonian. Takie obwody mogą być głębokie, co czyni to podejście lepszym do zastosowań w modelach sieciowych. Wektory stanu przygotowane przez te obwody tworzą bazę dla podprzestrzeni Krylova, dzięki czemu algorytm w sposób dowiedziony i efektywny zbiega do stanu podstawowego, przy odpowiednich założeniach.

Podejście zastosowane w tym samouczku można traktować jako połączenie technik używanych w SQD oraz kwantowej diagonalizacji Krylova (KQD). To połączone podejście jest czasem określane jako kwantowa diagonalizacja Krylova oparta na próbkowaniu (SQKD). Zapoznaj się z samouczkiem Kwantowa diagonalizacja Krylova hamiltonianów sieciowych, aby dowiedzieć się więcej o metodzie KQD.

Ten samouczek opiera się na pracy "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", do której można sięgnąć po więcej szczegółów.

Jednodomieszkowymodel Andersona (SIAM)

Jednowymiarowy hamiltonian SIAM jest sumą trzech wyrazów:

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

gdzie

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

Tutaj $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ są fermionowymi operatorami kreacji/anihilacji dla $\mathbf{j}^{\textrm{th}}$ węzła kąpielowego ze spinem $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ są operatorami kreacji/anihilacji dla trybu domieszki, a $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ i $V$ są liczbami rzeczywistymi opisującymi oddziaływania przeskokowe, jednocentrowe i hybrydyzacyjne, a $\varepsilon$ jest liczbą rzeczywistą określającą potencjał chemiczny.

Zauważ, że hamiltonian jest szczególnym przypadkiem ogólnego hamiltonianu elektronów oddziałujących,

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

gdzie $H_1$ składa się z wyrazów jednocząsteczkowych, kwadratowych względem operatorów kreacji i anihilacji fermionów, a $H_2$ składa się z wyrazów dwucząsteczkowych, czwartego stopnia. Dla SIAM,

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

a $H_1$ zawiera pozostałe wyrazy hamiltonianu. Aby reprezentować hamiltonian programistycznie, przechowujemy macierz $h_{pq}$ oraz tensor $h_{pqrs}$.

Bazy położeniowa i pędowa

Ze względu na przybliżoną symetrię translacyjną w $H_\textrm{bath}$, nie oczekujemy, że stan podstawowy będzie rzadki w bazie położeniowej (bazie orbitalnej, w której hamiltonian jest określony powyżej). Wydajność SQD jest gwarantowana tylko wtedy, gdy stan podstawowy jest rzadki, tzn. ma znaczący udział tylko w niewielkiej liczbie stanów bazy obliczeniowej. Aby poprawić rzadkość stanu podstawowego, przeprowadzamy symulację w bazie orbitalnej, w której $H_\textrm{bath}$ jest diagonalny. Tę bazę nazywamy bazą pędową. Ponieważ $H_\textrm{bath}$ jest kwadratowym hamiltonianem fermionowym, można go efektywnie zdiagonalizować przez obrót orbitalny.

Przybliżona ewolucja czasowa przez hamiltonian

Do przybliżenia ewolucji czasowej przez hamiltonian używamy rozkładu Trottera-Suzukiego drugiego rzędu,

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

W ramach transformacji Jordana-Wignera, ewolucja czasowa przez $H_2$ sprowadza się do pojedynczej bramki CPhase między orbitalami spinowymi w górę i w dół na miejscu domieszki. Ponieważ $H_1$ jest kwadratowym hamiltonianem fermionowym, ewolucja czasowa przez $H_1$ sprowadza się do obrotu orbitalnego.

Stany bazowe Krylova $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, gdzie $D$ jest wymiarem podprzestrzeni Krylova, są tworzone przez wielokrotne zastosowanie pojedynczego kroku Trottera, więc

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

W następującym schemacie pracy opartym na SQD będziemy próbkować z tego zestawu obwodów i przetwarzać końcowy zbiór łańcuchów bitów za pomocą SQD. To podejście kontrastuje z podejściem zastosowanym w powiązanym samouczku Kwantowa diagonalizacja oparta na próbkowaniu dla hamiltonianu chemicznego, gdzie próbki były pobierane z jednego heurystycznego obwodu wariacyjnego.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane następujące pakiety:

• Qiskit SDK v1.0 lub nowszy, z obsługą wizualizacji
• Qiskit Runtime v0.22 lub nowszy (pip install qiskit-ibm-runtime)
• SQD Qiskit addon v0.11 lub nowszy (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Krok 1: Odwzorowanie problemu na obwód kwantowy

Najpierw generujemy hamiltonian SIAM w bazie położeniowej. Hamiltonian jest reprezentowany przez macierz $h_{pq}$ oraz tensor $h_{pqrs}$. Następnie obracamy go do bazy pędowej. W bazie położeniowej umieszczamy domieszkę na pierwszym miejscu. Jednak gdy przechodzimy do bazy pędowej, przenosimy domieszkę na centralne miejsce, aby ułatwić oddziaływania z innymi orbitalami.

# 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

Następnie generujemy obwody do wytworzenia stanów bazowych Krylova. Dla każdego gatunku spinowego stan początkowy $\ket{\psi_0}$ jest dany przez superpozycję wszystkich możliwych wzbudzeń trzech elektronów najbliższych poziomowi Fermiego do 4 najbliższych pustych modów zaczynając od stanu $|00\cdots 0011 \cdots 11\rangle$, i realizowany przez zastosowanie siedmiu bramek XXPlusYYGate. Stany ewoluowane w czasie są produkowane przez kolejne zastosowania kroku Trottera drugiego rzędu.

Aby uzyskać bardziej szczegółowy opis tego modelu i sposobu projektowania obwodów, zapoznaj się z "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)

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Teraz, gdy mamy już utworzone obwody, możemy je zoptymalizować pod kątem docelowego sprzętu. Wybieramy najmniej obciążony QPU z co najmniej 127 qubitami. Więcej informacji znajdziesz w dokumentacji Qiskit IBM® Runtime.

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

Teraz używamy Qiskit, aby skompilować obwody (transpile) pod kątem docelowego backendu.

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)

Krok 3: Wykonanie przy użyciu prymitywów Qiskit

Po optymalizacji obwodów pod kątem wykonania na sprzęcie jesteśmy gotowi uruchomić je na docelowym sprzęcie i zebrać próbki do estymacji energii stanu podstawowego. Po użyciu prymitywu Sampler do próbkowania ciągów bitów z każdego obwodu łączymy wszystkie wyniki w jeden słownik zliczeń i rysujemy 20 najczęściej próbkowanych ciągów bitów.

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)

Krok 4: Przetwarzanie końcowe i zwrócenie wyniku w pożądanym formacie klasycznym

Teraz uruchamiamy algorytm SQD, używając funkcji diagonalize_fermionic_hamiltonian. Wyjaśnienia argumentów tej funkcji znajdziesz w dokumentacji API.

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

Poniższa komórka kodu rysuje wyniki. Pierwszy wykres przedstawia obliczoną energię w funkcji liczby iteracji odtwarzania konfiguracji, a drugi wykres pokazuje średnie obsadzenie każdego orbitalu przestrzennego po ostatniej iteracji. Jako energię referencyjną używamy wyników obliczeń DMRG przeprowadzonych osobno.

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

Weryfikacja energii

Energia zwrócona przez SQD jest z gwarancją górnym ograniczeniem prawdziwej energii stanu podstawowego. Wartość energii można zweryfikować, ponieważ SQD zwraca również współczynniki wektora stanu przybliżającego stan podstawowy. Możesz obliczyć energię z wektora stanu, korzystając z jego 1- i 2-cząsteczkowych zredukowanych macierzy gęstości, co zostało zademonstrowane w poniższej komórce kodu.

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

Odnośniki

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