Diagonalizzazione Quantistica di Krylov basata su Campionamento (SKQD)

Questa lezione sulla diagonalizzazione quantistica di Krylov basata su campionamento (SKQD) combina i metodi spiegati nelle lezioni precedenti. Consiste in un singolo esempio che sfrutta il framework dei pattern di Qiskit:

• Passo 1: Mappare il problema in circuiti quantistici e operatori
• Passo 2: Ottimizzare per l'hardware target
• Passo 3: Eseguire usando i Primitivi di Qiskit
• Passo 4: Post-elaborazione

Un passaggio importante nel metodo di diagonalizzazione quantistica basata su campionamento è generare vettori di qualità per il sottospazio. Nella lezione precedente, abbiamo usato l'ansatz LUCJ per generare vettori di sottospazio per un Hamiltoniano di chimica. In questa lezione, useremo gli stati quantistici di Krylov[1] come discusso nella lezione 2. Prima esamineremo come creare lo spazio di Krylov su un computer quantistico usando operazioni di evoluzione temporale. Poi campioneremo da esso. Proietteremo l'Hamiltoniano del sistema sul sottospazio campionato e lo diagonalizzeremo per stimare l'energia dello stato fondamentale. L'algoritmo converge provabilmente ed efficientemente allo stato fondamentale, sotto le ipotesi descritte nella lezione 2.

0. Lo spazio di Krylov

Ricorda che uno spazio di Krylov $\mathcal{K}^r$ di ordine $r$ è lo spazio generato dai vettori ottenuti moltiplicando potenze crescenti di una matrice $A$, fino a $r-1$, con un vettore di riferimento $\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\}$$

Se la matrice $A$ è l'Hamiltoniano $H$, lo spazio corrispondente è chiamato spazio di Krylov potenziale $\mathcal{K}_P$. Nel caso in cui $A$ sia l'operatore di evoluzione temporale generato dall'Hamiltoniano $U=e^{-iH(dt)}$, lo spazio è detto spazio di Krylov unitario $\mathcal{K}_U$. Il sottospazio di Krylov potenziale non può essere generato direttamente su un computer quantistico poiché $H$ non è un operatore unitario. Possiamo invece usare l'operatore di evoluzione temporale $U = e^{-iH(dt)}$, che può essere dimostrato fornire garanzie di convergenza analoghe a quelle dello spazio di Krylov potenziale. Le potenze di $U$ diventano quindi diversi passi temporali $U^k = e^{-iH(k dt)}$ dove $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. Mappare il problema in circuiti quantistici e operatori

In questa lezione consideriamo l'Hamiltoniano per la catena di spin-1/2 XX-Z antiferromagnetica con $L = 22$ siti con condizione al contorno periodica:

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

Per costruire lo spazio di Krylov, abbiamo bisogno di tre ingredienti principali:

1. Una scelta della dimensione di Krylov ($r$) e del passo temporale ($dt$).
2. Uno stato (vettore) di riferimento iniziale $\vert v \rangle$ con sovrapposizione polinomiale con lo stato target (fondamentale), dove lo stato target è sparso. Questo requisito di sovrapposizione polinomiale è lo stesso dell'algoritmo di stima di fase quantistica.
3. Operatori di evoluzione temporale $U^{k}=e^{-iH(k * dt)}$ ($k = 0, 1, 2, ..., r-1$).

Per un valore scelto di $r$ (e $dt$), creeremo $r$ circuiti quantistici separati e campioneremo da essi. Ogni circuito quantistico è creato unendo la rappresentazione circuitale quantistica dello stato di riferimento e l'operatore di evoluzione temporale per un valore di $k$.

Una dimensione di Krylov maggiore migliora la convergenza dell'energia stimata. In questa lezione impostiamo la dimensione a $5$ per illustrare il trend di convergenza.

Il Rif. [2] ha mostrato che un passo temporale sufficientemente piccolo per KQD è $\pi / \vert \vert H \vert \vert$, e che è preferibile sottostimare questo valore piuttosto che sovrastimarlo. D'altro canto, scegliere $dt$ troppo piccolo porta a un condizionamento peggiore del sottospazio di Krylov, poiché i vettori della base di Krylov differiscono di meno da un passo temporale all'altro. Inoltre, sebbene questa scelta di $dt$ sia provabilmente adeguata per la convergenza di SKQD, in questo contesto basato sul campionamento la scelta ottimale di $dt$ in pratica è un argomento di studio in corso. In questa lezione impostiamo $dt = 0.15$.

Oltre alla dimensione di Krylov e al passo temporale, dobbiamo impostare il numero di passi di Trotter per l'evoluzione temporale. Usare troppo pochi passi porta a errori di Trotterizzazione maggiori, mentre troppi passi portano a circuiti più profondi. In questa lezione impostiamo il numero di passi di Trotter a $6$.

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

Successivamente, dobbiamo scegliere uno stato di riferimento $\vert \psi \rangle$ che abbia una qualche sovrapposizione con lo stato fondamentale. Per questo Hamiltoniano, usiamo lo stato Neel con 1 e 0 alternati $\vert ...101...010...101 \rangle$ come stato di riferimento.

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

Infine, dobbiamo mappare l'operatore di evoluzione temporale in un circuito quantistico. Questo è stato fatto nella lezione 2, ma qui sfrutteremo i metodi di Qiskit, in particolare un metodo chiamato synthesis. Esistono diversi metodi per sintetizzare operatori matematici in circuiti quantistici con gate quantistici. Molte di queste tecniche sono disponibili nel modulo di sintesi di Qiskit. Useremo l'approccio LieTrotter per la sintesi [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. Ottimizzare per l'hardware target

Ora che abbiamo creato i circuiti, possiamo ottimizzarli per un hardware target. Scegliamo una QPU di scala 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")

Ora traspiliamo i circuiti verso il backend target usando un pass manager preimpostato.

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

3. Eseguire sull'hardware target

Dopo aver ottimizzato i circuiti per l'esecuzione su hardware, siamo pronti a eseguirli sull'hardware target e raccogliere campioni per la stima dell'energia dello stato fondamentale.

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-elaborare i risultati

Successivamente, aggreghiamo i conteggi per dimensioni di Krylov crescenti in modo cumulativo. Usando i conteggi cumulativi, costruiremo sottospazi per dimensioni di Krylov crescenti e analizzeremo il comportamento di convergenza.

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)

Per proiettare e diagonalizzare l'Hamiltoniano, usiamo le funzionalità di qiskit-addon-sqd. L'addon offre funzionalità per proiettare Hamiltoniani basati su stringhe di Pauli su un sottospazio e risolvere gli autovalori usando SciPy.

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

In linea di principio, possiamo filtrare le stringhe di bit con pattern errati prima di costruire il sottospazio. Ad esempio, lo stato fondamentale per l'Hamiltoniano antiferromagnetico in questa lezione ha tipicamente un numero uguale di spin "up" e "down", ovvero il numero di "1" nella stringa di bit deve essere esattamente la metà del numero totale di bit (spin) nel sistema. La funzione seguente filtra le stringhe di bit con un numero errato di "1" dai conteggi.

# 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

Usando le stringhe di bit con il numero corretto di elettroni up/down, costruiamo sottospazi e calcoliamo gli autovalori per dimensioni di Krylov crescenti. A seconda della dimensione del problema e delle risorse classiche disponibili, potrebbe essere necessario adottare il sottocampionamento (simile alla lezione su SQD) per mantenere la dimensione del sottospazio sotto controllo. Inoltre, possiamo applicare la nozione di recupero della configurazione simile alla Lezione 4. Possiamo calcolare l'occupazione degli elettroni per sito dagli autostati ricostruiti e usare l'informazione per correggere le stringhe di bit con un numero errato di elettroni up/down. Lasciamo questo come esercizio per i lettori interessati.

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)

Successivamente, tracciamo l'energia calcolata in funzione della dimensione di Krylov e la confrontiamo con l'energia esatta. L'energia esatta è calcolata separatamente usando un metodo classico a forza bruta. Possiamo vedere che l'energia dello stato fondamentale stimata converge con l'aumentare della dimensione dello spazio di Krylov. Sebbene la dimensione di Krylov di $5$ sia limitante, i risultati mostrano comunque una convergenza notevole, che dovrebbe migliorare con una dimensione di Krylov maggiore [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()

Verifica la tua comprensione

Leggi le domande qui sotto, pensa alle tue risposte, poi clicca sui triangoli per scoprire le soluzioni.

Cosa si potrebbe fare per migliorare la convergenza nel grafico sopra?

Risposta:

Aumentare la dimensione di Krylov. In generale, si potrebbe anche aumentare il numero di campioni (shots), ma questo è già piuttosto elevato nel calcolo precedente.

Quali sono i principali vantaggi di SKQD rispetto a (a) SQD e (b) KQD?

Risposta:

Potrebbero esserci altre risposte valide, ma le risposte complete dovrebbero includere quanto segue:

(a) SKQD è dotato di garanzie di convergenza che SQD non ha. In SQD è necessario fare un'ottima stima per l'ansatz con un'eccellente sovrapposizione con il supporto dello stato fondamentale nella base computazionale, oppure introdurre una componente variazionale nel calcolo per campionare una famiglia di ansatz.

(b) SKQD richiede molto meno tempo QPU, perché evita il costoso calcolo degli elementi di matrice tramite il test di Hadamard.

5. Riepilogo

• La stima dell'energia dello stato fondamentale tramite il campionamento degli stati della base di Krylov è molto adatta ai modelli su reticolo, inclusi i sistemi di spin, i problemi di materia condensata e le teorie di gauge su reticolo. Questo approccio scala molto meglio di VQE, poiché non richiede l'ottimizzazione su molti parametri in un ansatz variazionale come in VQE, né in SQD basato su ansatz euristico (ad esempio, il problema di chimica nella lezione precedente).
  • Per mantenere basse le profondità dei circuiti, è saggio affrontare problemi su reticolo adatti all'hardware pre-fault tolerant.
• SKQD non incorre nel problema di misura quantistica come in VQE. Non ci sono gruppi di operatori di Pauli commutanti da stimare.
• SKQD è robusto rispetto a campioni rumorosi, poiché si può usare una routine di post-selezione specifica del problema (ad esempio, filtrare le stringhe di bit che non rispettano pattern specifici del problema) o sostenere un overhead di diagonalizzazione classica (cioè diagonalizzare in un sottospazio più grande) per rimuovere efficacemente l'effetto del rumore.

Riferimenti

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