Krylovova kvantová diagonalizace založená na vzorkování (SKQD)

Tato lekce o Krylovově kvantové diagonalizaci založené na vzorkování (SKQD) kombinuje metody vysvětlené v předchozích lekcích. Skládá se z jednoho příkladu využívajícího framework Qiskit patterns:

• Krok 1: Namapování problému na kvantové Circuit a operátory
• Krok 2: Optimalizace pro cílový hardware
• Krok 3: Spuštění pomocí Qiskit Primitives
• Krok 4: Postprocessing

Důležitým krokem v metodě kvantové diagonalizace založené na vzorkování je generování kvalitních vektorů pro podprostor. V předchozí lekci jsme používali ansatz LUCJ k generování vektorů podprostoru pro chemický hamiltonián. V této lekci budeme používat kvantové Krylovovy stavy[1], jak bylo probráno v lekci 2. Nejprve si zopakujeme, jak vytvořit Krylovův prostor na kvantovém počítači pomocí operací časového vývoje. Poté z něj budeme vzorkovat. Promítneme hamiltonián systému na vzorkovaný podprostor a diagonalizujeme ho, abychom odhadli energii základního stavu. Algoritmus prokazatelně a efektivně konverguje k základnímu stavu za předpokladů popsaných v lekci 2.

0. Krylovův prostor

Připomeň si, že Krylovův prostor $\mathcal{K}^r$ řádu $r$ je prostor generovaný vektory získanými násobením vyšších mocnin matice $A$, až do $r-1$, referenčním vektorem $\vert v \rangle$.

$$\mathcal{K}^r = \left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}$$

Je-li matice $A$ hamiltonián $H$, nazývá se odpovídající prostor mocninným Krylovým prostorem $\mathcal{K}_P$. V případě, kdy $A$ je operátor časového vývoje generovaný hamiltoniánem $U=e^{-iH(dt)}$, označujeme tento prostor jako unitární Krylovův prostor $\mathcal{K}_U$. Mocninný Krylovův podprostor nelze na kvantovém počítači generovat přímo, protože $H$ není unitární operátor. Místo toho můžeme použít operátor časového vývoje $U = e^{-iH(dt)}$, u nějž lze ukázat, že poskytuje podobné záruky konvergence jako mocninný Krylovův prostor. Mocniny $U$ pak odpovídají různým časovým krokům $U^k = e^{-iH(k dt)}$, kde $k = 0, 1, 2, ..., (r-1)$.

$$\mathcal{K}_U^r = \left\{ \vert \psi \rangle, U \vert \psi \rangle, U^2 \vert \psi \rangle, ..., U^{r-1} \vert \psi \rangle \right\}$$

1. Namapování problému na kvantové Circuit a operátory

V této lekci uvažujeme hamiltonián antiferomagnetického XX-Z spinového řetězce se spinem 1/2 o $L = 22$ místech s periodickou okrajovou podmínkou:

$$H = \sum_{i, j}^{N} J_{xy} (X_{i} X_{j} + Y_{i} Y_{j}) + Z_{i} Z_{j}$$

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

from qiskit.transpiler import CouplingMap
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

num_spins = 22
coupling_map = CouplingMap.from_ring(num_spins)
H_op = generate_xyz_hamiltonian(coupling_map, coupling_constants=(0.3, 0.3, 1.0))

Pro sestavení Krylovova prostoru potřebujeme tři hlavní ingredience:

1. Volbu Krylovovy dimenze ($r$) a časového kroku ($dt$).
2. Počáteční (referenční) stav (vektor $\vert v \rangle$ výše) s polynomiálním překryvem s cílovým (základním) stavem, přičemž cílový stav je řídký. Tento požadavek na polynomiální překryv je stejný jako v algoritmu kvantové fázové estimace.
3. Operátory časového vývoje $U^{k}=e^{-iH(k * dt)}$ ($k = 0, 1, 2, ..., r-1$).

Pro zvolenou hodnotu $r$ (a $dt$) vytvoříme $r$ samostatných kvantových Circuit a budeme z nich vzorkovat. Každý kvantový Circuit se vytvoří spojením kvantového Circuit reprezentujícího referenční stav a operátoru časového vývoje pro danou hodnotu $k$.

Větší Krylovova dimenze zlepšuje konvergenci odhadované energie. V této lekci nastavíme dimenzi na $5$, abychom ilustrovali trend konvergence.

Ref [2] ukázal, že dostatečně malý časový krok pro KQD je $\pi / \vert \vert H \vert \vert$ a že je výhodnější tuto hodnotu podestimovat než přestimovat. Na druhou stranu příliš malé $dt$ vede k horšímu podmínění Krylovova podprostoru, protože se vektory Krylovovy báze od časového kroku k časovému kroku liší méně. Přestože tato volba $dt$ je prokazatelně dostatečná pro konvergenci SKQD, v tomto kontextu vzorkování je optimální volba $dt$ v praxi předmětem probíhajícího výzkumu. V této lekci nastavíme $dt = 0.15$.

Kromě Krylovovy dimenze a časového kroku musíme nastavit počet Trotterových kroků pro časový vývoj. Příliš málo kroků vede k větším chybám trotterizace, zatímco příliš mnoho kroků vede k hlubším Circuit. V této lekci nastavíme počet Trotterových kroků na $6$.

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of krylov subspace
dt = 0.15
num_trotter_steps = 6

Dále musíme zvolit referenční stav $\vert \psi \rangle$, který má určitý překryv se základním stavem. Pro tento hamiltonián používáme jako referenční stav Neelův stav se střídajícími se jedničkami a nulami $\vert ...101...010...101 \rangle$.

# Prep `Neel` state as the reference state for evolution
from qiskit import QuantumCircuit

qc_state_prep = QuantumCircuit(num_spins)
for i in range(num_spins):
    if i % 2 == 0:
        qc_state_prep.x(i)

Nakonec musíme namapovat operátor časového vývoje na kvantový Circuit. To bylo provedeno v lekci 2, ale zde využijeme metody v Qiskitu, konkrétně metodu nazvanou synthesis. Existují různé metody syntézy matematických operátorů na kvantové Circuit s kvantovými Gate. Mnoho takových technik je dostupných v modulu syntézy Qiskit. Použijeme přístup LieTrotter pro syntézu [3] [4].

from qiskit.circuit import QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

evol_gate = PauliEvolutionGate(
    H_op, time=(dt / num_trotter_steps), synthesis=LieTrotter(reps=num_trotter_steps)
)  # `U` operator

qr = QuantumRegister(num_spins)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

circuits = []
for rep in range(krylov_dim):
    circ = qc_state_prep.copy()

    # Repeating the `U` operator to implement U^0, U^1, U^2, and so on, for power Krylov space
    for _ in range(rep):
        circ.compose(other=qc_evol, inplace=True)

    circ.measure_all()
    circuits.append(circ)

circuits[1].decompose().draw("mpl", fold=-1)

circuits[2].decompose().draw("mpl", fold=-1)

2. Optimalizace pro cílový hardware

Nyní, když jsme vytvořili Circuit, je můžeme optimalizovat pro cílový hardware. Vybereme QPU ve škále utility.

import warnings

from qiskit import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService

warnings.filterwarnings("ignore")

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Nyní transpilujeme Circuit na cílový Backend pomocí předpřipraveného správce průchodů.

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_circuits = pm.run(circuits=circuits)

3. Spuštění na cílovém hardware

Po optimalizaci Circuit pro spuštění na hardware jsme připraveni je spustit na cílovém hardware a shromáždit vzorky pro odhad energie základního stavu.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
job = sampler.run(isa_circuits, shots=100_000)  # Takes approximately 2m 58s of QPU time

counts_all = [job.result()[k].data.meas.get_counts() for k in range(krylov_dim)]

4. Postprocessing výsledků

Dále agregujeme počty pro rostoucí Krylovovy dimenze kumulativním způsobem. Pomocí kumulativních počtů budeme generovat podprostory pro rostoucí Krylovovu dimenzi a analyzovat chování konvergence.

from collections import Counter

counts_cumulative = []
for i in range(krylov_dim):
    counter = Counter()
    for d in counts_all[: i + 1]:
        counter.update(d)

    counts = dict(counter)
    counts_cumulative.append(counts)

Pro promítnutí a diagonalizaci hamiltoniánu využíváme funkce z qiskit-addon-sqd. Doplněk nabízí funkce pro promítnutí hamiltoniánů založených na Pauliho řetězcích na podprostor a řeší vlastní hodnoty pomocí SciPy.

from qiskit_addon_sqd.counts import counts_to_arrays
from qiskit_addon_sqd.qubit import solve_qubit

V principu můžeme před generováním podprostoru odfiltrovat bitové řetězce s nesprávným vzorem. Například základní stav antiferomagnetického hamiltoniánu v této lekci má typicky stejný počet spinů „nahoru" a „dolů", tj. počet jedniček v bitovém řetězci by měl být přesně polovina celkového počtu bitů (spinů) v systému. Následující funkce odfiltruje z počtů bitové řetězce s nesprávným počtem jedniček.

# Filters out bitstrings that do not have specified number (`num_ones`) of `1` bits.
def postselect_counts(counts, num_ones):
    filtered_counts = {}
    for bitstring, freq in counts.items():
        if bitstring.count("1") == num_ones:
            filtered_counts[bitstring] = freq

    return filtered_counts

Pomocí bitových řetězců se správným počtem elektronů se spinem nahoru/dolů generujeme podprostory a počítáme vlastní hodnoty pro rostoucí Krylovovu dimenzi. V závislosti na velikosti problému a dostupných klasických zdrojích může být nutné použít podvzorkování (podobně jako v lekci o SQD), aby se dimenze podprostoru udržela v rozumných mezích. Navíc můžeme aplikovat pojem rekonstrukce konfigurace podobně jako v lekci 4. Můžeme vypočítat obsazenost elektronů na místo z rekonstruovaných vlastních stavů a tyto informace použít k opravě bitových řetězců s nesprávným počtem elektronů se spinem nahoru/dolů. Toto ponecháváme jako cvičení pro zájemce.

import numpy as np

num_batches = 10
rand_seed = 0
scipy_kwargs = {"k": 2, "which": "SA"}

ground_state_energies = []
for idx, counts in enumerate(counts_cumulative):
    counts = postselect_counts(counts, num_ones=num_spins // 2)
    bitstring_matrix, probs = counts_to_arrays(counts=counts)

    eigenvals, eigenstates = solve_qubit(
        bitstring_matrix, H_op, verbose=False, **scipy_kwargs
    )
    gs_en = np.min(eigenvals)
    ground_state_energies.append(gs_en)

Dále vykreslíme vypočtenou energii jako funkci Krylovovy dimenze a porovnáme ji s přesnou energií. Přesná energie je vypočtena samostatně pomocí klasické metody hrubé síly. Vidíme, že odhadovaná energie základního stavu konverguje s rostoucí dimenzí Krylovova prostoru. Ačkoli Krylovova dimenze $5$ je omezující, výsledky stále vykazují pozoruhodnou konvergenci, která se očekává zlepší s větší Krylovovou dimenzí [1].

import matplotlib.pyplot as plt

exact_gs_en = -23.934184
plt.plot(
    range(1, krylov_dim + 1),
    ground_state_energies,
    color="blue",
    linestyle="-.",
    label="estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [exact_gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.ylim([-24, -22.50])
plt.title(
    "Estimating Ground state energy with Sample-based Krylov Quantum Diagonalization"
)
plt.show()

Ověř si své znalosti

Přečti si níže uvedené otázky, zamysli se nad odpověďmi a poté klikni na trojúhelníky pro odhalení řešení.

Co by se dalo udělat pro zlepšení konvergence ve výše uvedeném grafu?

Odpověď:

Zvýšit Krylovovu dimenzi. Obecně by bylo možné také zvýšit počet vzorků, ale ten je ve výše uvedeném výpočtu již poměrně vysoký.

Jaké jsou hlavní výhody SKQD oproti (a) SQD a (b) KQD?

Odpověď:

Mohou existovat i jiné platné odpovědi, ale úplné odpovědi by měly obsahovat následující:

(a) SKQD přichází se zárukami konvergence, které SQD nemá. V SQD musíš buď velmi dobře odhadnout svůj ansatz, který má výborný překryv s podporou základního stavu v výpočetní bázi, nebo musíš do výpočtu zavést variační složku pro vzorkování rodiny ansatzů.

(b) SKQD vyžaduje mnohem méně QPU času, protože se vyhýbá nákladnému výpočtu maticových elementů pomocí Hadamardova testu.

5. Shrnutí

• Odhady energie základního stavu pomocí vzorkování Krylovových bázových stavů jsou velmi vhodné pro mřížkové modely včetně spinových systémů, problémů kondenzované hmoty a mřížkových kalibrových teorií. Tento přístup se škáluje mnohem lépe než VQE, protože nevyžaduje optimalizaci přes mnoho parametrů ve variačním ansatzu jako VQE, ani heuristické SQD založené na ansatzu (například chemický problém v předchozí lekci).
  • Pro udržení malé hloubky Circuit je vhodné řešit mřížkové problémy, které jsou vhodné pro hardware před dosažením plné odolnosti vůči chybám.
• SKQD nepřináší problém kvantového měření jako VQE. Není nutné odhadovat skupiny komutujících Pauliho operátorů.
• SKQD je robustní vůči hlučným vzorkům, protože lze použít rutinu post-selekce specifickou pro daný problém (například odfiltrovat bitové řetězce, které se neřídí vzory specifickými pro daný problém) nebo přijmout klasickou diagonalizační režii (tj. diagonalizovat ve větším podprostoru), čímž se efektivně odstraní vliv šumu.

Reference

[1] Jeffery Yu et al., "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization" (2025). arxiv:quant-ph/2501.09702.

[2] Ethan N. Epperly, Lin Lin, and Yuji Nakatsukasa. "A theory of quantum subspace diagonalization". SIAM Journal on Matrix Analysis and Applications 43, 1263–1290 (2022).

[2] N. Hatano and M. Suzuki, "Finding Exponential Product Formulas of Higher Orders" (2005). arXiv:math-ph/0506007.

[4] D. Berry, G. Ahokas, R. Cleve and B. Sanders, "Efficient quantum algorithms for simulating sparse Hamiltonians" (2006). arXiv:quant-ph/0508139.