Sample-based Krylov Quantum Diagonalization (SKQD)

Ta lekcja poświęcona Sample-based Krylov quantum diagonalization (SKQD) łączy metody opisane w poprzednich lekcjach. Składa się z jednego przykładu wykorzystującego framework Qiskit patterns:

• Krok 1: Mapowanie problemu na obwody kwantowe i operatory
• Krok 2: Optymalizacja pod docelowy sprzęt
• Krok 3: Wykonanie z użyciem Qiskit Primitives
• Krok 4: Post-processing

Ważnym krokiem w metodzie sample-based quantum diagonalization jest generowanie wektorów wysokiej jakości dla podprzestrzeni. W poprzedniej lekcji używaliśmy ansatzu LUCJ do generowania wektorów podprzestrzeni dla chemicznego Hamiltonianu. W tej lekcji użyjemy kwantowych stanów Krylov[1], tak jak omówiono w lekcji 2. Najpierw przypomnimy, jak tworzyć przestrzeń Krylov na komputerze kwantowym z użyciem operacji ewolucji czasowej. Następnie będziemy z niej próbkować. Rzutujemy Hamiltonian układu na spróbkowaną podprzestrzeń i dokonamy jego diagonalizacji, aby oszacować energię stanu podstawowego. Algorytm zbieżnie i efektywnie dąży do stanu podstawowego, przy założeniach opisanych w lekcji 2.

0. The Krylov space

Przypomnijmy, że przestrzeń Krylov $\mathcal{K}^r$ rzędu $r$ to przestrzeń rozpiętą przez wektory otrzymane przez mnożenie kolejnych potęg macierzy $A$ — do $r-1$ — przez wektor referencyjny $\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 macierz $A$ jest Hamiltonianem $H$, odpowiadająca mu przestrzeń nazywana jest potęgową przestrzenią Krylov $\mathcal{K}_P$. W przypadku gdy $A$ jest operatorem ewolucji czasowej generowanym przez Hamiltonian $U=e^{-iH(dt)}$, przestrzeń ta nosi nazwę unitarnej przestrzeni Krylov $\mathcal{K}_U$. Potęgowa podprzestrzeń Krylov nie może być generowana bezpośrednio na komputerze kwantowym, ponieważ $H$ nie jest operatorem unitarnym. Zamiast tego możemy użyć operatora ewolucji czasowej $U = e^{-iH(dt)}$, który, jak można pokazać, daje podobne gwarancje zbieżności jak potęgowa przestrzeń Krylov. Potęgi $U$ stają się wówczas różnymi krokami czasowymi $U^k = e^{-iH(k dt)}$, gdzie $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. Map problem to quantum circuits and operators

W tej lekcji rozważamy Hamiltonian dla antysferromagnetycznego łańcucha spinowego XX-Z o spinie 1/2 z $L = 22$ węzłami i periodycznym warunkiem brzegowym:

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

Aby zbudować przestrzeń Krylov, potrzebujemy trzech głównych składników:

1. Wyboru wymiaru Krylov ($r$) i kroku czasowego ($dt$).
2. Stanu początkowego (referencyjnego) (wektor $\vert v \rangle$ powyżej) o wielomianowym pokryciu z docelowym (podstawowym) stanem, gdzie stan docelowy jest rzadki. Wymaganie wielomianowego pokrycia jest takie samo jak w algorytmie kwantowej estymacji fazy.
3. Operatorów ewolucji czasowej $U^{k}=e^{-iH(k * dt)}$ ($k = 0, 1, 2, ..., r-1$).

Dla wybranej wartości $r$ (i $dt$) tworzymy $r$ osobnych obwodów kwantowych i próbkujemy z nich. Każdy obwód kwantowy powstaje przez połączenie kwantowej reprezentacji obwodowej stanu referencyjnego z operatorem ewolucji czasowej dla danej wartości $k$.

Większy wymiar Krylov poprawia zbieżność estymowanej energii. W tej lekcji ustawiamy wymiar na $5$, aby zilustrować trend zbieżności.

Ref [2] wykazał, że wystarczająco mały krok czasowy dla KQD wynosi $\pi / \vert \vert H \vert \vert$, i że lepiej go niedoszacować niż przeszacować. Z drugiej strony, zbyt małe $dt$ prowadzi do gorszego uwarunkowania podprzestrzeni Krylov, ponieważ wektory bazy Krylov mniej różnią się między kolejnymi krokami czasowymi. Ponadto, choć ten wybór $dt$ jest provably wystarczający dla zbieżności SKQD, w tym kontekście opartym na próbkowaniu optymalny dobór $dt$ w praktyce jest nadal przedmiotem badań. W tej lekcji ustawiamy $dt = 0.15$.

Oprócz wymiaru Krylov i kroku czasowego musimy ustalić liczbę kroków Trottera dla ewolucji czasowej. Zbyt mała liczba kroków prowadzi do większych błędów troterizacji, podczas gdy zbyt duża — do głębszych obwodów. W tej lekcji ustawiamy liczbę kroków Trottera na $6$.

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

Następnie musimy wybrać stan referencyjny $\vert \psi \rangle$, który ma pewne pokrycie ze stanem podstawowym. Dla tego Hamiltonianu jako stan referencyjny używamy stanu Neel z naprzemiennymi jedynkami i zerami $\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)

Na koniec musimy odwzorować operator ewolucji czasowej na obwód kwantowy. Zostało to zrobione w lekcji 2, ale tutaj skorzystamy z metod dostępnych w Qiskit, a konkretnie z metody o nazwie synthesis. Istnieją różne metody syntezowania operatorów matematycznych na obwody kwantowe z bramkami kwantowymi. Wiele takich technik dostępnych jest w module syntezy Qiskit. Użyjemy podejścia LieTrotter do syntezy [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. Optimize for target hardware

Mając już gotowe Circuit, możemy je zoptymalizować pod kątem docelowego sprzętu. Wybieramy QPU w skali użytkowej.

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

Teraz wykonujemy transpilację Circuit do docelowego backendu przy użyciu wstępnie zdefiniowanego menedżera przejść.

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

3. Execute on target hardware

Po zoptymalizowaniu Circuit do uruchomienia na sprzęcie jesteśmy gotowi, aby uruchomić je na docelowym sprzęcie i zebrać próbki potrzebne do estymacji energii stanu podstawowego.

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. Post-process results

Następnie agregujemy zliczenia dla rosnących wymiarów Kryłowa w sposób kumulatywny. Korzystając ze skumulowanych zliczeń, będziemy rozpiąć podprzestrzenie dla rosnących wymiarów Kryłowa i analizować zachowanie zbieżności.

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)

Aby rzutować i diagonalizować Hamiltonian, korzystamy z możliwości biblioteki qiskit-addon-sqd. Addon oferuje funkcje do rzutowania Hamiltonianów opartych na ciągach Pauli na podprzestrzeń i rozwiązuje problem własny przy użyciu SciPy.

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

Co do zasady możemy odfiltrować ciągi bitów z nieprawidłowym wzorcem przed rozwinięciem podprzestrzeni. Na przykład stan podstawowy dla Hamiltonianu antysferromagnetycznego w tej lekcji ma zazwyczaj tyle samo spinów „w górę", co „w dół" — liczba jedynek „1" w ciągu bitów powinna być dokładnie połową łącznej liczby bitów (spinów) w układzie. Poniższa funkcja odfiltrowuje ciągi bitów z nieprawidłową liczbą jedynek „1" ze zliczeń.

# 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

Korzystając z ciągów bitów z prawidłową liczbą elektronów „w górę"/„w dół", rozwijamy podprzestrzenie i obliczamy wartości własne dla rosnącego wymiaru Kryłowa. W zależności od rozmiaru problemu i dostępnych zasobów klasycznych możemy potrzebować podpróbkowania (podobnie jak w lekcji o SQD), aby utrzymać wymiar podprzestrzeni w ryzach. Ponadto możemy zastosować koncepcję odtwarzania konfiguracji podobną do tej z lekcji 4. Możemy obliczyć obsadzenie elektronowe na węzeł z odtworzonych stanów własnych i wykorzystać tę informację do korekcji ciągów bitów z nieprawidłową liczbą elektronów „w górę"/„w dół". Pozostawiamy to jako ćwiczenie dla zainteresowanych czytelników.

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)

Następnie wykreślamy obliczoną energię w funkcji wymiaru Kryłowa i porównujemy z energią dokładną. Energia dokładna jest obliczona osobno metodą klasyczną brute-force. Widzimy, że estymowana energia stanu podstawowego zbiega się wraz ze wzrostem wymiaru przestrzeni Kryłowa. Choć wymiar Kryłowa równy $5$ jest ograniczający, wyniki nadal wykazują imponującą zbieżność, która powinna się poprawiać wraz ze wzrostem wymiaru Kryłowa [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()

Check your understanding

Przeczytaj poniższe pytania, zastanów się nad odpowiedziami, a następnie kliknij trójkąty, aby zobaczyć rozwiązania.

Co można zrobić, aby poprawić zbieżność na powyższym wykresie?

Odpowiedź:

Zwiększyć wymiar Kryłowa. Ogólnie rzecz biorąc, można by też zwiększyć liczbę pomiarów (shots), ale w powyższych obliczeniach jest ona już dość wysoka.

Jakie są główne zalety SKQD w porównaniu z (a) SQD i (b) KQD?

Odpowiedź:

Mogą istnieć inne poprawne odpowiedzi, ale kompletna odpowiedź powinna zawierać następujące punkty:

(a) SKQD posiada gwarancje zbieżności, których SQD nie daje. W SQD albo trzeba bardzo dobrze dobrać ansatz, który ma doskonałe pokrycie z nośnikiem stanu podstawowego w bazie obliczeniowej, albo trzeba wprowadzić do obliczeń element wariacyjny, aby próbkować rodzinę ansatze.

(b) SKQD wymaga znacznie mniej czasu QPU, ponieważ unika kosztownego obliczania elementów macierzowych za pomocą testu Hadamarda.

5. Summary

• Estymacja energii stanu podstawowego poprzez próbkowanie stanów bazy Kryłowa doskonale nadaje się do modeli sieciowych, w tym układów spinowych, problemów fizyki materii skondensowanej oraz sieciowych teorii cechowania. To podejście skaluje się znacznie lepiej niż VQE, ponieważ nie wymaga optymalizacji po wielu parametrach wariacyjnego ansatzu jak w VQE ani stosowania heurystycznego ansatzu w SQD (na przykład problemu chemicznego z poprzedniej lekcji).
  • Aby utrzymać głębokość Circuit na niskim poziomie, warto podejmować problemy sieciowe podatne na sprzęt przed osiągnięciem tolerancji na błędy.
• SKQD nie napotyka kwantowego problemu pomiarowego jak w VQE. Nie ma grup komutujących operatorów Pauli, które trzeba estymować.
• SKQD jest odporne na zaszumione próbki, ponieważ można zastosować selekcję wstępną specyficzną dla problemu (na przykład odfiltrowanie ciągów bitów niespełniających wzorców charakterystycznych dla problemu) lub ponieść klasyczny narzut diagonalizacyjny (czyli diagonalizować w większej podprzestrzeni), aby skutecznie wyeliminować wpływ szumu.

References

[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).

[3] 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.