Diagonalizzazione quantistica di Krylov di Hamiltoniane su reticolo

Stima di utilizzo: 20 minuti su un Heron r2 (NOTA: Questa è solo una stima. Il vostro tempo di esecuzione potrebbe variare.)

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

# This cell is hidden from users – it disables some lint rules
# ruff: noqa: E402 E722 F601

Contesto

Questo tutorial dimostra come implementare l'Algoritmo di Diagonalizzazione Quantistica di Krylov (KQD) nel contesto dei pattern Qiskit. Imparerete prima la teoria alla base dell'algoritmo e poi vedrete una dimostrazione della sua esecuzione su una QPU.

In diverse discipline, siamo interessati a conoscere le proprietà dello stato fondamentale dei sistemi quantistici. Gli esempi includono la comprensione della natura fondamentale delle particelle e delle forze, la previsione e la comprensione del comportamento di materiali complessi e la comprensione delle interazioni e reazioni biochimiche. A causa della crescita esponenziale dello spazio di Hilbert e della correlazione che si manifesta nei sistemi entangled, gli algoritmi classici faticano a risolvere questo problema per sistemi quantistici di dimensioni crescenti. A un'estremità dello spettro vi è l'approccio esistente che sfrutta l'hardware quantistico concentrandosi sui metodi quantistici variazionali (ad esempio, variational quantum eigensolver). Queste tecniche affrontano sfide con i dispositivi attuali a causa dell'elevato numero di chiamate di funzione richieste nel processo di ottimizzazione, che aggiunge un grande overhead in risorse una volta introdotte tecniche avanzate di mitigazione degli errori, limitando così la loro efficacia a sistemi di piccole dimensioni. All'altra estremità dello spettro, vi sono metodi quantistici fault-tolerant con garanzie di prestazioni (ad esempio, quantum phase estimation), che richiedono circuiti profondi che possono essere eseguiti solo su un dispositivo fault-tolerant. Per questi motivi, introduciamo qui un algoritmo quantistico basato su metodi di sottospazio (come descritto in questo review paper), l'algoritmo di diagonalizzazione quantistica di Krylov (KQD). Questo algoritmo funziona bene su larga scala [1] sull'hardware quantistico esistente, condivide garanzie di prestazioni simili alla stima di fase, è compatibile con tecniche avanzate di mitigazione degli errori e potrebbe fornire risultati che sono classicamente inaccessibili.

Requisiti

Prima di iniziare questo tutorial, assicuratevi di avere installato quanto segue:

• Qiskit SDK v2.0 o successivo, con supporto per la visualizzazione
• Qiskit Runtime v0.22 o successivo ( pip install qiskit-ibm-runtime )

Configurazione

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import itertools as it
import copy
from sympy import Matrix
import warnings

warnings.filterwarnings("ignore")

from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState
from qiskit.circuit import Parameter, IfElseOp
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter
from qiskit.transpiler import Target, CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    EstimatorV2 as Estimator,
)

def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in the Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of the Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the regularized generalized eigenvalue problem

    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array(
        [vec for val, vec in zip(s_vals, s_vecs) if val > threshold]
    )
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]

def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]  # list all of the possible sets of n_exc indices of 1s in n_exc-particle states

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

Step 1: Mappare gli input classici a un problema quantistico

Lo spazio di Krylov

Lo spazio di Krylov $\mathcal{K}^r$ di ordine $r$ è lo spazio generato dai vettori ottenuti moltiplicando potenze più elevate 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'Hamiltoniana $H$, faremo riferimento allo spazio corrispondente come spazio di Krylov di potenza $\mathcal{K}_P$. Nel caso in cui $A$ sia l'operatore di evoluzione temporale generato dall'Hamiltoniana $U=e^{-iHt}$, faremo riferimento allo spazio come spazio di Krylov unitario $\mathcal{K}_U$. Il sottospazio di Krylov di potenza che utilizziamo classicamente non può essere generato direttamente su un computer quantistico poiché $H$ non è un operatore unitario. Invece, possiamo usare l'operatore di evoluzione temporale $U = e^{-iHt}$ che può essere dimostrato fornire garanzie di convergenza simili al metodo delle potenze. Le potenze di $U$ diventano quindi diversi passi temporali $U^k = e^{-iH(kt)}$.

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

Consultate l'Appendice per una derivazione dettagliata di come lo spazio di Krylov unitario permetta di rappresentare accuratamente gli autostati a bassa energia.

Algoritmo di diagonalizzazione quantistica di Krylov

Data un'Hamiltoniana $H$ che desideriamo diagonalizzare, consideriamo prima il corrispondente spazio di Krylov unitario $\mathcal{K}_U$. L'obiettivo è trovare una rappresentazione compatta dell'Hamiltoniana in $\mathcal{K}_U$, che chiameremo $\tilde{H}$. Gli elementi di matrice di $\tilde{H}$, la proiezione dell'Hamiltoniana nello spazio di Krylov, possono essere calcolati calcolando i seguenti valori di aspettazione

$$\tilde{H}_{mn} = \langle \psi_m \vert H \vert \psi_n \rangle =$$ $$= \langle \psi \vert e^{i H t_m} H e^{-i H t_n} \vert \psi \rangle$$ $$= \langle \psi \vert e^{i H m dt} H e^{-i H n dt} \vert \psi \rangle$$

Dove $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ sono i vettori dello spazio di Krylov unitario e $t_n = n dt$ sono i multipli del passo temporale $dt$ scelto. Su un computer quantistico, il calcolo di ciascun elemento di matrice può essere effettuato con qualsiasi algoritmo che permetta di ottenere la sovrapposizione tra stati quantistici. Questo tutorial si concentra sul test di Hadamard. Dato che $\mathcal{K}_U$ ha dimensione $r$, l'Hamiltoniana proiettata nel sottospazio avrà dimensioni $r \times r$. Con $r$ sufficientemente piccolo (generalmente $r<<100$ è sufficiente per ottenere la convergenza delle stime delle autoenegie) possiamo quindi facilmente diagonalizzare l'Hamiltoniana proiettata $\tilde{H}$. Tuttavia, non possiamo diagonalizzare direttamente $\tilde{H}$ a causa della non-ortogonalità dei vettori dello spazio di Krylov. Dovremo misurare le loro sovrapposizioni e costruire una matrice $\tilde{S}$

$$\tilde{S}_{mn} = \langle \psi_m \vert \psi_n \rangle$$

Questo ci permette di risolvere il problema agli autovalori in uno spazio non ortogonale (chiamato anche problema agli autovalori generalizzato)

$$\tilde{H} \ \vec{c} = E \ \tilde{S} \ \vec{c}$$

Si possono quindi ottenere stime degli autovalori e degli autostati di $H$ osservando quelli di $\tilde{H}$. Ad esempio, la stima dell'energia dello stato fondamentale si ottiene prendendo l'autovalore più piccolo $c$ e lo stato fondamentale dal corrispondente autovettore $\vec{c}$. I coefficienti in $\vec{c}$ determinano il contributo dei diversi vettori che generano $\mathcal{K}_U$.

La figura mostra una rappresentazione circuitale del test di Hadamard modificato, un metodo utilizzato per calcolare la sovrapposizione tra diversi stati quantistici. Per ogni elemento di matrice $\tilde{H}_{i,j}$, viene eseguito un test di Hadamard tra lo stato $\vert \psi_i \rangle$, $\vert \psi_j \rangle$. Questo è evidenziato nella figura dallo schema di colori per gli elementi di matrice e le corrispondenti operazioni $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$. Pertanto, è necessario un insieme di test di Hadamard per tutte le possibili combinazioni di vettori dello spazio di Krylov per calcolare tutti gli elementi di matrice dell'Hamiltoniana proiettata $\tilde{H}$. Il filo superiore nel circuito del test di Hadamard è un qubit ancilla che viene misurato nella base X o Y, il suo valore di aspettazione determina il valore della sovrapposizione tra gli stati. Il filo inferiore rappresenta tutti i qubit dell'Hamiltoniana del sistema. L'operazione $\text{Prep} \; \psi_i$ prepara il qubit del sistema nello stato $\vert \psi_i \rangle$ controllato dallo stato del qubit ancilla (similmente per $\text{Prep} \; \psi_j$) e l'operazione $P$ rappresenta la decomposizione di Pauli dell'Hamiltoniana del sistema $H = \sum_i P_i$. Una derivazione più dettagliata delle operazioni calcolate dal test di Hadamard è fornita di seguito.

Definire l'Hamiltoniana

Consideriamo l'Hamiltoniana di Heisenberg per $N$ qubit su una catena lineare: $H= \sum_{i,j}^N X_i X_j + Y_i Y_j - J Z_i Z_j$

# Define problem Hamiltonian.
n_qubits = 30
J = 1  # coupling strength for ZZ interaction

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [(term, J) if term.count("Z") == 2 else (term, 1) for term in H_int]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)

[('ZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 1), ('XXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 1), ('YYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 1)]

Impostare i parametri per l'algoritmo

Scegliamo euristicamente un valore per il passo temporale dt (basato su limiti superiori della norma dell'Hamiltoniana). Il riferimento [2] ha dimostrato che un passo temporale sufficientemente piccolo è $\pi/\vert \vert H \vert \vert$, e che è preferibile fino a un certo punto sottostimare questo valore piuttosto che sovrastimarlo, poiché sovrastimarlo può permettere ai contributi degli stati ad alta energia di corrompere anche lo stato ottimale nello spazio di Krylov. D'altra parte, scegliere $dt$ troppo piccolo porta a un peggiore condizionamento del sottospazio di Krylov, poiché i vettori di base di Krylov differiscono meno da un passo temporale all'altro.

# Get Hamiltonian restricted to single-particle states
single_particle_H = np.zeros((n_qubits, n_qubits))
for i in range(n_qubits):
    for j in range(i + 1):
        for p, coeff in H_op.to_list():
            p_x = Pauli(p).x
            p_z = Pauli(p).z
            if all(
                p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)
            ):
                sgn = (
                    (-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))
                ) * ((-1) ** p_z[i])
            else:
                sgn = 0
            single_particle_H[i, j] += sgn * coeff
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        single_particle_H[i, j] = np.conj(single_particle_H[j, i])

# Set dt according to spectral norm
dt = np.pi / np.linalg.norm(single_particle_H, ord=2)
dt

np.float64(0.10833078115826875)

E impostiamo altri parametri dell'algoritmo. Per il bene di questo tutorial, ci limiteremo a utilizzare uno spazio di Krylov con solo cinque dimensioni, il che è abbastanza limitante.

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

Preparazione dello stato

Scegliete uno stato di riferimento $\vert \psi \rangle$ che abbia una certa sovrapposizione con lo stato fondamentale. Per questa Hamiltoniana, utilizziamo uno stato con un'eccitazione nel qubit centrale $\vert 00..010...00 \rangle$ come nostro stato di riferimento.

qc_state_prep = QuantumCircuit(n_qubits)
qc_state_prep.x(int(n_qubits / 2) + 1)
qc_state_prep.draw("mpl", scale=0.5)

Evoluzione temporale

Possiamo realizzare l'operatore di evoluzione temporale generato da una data Hamiltoniana: $U=e^{-iHt}$ tramite l'approssimazione di Lie-Trotter.

t = Parameter("t")

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)
)

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

<qiskit.circuit.instructionset.InstructionSet at 0x11eef9be0>

Test di Hadamard

$$\begin{equation*} |0\rangle|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle + |1\rangle \Big)|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle|0\rangle^N+|1\rangle |\psi_i\rangle\Big) \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle |0\rangle^N+|1\rangle P |\psi_i\rangle\Big) \quad\longrightarrow\quad\frac{1}{\sqrt{2}}\Big(|0\rangle |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{equation*}$$

Dove $P$ è uno dei termini nella decomposizione dell'Hamiltoniana $H=\sum P$ e $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ sono operazioni controllate che preparano $|\psi_i\rangle$, $|\psi_j\rangle$ vettori dello spazio di Krylov unitario, con $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Per misurare $X$, applicare prima $H$...

$$\begin{equation*} \longrightarrow\quad\frac{1}{2}|0\rangle\Big( |\psi_j\rangle + P|\psi_i\rangle\Big) + \frac{1}{2}|1\rangle\Big(|\psi_j\rangle - P|\psi_i\rangle\Big) \end{equation*}$$

... poi misurare:

$$\begin{equation*} \begin{split} \Rightarrow\quad\langle X\rangle &= \frac{1}{4}\Bigg(\Big\|| \psi_j\rangle + P|\psi_i\rangle \Big\|^2-\Big\||\psi_j\rangle - P|\psi_i\rangle\Big\|^2\Bigg) \\ &= \text{Re}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{split} \end{equation*}$$

Dall'identità $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Analogamente, misurare $Y$ produce

$$\begin{equation*} \langle Y\rangle = \text{Im}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{equation*}$$

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)

## Create the time-evo op dagger circuit
evol_gate_d = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)
evol_gate_d = evol_gate_d.inverse()

# Put pieces together
qc_reg = QuantumRegister(n_qubits)
qc_temp = QuantumCircuit(qc_reg)
qc_temp.compose(qc_state_prep, inplace=True)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate, qargs=qc_reg)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate_d, qargs=qc_reg)
qc_temp.compose(qc_state_prep.inverse(), inplace=True)

# Create controlled version of the circuit
controlled_U = qc_temp.to_gate().control(1)

# Create hadamard test circuit for real part
qr = QuantumRegister(n_qubits + 1)
qc_real = QuantumCircuit(qr)
qc_real.h(0)
qc_real.append(controlled_U, list(range(n_qubits + 1)))
qc_real.h(0)

print(
    "Circuit for calculating the real part of the overlap in S via Hadamard test"
)
qc_real.draw("mpl", fold=-1, scale=0.5)

Circuit for calculating the real part of the overlap in S via Hadamard test

Il circuito del test di Hadamard può essere un circuito profondo una volta decomposto in porte native (che aumenterà ancora di più se teniamo conto della topologia del dispositivo)

print(
    "Number of layers of 2Q operations",
    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),
)

Number of layers of 2Q operations 112753

Step 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Test di Hadamard efficiente

Possiamo ottimizzare i circuiti profondi per il test di Hadamard che abbiamo ottenuto introducendo alcune approssimazioni e facendo affidamento su alcune assunzioni sull'Hamiltoniana del modello. Ad esempio, considerate il seguente circuito per il test di Hadamard:

Supponiamo di poter calcolare classicamente $E_0$, l'autovalore di $|0\rangle^N$ sotto l'Hamiltoniana $H$. Questa condizione è soddisfatta quando l'Hamiltoniana preserva la simmetria U(1). Sebbene questa possa sembrare un'assunzione forte, ci sono molti casi in cui è sicuro assumere che esista uno stato di vuoto (in questo caso corrisponde allo stato $|0\rangle^N$) che non è influenzato dall'azione dell'Hamiltoniana. Questo è vero ad esempio per le Hamiltoniane di chimica che descrivono molecole stabili (dove il numero di elettroni è conservato). Dato che il gate $\text{Prep} \; \psi$ prepara lo stato di riferimento desiderato $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, ad esempio, per preparare lo stato HF per la chimica $\text{Prep} \; \psi$ sarebbe un prodotto di NOT a singolo qubit, quindi controlled-$\text{Prep} \; \psi$ è semplicemente un prodotto di CNOT. Allora il circuito sopra implementa il seguente stato prima della misura:

$$\begin{equation} \begin{split} \ket{0} \ket{0}^N\xrightarrow{H}&\frac{1}{\sqrt{2}} \left( \ket{0}\ket{0}^N+ \ket{1} \ket{0}^N \right)\\ \xrightarrow{\text{1-ctrl-init}}&\frac{1}{\sqrt{2}}\left(|0\rangle|0\rangle^N+|1\rangle|\psi\rangle\right)\\ \xrightarrow{U}&\frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0}\ket{0}^N+\ket{1} U\ket{\psi}\right)\\ \xrightarrow{\text{0-ctrl-init}}&\frac{1}{\sqrt{2}} \left( e^{i\phi}\ket{0} \ket{\psi} +\ket{1} U\ket{\psi} \right)\\ =&\frac{1}{2} \left( \ket{+}\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) +\ket{-}\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \right)\\ =&\frac{1}{2} \left( \ket{+i}\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) +\ket{-i}\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \right) \end{split} \end{equation}$$

dove abbiamo usato la fase shift classicamente simulabile $U\ket{0}^N = e^{i\phi}\ket{0}^N$ nella terza riga. Pertanto i valori di aspettazione si ottengono come

$$\begin{equation} \begin{split} \langle X\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \Big)\\ &=\text{Re}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right], \end{split} \end{equation}$$ $$\begin{equation} \begin{split} \langle Y\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \Big)\\ &=\text{Im}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right]. \end{split} \end{equation}$$

Utilizzando queste assunzioni siamo stati in grado di scrivere i valori di aspettazione degli operatori di interesse con meno operazioni controllate. Infatti, dobbiamo solo implementare la preparazione di stato controllata $\text{Prep} \; \psi$ e non le evoluzioni temporali controllate. Riformulare il nostro calcolo come sopra ci permetterà di ridurre notevolmente la profondità dei circuiti risultanti.

Decomporre l'operatore di evoluzione temporale con la decomposizione di Trotter

Invece di implementare l'operatore di evoluzione temporale esattamente, possiamo usare la decomposizione di Trotter per implementarne un'approssimazione. Ripetere più volte una decomposizione di Trotter di un certo ordine ci dà un'ulteriore riduzione dell'errore introdotto dall'approssimazione. Nel seguito, costruiamo direttamente l'implementazione di Trotter nel modo più efficiente per il grafo di interazione dell'Hamiltoniana che stiamo considerando (solo interazioni tra vicini più prossimi). In pratica inseriamo rotazioni di Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ con un angolo parametrizzato $t$ che corrispondono all'implementazione approssimata di $e^{-i (XX + YY + ZZ) t}$. Data la differenza nella definizione delle rotazioni di Pauli e dell'evoluzione temporale che stiamo cercando di implementare, dovremo usare il parametro $2*dt$ per ottenere un'evoluzione temporale di $dt$. Inoltre, invertiamo l'ordine delle operazioni per i numeri dispari di ripetizioni degli step di Trotter, il che è funzionalmente equivalente ma consente di sintetizzare operazioni adiacenti in un singolo unitario $SU(2)$. Questo fornisce un circuito molto più superficiale di quello ottenuto usando la funzionalità generica PauliEvolutionGate().

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(t, 0, 1)
Rxyz_circ.ryy(t, 0, 1)
Rxyz_circ.rzz(t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="RXX+YY+ZZ")

interaction_list = [
    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],
    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],
]  # linear chain

qr = QuantumRegister(n_qubits)
trotter_step_circ = QuantumCircuit(qr)
for i, color in enumerate(interaction_list):
    for interaction in color:
        trotter_step_circ.append(Rxyz_instr, interaction)
    if i < len(interaction_list) - 1:
        trotter_step_circ.barrier()
reverse_trotter_step_circ = trotter_step_circ.reverse_ops()

qc_evol = QuantumCircuit(qr)
for step in range(num_trotter_steps):
    if step % 2 == 0:
        qc_evol = qc_evol.compose(trotter_step_circ)
    else:
        qc_evol = qc_evol.compose(reverse_trotter_step_circ)

qc_evol.decompose().draw("mpl", fold=-1, scale=0.5)

Utilizzare un circuito ottimizzato per la preparazione dello stato

control = 0
excitation = int(n_qubits / 2) + 1
controlled_state_prep = QuantumCircuit(n_qubits + 1)
controlled_state_prep.cx(control, excitation)
controlled_state_prep.draw("mpl", fold=-1, scale=0.5)

Circuiti template per calcolare gli elementi di matrice di $\tilde{S}$ e $\tilde{H}$ tramite il test di Hadamard

L'unica differenza tra i circuiti utilizzati nel test di Hadamard sarà la fase nell'operatore di evoluzione temporale e gli osservabili misurati. Pertanto possiamo preparare un circuito template che rappresenta il circuito generico per il test di Hadamard, con segnaposto per i gate che dipendono dall'operatore di evoluzione temporale.

# Parameters for the template circuits
parameters = []
for idx in range(1, krylov_dim):
    parameters.append(2 * dt_circ * (idx))

# Create modified hadamard test circuit
qr = QuantumRegister(n_qubits + 1)
qc = QuantumCircuit(qr)
qc.h(0)
qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)
qc.barrier()
qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)
qc.barrier()
qc.x(0)
qc.compose(
    controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True
)
qc.x(0)

qc.decompose().draw("mpl", fold=-1)

print(
    "The optimized circuit has 2Q gates depth: ",
    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),
)

The optimized circuit has 2Q gates depth:  74

Abbiamo ridotto considerevolmente la profondità del test di Hadamard con una combinazione di approssimazione di Trotter e unitari non controllati

Step 3: Eseguire utilizzando le primitive di Qiskit

Istanziare il backend e impostare i parametri di runtime

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
if (
    "if_else" not in backend.target.operation_names
):  # Needed as "op_name" could be "if_else"
    backend.target.add_instruction(IfElseOp, name="if_else")
print(backend.name)

Transpilare su una QPU

Prima, scegliamo sottoinsiemi della mappa di accoppiamento con qubit che hanno "buone" prestazioni (dove "buone" è abbastanza arbitrario qui, vogliamo principalmente evitare qubit che si comportano davvero male) e creiamo un nuovo target per la transpilazione

target = backend.target
cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    t2 = target.qubit_properties[q].t2 * 1e6
    if meas_err > 0.02 or t2 < 100:
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except:
                    continue

for q in cmap_list:
    op_name = list(target.operation_names_for_qargs(q))[0]
    twoq_gate_err = target[f"{op_name}"][q].error
    if twoq_gate_err > 0.005:
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q)
                except:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)
cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates,
    coupling_map=cust_cmap,
)

Poi transpiliamo il circuito virtuale al miglior layout fisico in questo nuovo target

basis_gates = list(target.operation_names)
pm = generate_preset_pass_manager(
    optimization_level=3,
    target=cust_target,
    basis_gates=basis_gates,
)

qc_trans = pm.run(qc)

print("depth", qc_trans.depth(lambda x: x[0].num_qubits == 2))
print("num 2q ops", qc_trans.count_ops())
print(
    "physical qubits",
    sorted(
        [
            idx
            for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()
            if qb._register.name != "ancilla"
        ]
    ),
)

depth 52
num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])
physical qubits [91, 92, 93, 94, 95, 98, 99, 108, 109, 110, 111, 113, 114, 115, 119, 127, 132, 133, 134, 135, 137, 139, 147, 148, 149, 150, 151, 152, 153, 154, 155]

Creare PUB per l'esecuzione con Estimator

# Define observables to measure for S
observable_S_real = "I" * (n_qubits) + "X"
observable_S_imag = "I" * (n_qubits) + "Y"

observable_op_real = SparsePauliOp(
    observable_S_real
)  # define a sparse pauli operator for the observable
observable_op_imag = SparsePauliOp(observable_S_imag)

layout = qc_trans.layout  # get layout of transpiled circuit
observable_op_real = observable_op_real.apply_layout(
    layout
)  # apply physical layout to the observable
observable_op_imag = observable_op_imag.apply_layout(layout)
observable_S_real = (
    observable_op_real.paulis.to_labels()
)  # get the label of the physical observable
observable_S_imag = observable_op_imag.paulis.to_labels()

observables_S = [[observable_S_real], [observable_S_imag]]

# Define observables to measure for H
# Hamiltonian terms to measure
observable_list = []
for pauli, coeff in zip(H_op.paulis, H_op.coeffs):
    # print(pauli)
    observable_H_real = pauli[::-1].to_label() + "X"
    observable_H_imag = pauli[::-1].to_label() + "Y"
    observable_list.append([observable_H_real])
    observable_list.append([observable_H_imag])

layout = qc_trans.layout

observable_trans_list = []
for observable in observable_list:
    observable_op = SparsePauliOp(observable)
    observable_op = observable_op.apply_layout(layout)
    observable_trans_list.append([observable_op.paulis.to_labels()])

observables_H = observable_trans_list

# Define a sweep over parameter values
params = np.vstack(parameters).T

# Estimate the expectation value for all combinations of
# observables and parameter values, where the pub result will have
# shape (# observables, # parameter values).
pub = (qc_trans, observables_S + observables_H, params)

Eseguire i circuiti

I circuiti per $t=0$ sono calcolabili classicamente

qc_cliff = qc.assign_parameters({t: 0})

# Get expectation values from experiment
S_expval_real = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "X")
)
S_expval_imag = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "Y")
)

# Get expectation values
S_expval = S_expval_real + 1j * S_expval_imag

H_expval = 0
for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Get expectation values from experiment
    expval_real = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "X")
    )
    expval_imag = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "Y")
    )
    expval = expval_real + 1j * expval_imag

    # Fill-in matrix elements
    H_expval += coeff * expval

print(H_expval)

(25+0j)

Eseguire i circuiti per $S$ e $\tilde{H}$ con l'Estimator

# Experiment options
num_randomizations = 300
num_randomizations_learning = 30
shots_per_randomization = 100
noise_factors = [1, 1.2, 1.4]
learning_pair_depths = [0, 4, 24, 48]

experimental_opts = {}
experimental_opts["resilience"] = {
    "measure_mitigation": True,
    "measure_noise_learning": {
        "num_randomizations": num_randomizations_learning,
        "shots_per_randomization": shots_per_randomization,
    },
    "zne_mitigation": True,
    "zne": {"noise_factors": noise_factors},
    "layer_noise_learning": {
        "max_layers_to_learn": 10,
        "layer_pair_depths": learning_pair_depths,
        "shots_per_randomization": shots_per_randomization,
        "num_randomizations": num_randomizations_learning,
    },
    "zne": {
        "amplifier": "pea",
        "extrapolated_noise_factors": [0] + noise_factors,
    },
}
experimental_opts["twirling"] = {
    "num_randomizations": num_randomizations,
    "shots_per_randomization": shots_per_randomization,
    "strategy": "all",
}

estimator = Estimator(mode=backend, options=experimental_opts)

job = estimator.run([pub])

Step 4: Post-elaborazione e restituzione del risultato in formato classico desiderato

results = job.result()[0]

Calcolare le matrici Hamiltoniana Effettiva e di Overlap

Prima calcoliamo la fase accumulata dallo stato $\vert 0 \rangle$ durante l'evoluzione temporale non controllata

prefactors = [
    np.exp(-1j * sum([c for p, c in H_op.to_list() if "Z" in p]) * i * dt)
    for i in range(1, krylov_dim)
]

Una volta ottenuti i risultati delle esecuzioni del circuito possiamo post-elaborare i dati per calcolare gli elementi di matrice di $S$

# Assemble S, the overlap matrix of dimension D:
S_first_row = np.zeros(krylov_dim, dtype=complex)
S_first_row[0] = 1 + 0j

# Add in ancilla-only measurements:
for i in range(krylov_dim - 1):
    # Get expectation values from experiment
    expval_real = results.data.evs[0][0][
        i
    ]  # automatic extrapolated evs if ZNE is used
    expval_imag = results.data.evs[1][0][
        i
    ]  # automatic extrapolated evs if ZNE is used

    # Get expectation values
    expval = expval_real + 1j * expval_imag
    S_first_row[i + 1] += prefactors[i] * expval

S_first_row_list = S_first_row.tolist()  # for saving purposes

S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        S_circ[j, i] = S_first_row[i - j]
    else:
        S_circ[j, i] = np.conj(S_first_row[j - i])

Matrix(S_circ)

$$\displaystyle \left[\begin{matrix}1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i & 0.0012070853532697 + 0.312052218182462 i\\-0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i\\0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i\\-0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i\\0.0012070853532697 - 0.312052218182462 i & -0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0\end{matrix}\right]$$

E gli elementi di matrice di $\tilde{H}$

# Assemble S, the overlap matrix of dimension D:
H_first_row = np.zeros(krylov_dim, dtype=complex)
H_first_row[0] = H_expval

for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Add in ancilla-only measurements:
    for i in range(krylov_dim - 1):
        # Get expectation values from experiment
        expval_real = results.data.evs[2 + 2 * obs_idx][0][
            i
        ]  # automatic extrapolated evs if ZNE is used
        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][
            i
        ]  # automatic extrapolated evs if ZNE is used

        # Get expectation values
        expval = expval_real + 1j * expval_imag
        H_first_row[i + 1] += prefactors[i] * coeff * expval

H_first_row_list = H_first_row.tolist()

H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        H_eff_circ[j, i] = H_first_row[i - j]
    else:
        H_eff_circ[j, i] = np.conj(H_first_row[j - i])

Matrix(H_eff_circ)

$$\displaystyle \left[\begin{matrix}25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i & 1.98818301405581 + 5.8897614762563 i\\-14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i\\10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i\\-5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i\\1.98818301405581 - 5.8897614762563 i & -5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0\end{matrix}\right]$$

Infine, possiamo risolvere il problema agli autovalori generalizzato per $\tilde{H}$:

$\tilde{H} \vec{c} = c S \vec{c}$

e ottenere una stima dell'energia dello stato fondamentale $c_{min}$

gnd_en_circ_est_list = []
for d in range(1, krylov_dim + 1):
    # Solve generalized eigenvalue problem for different size of the Krylov space
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-1
    )
    gnd_en_circ_est_list.append(gnd_en_circ_est)
    print("The estimated ground state energy is: ", gnd_en_circ_est)

The estimated ground state energy is:  25.0
The estimated ground state energy is:  22.572154819954875
The estimated ground state energy is:  21.691509219286587
The estimated ground state energy is:  21.23882298756386
The estimated ground state energy is:  20.965499325470294

Per un settore a singola particella, possiamo calcolare efficientemente lo stato fondamentale di questo settore dell'Hamiltoniana in modo classico

gs_en = single_particle_gs(H_op, n_qubits)

n_sys_qubits 30
n_exc 1 , subspace dimension 31
single particle ground state energy:  21.021912418526906

plt.plot(
    range(1, krylov_dim + 1),
    gnd_en_circ_est_list,
    color="blue",
    linestyle="-.",
    label="KQD estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [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.title(
    "Estimating Ground state energy with Krylov Quantum Diagonalization"
)
plt.show()

Appendice: Sottospazio di Krylov da evoluzioni temporali reali

Lo spazio di Krylov unitario è definito come

$$\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ |\psi\rangle, e^{-iH\,dt} |\psi\rangle, \dots, e^{-irH\,dt} |\psi\rangle \right\}$$

per un certo passo temporale $dt$ che determineremo più avanti. Supponiamo temporaneamente che $r$ sia pari: allora definiamo $d=r/2$. Notate che quando proiettiamo l'Hamiltoniana nello spazio di Krylov sopra, essa è indistinguibile dallo spazio di Krylov

$$\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ e^{i\,d\,H\,dt}|\psi\rangle, e^{i(d-1)H\,dt} |\psi\rangle, \dots, e^{-i(d-1)H\,dt} |\psi\rangle, e^{-i\,d\,H\,dt} |\psi\rangle \right\},$$

cioè dove tutte le evoluzioni temporali sono spostate indietro di $d$ passi temporali. Il motivo per cui è indistinguibile è che gli elementi di matrice

$$\tilde{H}_{j,k} = \langle\psi|e^{i\,j\,H\,dt}He^{-i\,k\,H\,dt}|\psi\rangle=\langle\psi|He^{i(j-k)H\,dt}|\psi\rangle$$

sono invarianti sotto spostamenti complessivi del tempo di evoluzione, poiché le evoluzioni temporali commutano con l'Hamiltoniana. Per $r$ dispari, possiamo usare l'analisi per $r-1$.

Vogliamo mostrare che da qualche parte in questo spazio di Krylov è garantito che ci sia uno stato a bassa energia. Lo facciamo attraverso il seguente risultato, che è derivato dal Teorema 3.1 in [3]:

Affermazione 1: esiste una funzione $f$ tale che per energie $E$ nell'intervallo spettrale dell'Hamiltoniana (cioè tra l'energia dello stato fondamentale e l'energia massima)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ per tutti i valori di $E$ che si trovano a distanza $\ge\delta$ da $E_0$, cioè è esponenzialmente soppressa
3. $f(E)$ è una combinazione lineare di $e^{ijE\,dt}$ per $j=-d,-d+1,...,d-1,d$

Forniamo una dimostrazione qui sotto, ma può essere tranquillamente saltata a meno che non si voglia comprendere l'argomento completo e rigoroso. Per ora ci concentriamo sulle implicazioni dell'affermazione sopra. Dalla proprietà 3 sopra, possiamo vedere che lo spazio di Krylov spostato sopra contiene lo stato $f(H)|\psi\rangle$. Questo è il nostro stato a bassa energia. Per capire perché, scriviamo $|\psi\rangle$ nella base degli autostati di energia:

$$|\psi\rangle = \sum_{k=0}^{N}\gamma_k|E_k\rangle,$$

dove $|E_k\rangle$ è il k-esimo autostato di energia e $\gamma_k$ è la sua ampiezza nello stato iniziale $|\psi\rangle$. Espresso in questi termini, $f(H)|\psi\rangle$ è dato da

$$f(H)|\psi\rangle = \sum_{k=0}^{N}\gamma_kf(E_k)|E_k\rangle,$$

usando il fatto che possiamo sostituire $H$ con $E_k$ quando agisce sull'autostato $|E_k\rangle$. L'errore di energia di questo stato è quindi

$$\text{energy error} = \frac{\langle\psi|f(H)(H-E_0)f(H)|\psi\rangle}{\langle\psi|f(H)^2|\psi\rangle}$$ $$= \frac{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.$$

Per trasformare questo in un limite superiore più facile da comprendere, separiamo prima la somma al numeratore in termini con $E_k-E_0\le\delta$ e termini con $E_k-E_0>\delta$:

$$\text{energy error} = \frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} + \frac{\sum_{E_k> E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.$$

Possiamo limitare superiormente il primo termine con $\delta$,

$$\frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} < \frac{\delta\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} \le \delta,$$

dove il primo passaggio segue dal fatto che $E_k-E_0\le\delta$ per ogni $E_k$ nella somma, e il secondo passaggio segue dal fatto che la somma al numeratore è un sottoinsieme della somma al denominatore. Per il secondo termine, prima limitiamo inferiormente il denominatore con $|\gamma_0|^2$, poiché $f(E_0)^2=1$: aggiungendo tutto insieme, questo dà

$$\text{energy error} \le \delta + \frac{1}{|\gamma_0|^2}\sum_{E_k>E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0).$$

Per semplificare ciò che resta, notate che per tutti questi $E_k$, dalla definizione di $f$ sappiamo che $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Inoltre limitando superiormente $E_k-E_0<2\|H\|$ e limitando superiormente $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ si ottiene

$$\text{energy error} \le \delta + \frac{8}{|\gamma_0|^2}\|H\|\left(1 + \delta\right)^{-2d}.$$

Questo vale per qualsiasi $\delta>0$, quindi se impostiamo $\delta$ uguale al nostro errore obiettivo, allora il limite di errore sopra converge verso quello esponenzialmente con la dimensione di Krylov $2d=r$. Notate inoltre che se $\delta<E_1-E_0$ allora il termine $\delta$ scompare completamente nel limite sopra.

Per completare l'argomento, notiamo prima che quanto sopra è solo l'errore di energia del particolare stato $f(H)|\psi\rangle$, piuttosto che l'errore di energia dello stato a energia più bassa nello spazio di Krylov. Tuttavia, dal principio variazionale (Rayleigh-Ritz), l'errore di energia dello stato a energia più bassa nello spazio di Krylov è limitato superiormente dall'errore di energia di qualsiasi stato nello spazio di Krylov, quindi quanto sopra è anche un limite superiore sull'errore di energia dello stato a energia più bassa, cioè l'output dell'algoritmo di diagonalizzazione quantistica di Krylov.

Un'analisi simile a quella sopra può essere condotta tenendo conto anche del rumore e della procedura di sogliatura discussa nel notebook. Vedere [2] e [4] per questa analisi.

Appendice: dimostrazione dell'Affermazione 1

Quanto segue è per lo più derivato da [3], Teorema 3.1: sia $0 < a < b$ e sia $\Pi^*_d$ lo spazio dei polinomi residui (polinomi il cui valore in 0 è 1) di grado al massimo $d$. La soluzione a

$$\beta(a, b, d) = \min_{p \in \Pi^*_d} \max_{x \in [a, b]} |p(x)| \quad$$

è

$$p^*(x) = \frac{T_d\left(\frac{b + a - 2x}{b - a}\right)}{T_d\left(\frac{b + a}{b - a}\right)}, \quad$$

e il corrispondente valore minimo è

$$\beta(a, b, d) = T_d^{-1}\left(\frac{b + a}{b - a}\right).$$

Vogliamo convertire questo in una funzione che possa essere espressa naturalmente in termini di esponenziali complessi, perché queste sono le evoluzioni temporali reali che generano lo spazio di Krylov quantistico. Per farlo, è conveniente introdurre la seguente trasformazione di energie all'interno dell'intervallo spettrale dell'Hamiltoniana in numeri nell'intervallo $[0,1]$: definiamo

$$g(E) = \frac{1-\cos\big((E-E_0)dt\big)}{2},$$

dove $dt$ è un passo temporale tale che $-\pi < E_0dt < E_\text{max}dt < \pi$. Notate che $g(E_0)=0$ e $g(E)$ cresce man mano che $E$ si allontana da $E_0$.

Ora usando il polinomio $p^*(x)$ con i parametri a, b, d impostati come $a = g(E_0 + \delta)$, $b = 1$, e d = int(r/2), definiamo la funzione:

$$f(E) = p^* \left( g(E) \right) = \frac{T_d\left(1 + 2\frac{\cos\big((E-E_0)dt\big) - \cos\big(\delta\,dt\big)}{1 +\cos\big(\delta\,dt\big)}\right)}{T_d\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)}$$

dove $E_0$ è l'energia dello stato fondamentale. Possiamo vedere inserendo $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ che $f(E)$ è un polinomio trigonometrico di grado $d$, cioè una combinazione lineare di $e^{ijE\,dt}$ per $j=-d,-d+1,...,d-1,d$. Inoltre, dalla definizione di $p^*(x)$ sopra abbiamo che $f(E_0)=p(0)=1$ e per qualsiasi $E$ nell'intervallo spettrale tale che $\vert E-E_0 \vert > \delta$ abbiamo

$$|f(E)| \le \beta(a, b, d) = T_d^{-1}\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)$$ $$\leq 2\left(1 + \delta\right)^{-d} = 2\left(1 + \delta\right)^{-\lfloor k/2\rfloor}.$$

References

[1] N. Yoshioka, M. Amico, W. Kirby et al. "Diagonalization of large many-body Hamiltonians on a quantum processor". arXiv:2407.14431

[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] Å. Björck. "Numerical methods in matrix computations". Texts in Applied Mathematics. Springer International Publishing. (2014).

[4] William Kirby. "Analysis of quantum Krylov algorithms with errors". Quantum 8, 1457 (2024).

Sondaggio sul tutorial

Vi preghiamo di completare questo breve sondaggio per fornire un feedback su questo tutorial. Le vostre opinioni ci aiuteranno a migliorare i nostri contenuti e l'esperienza utente.

Link to survey