Krylovova kvantová diagonalizace mřížkových Hamiltonů

Odhadovaný čas: 20 minut na Heron r2 (POZNÁMKA: Jedná se pouze o odhad. Skutečný čas se může lišit.)

# 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

Pozadí

Tento tutoriál ukazuje, jak implementovat algoritmus Krylovovy kvantové diagonalizace (KQD) v rámci vzorů Qiskitu. Nejprve se dozvíš o teorii za algoritmem a pak uvidíš ukázku jeho spuštění na QPU.

Napříč obory nás zajímají vlastnosti základního stavu kvantových systémů. Mezi příklady patří pochopení základní povahy částic a sil, předpovídání a porozumění chování složitých materiálů a pochopení bio-chemických interakcí a reakcí. Kvůli exponenciálnímu růstu Hilbertova prostoru a korelacím, které vznikají v provázaných systémech, mají klasické algoritmy problém toto řešit pro kvantové systémy rostoucí velikosti. Na jednom konci spektra jsou stávající přístupy využívající kvantový hardware zaměřené na variační kvantové metody (například variační kvantový eigensolver). Tyto techniky čelí výzvám na současných zařízeních kvůli vysokému počtu volání funkcí potřebných v optimalizačním procesu, což přidává velkou režii na zdroje při použití pokročilých technik potlačení chyb, a tak omezuje jejich účinnost na malé systémy. Na druhém konci spektra jsou metody odolné vůči chybám s výkonnostními zárukami (například kvantová fázová estimace), které vyžadují hluboké Circuit, spustitelné pouze na zařízení odolném vůči chybám. Z těchto důvodů zde představujeme kvantový algoritmus založený na metodách podprostoru (jak je popsáno v tomto přehledovém článku), algoritmus Krylovovy kvantové diagonalizace (KQD). Tento algoritmus funguje dobře ve velkém měřítku [1] na stávajícím kvantovém hardwaru, sdílí podobné výkonnostní záruky jako fázová estimace, je kompatibilní s pokročilými technikami potlačení chyb a mohl by poskytovat výsledky klasicky nedosažitelné.

Požadavky

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalováno následující:

• Qiskit SDK v2.0 nebo novější, s podporou vizualizace
• Qiskit Runtime v0.22 nebo novější ( pip install qiskit-ibm-runtime )

Nastavení

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

Krok 1: Mapování klasických vstupů na kvantový problém

Krylovův prostor

Krylovův prostor $\mathcal{K}^r$ řádu $r$ je prostor rozepnutý vektory získanými násobením vyšších mocnin matice $A$, až do $r-1$, s 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\}$$

Pokud je matice $A$ Hamiltonián $H$, budeme odpovídající prostor nazývat mocninový Krylovův prostor $\mathcal{K}_P$. V případě, kdy $A$ je operátor časového vývoje generovaný Hamiltonánem $U=e^{-iHt}$, budeme prostor nazývat unitární Krylovův prostor $\mathcal{K}_U$. Mocninový Krylovův podprostor, který klasicky používáme, nelze přímo generovat na kvantovém počítači, protože $H$ není unitární operátor. Místo toho můžeme použít operátor časového vývoje $U = e^{-iHt}$, u nějž lze ukázat, že poskytuje podobné záruky konvergence jako mocninová metoda. Mocniny $U$ se pak stávají různými časovými kroky $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\}$$

Viz Přílohu pro podrobné odvození toho, jak unitární Krylovův prostor umožňuje přesně reprezentovat nízkoenergetické vlastní stavy.

Algoritmus Krylovovy kvantové diagonalizace

Pro daný Hamiltonián $H$, který chceme diagonalizovat, nejprve uvažujeme odpovídající unitární Krylovův prostor $\mathcal{K}_U$. Cílem je najít kompaktní reprezentaci Hamiltoniánu v $\mathcal{K}_U$, kterou budeme označovat $\tilde{H}$. Maticové prvky $\tilde{H}$, projekce Hamiltoniánu do Krylovova prostoru, lze vypočítat výpočtem následujících střední hodnot

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

Kde $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ jsou vektory unitárního Krylovova prostoru a $t_n = n dt$ jsou násobky zvoleného časového kroku $dt$. Na kvantovém počítači lze výpočet každého maticového prvku provést libovolným algoritmem, který umožňuje získat překryv mezi kvantovými stavy. Tento tutoriál se zaměřuje na Hadamardův test. Vzhledem k tomu, že $\mathcal{K}_U$ má dimenzi $r$, bude mít Hamiltonián promítnutý do podprostoru rozměry $r \times r$. Při dostatečně malém $r$ (obecně $r<<100$ stačí k dosažení konvergence odhadů vlastních energií) pak snadno diagonalizujeme promítnutý Hamiltonián $\tilde{H}$. Nemůžeme však přímo diagonalizovat $\tilde{H}$ kvůli neortogonalitě vektorů Krylovova prostoru. Musíme změřit jejich překryvy a sestavit matici $\tilde{S}$

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

To nám umožňuje řešit problém vlastních hodnot v neortogonálním prostoru (také nazývaný zobecněný problém vlastních hodnot)

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

Lze pak získat odhady vlastních hodnot a vlastních stavů $H$ pohledem na vlastní hodnoty $\tilde{H}$. Například odhad energie základního stavu se získá jako nejmenší vlastní hodnota $c$ a základní stav z odpovídajícího vlastního vektoru $\vec{c}$. Koeficienty v $\vec{c}$ určují příspěvek různých vektorů, které rozpínají $\mathcal{K}_U$.

Obrázek ukazuje reprezentaci Circuit modifikovaného Hadamardova testu, metody používané k výpočtu překryvu mezi různými kvantovými stavy. Pro každý maticový prvek $\tilde{H}_{i,j}$ se provede Hadamardův test mezi stavy $\vert \psi_i \rangle$, $\vert \psi_j \rangle$. To je zvýrazněno v obrázku barevným schématem pro maticové prvky a odpovídající operace $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$. Tedy pro výpočet všech maticových prvků promítnutého Hamiltoniánu $\tilde{H}$ je potřeba sada Hadamardových testů pro všechny možné kombinace vektorů Krylovova prostoru. Horní drát v obvodu Hadamardova testu je ancilla Qubit, který se měří buď v základně X nebo Y, jeho střední hodnota určuje hodnotu překryvu mezi stavy. Dolní drát představuje všechny Qubit systémového Hamiltoniánu. Operace $\text{Prep} \; \psi_i$ připraví systémové Qubit ve stavu $\vert \psi_i \rangle$ řízeném stavem ancilla Qubit (podobně pro $\text{Prep} \; \psi_j$) a operace $P$ představuje Pauliho rozklad systémového Hamiltoniánu $H = \sum_i P_i$. Podrobnější odvození operací počítaných Hadamardovým testem je uvedeno níže.

Definice Hamiltoniánu

Uvažujme Heisenbergův Hamiltonián pro $N$ Qubit na lineárním řetězci: $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)]

Nastavení parametrů algoritmu

Heuristicky volíme hodnotu časového kroku dt (na základě horních odhadů normy Hamiltoniánu). Ref [2] ukázal, že dostatečně malý časový krok je $\pi/\vert \vert H \vert \vert$, a že je výhodnější tuto hodnotu spíše podhodnotit než přehodnotit, protože přehodnocení může umožnit příspěvkům z vysokoenergetických stavů zkreslit i optimální stav v Krylovově prostoru. Na druhé straně příliš malá hodnota $dt$ vede k horšímu podmínění Krylovova podprostoru, protože vektory Krylovovy základny se mezi časovými kroky méně liší.

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

A nastavíme další parametry algoritmu. Pro účely tohoto tutoriálu se omezíme na použití Krylovova prostoru s pouhými pěti dimenzemi, což je docela omezující.

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

Příprava stavu

Zvolíme referenční stav $\vert \psi \rangle$, který má nějaký překryv se základním stavem. Pro tento Hamiltonián použijeme jako referenční stav stav s excitací ve středním Qubit $\vert 00..010...00 \rangle$.

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

Časový vývoj

Operátor časového vývoje generovaný daným Hamiltonánem: $U=e^{-iHt}$ lze realizovat pomocí Lie-Trotterovy aproximace.

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>

Hadamardův test

$$\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 |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{equation*}$$

Kde $P$ je jeden z členů v rozkladu Hamiltoniánu $H=\sum P$ a $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ jsou řízené operace, které připravují vektory $|\psi_i\rangle$, $|\psi_j\rangle$ unitárního Krylovova prostoru, přičemž $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Pro měření $X$ nejprve aplikuj $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*}$$

... pak změř:

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

Z identity $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Podobně měření $Y$ dává

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

Circuit Hadamardova testu může být hluboký Circuit po rozkladu na nativní Gate (který se ještě více prohloubí, pokud vezmeme v úvahu topologii zařízení)

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

Krok 2: Optimalizace problému pro spouštění na kvantovém hardwaru

Efektivní Hadamardův test

Hluboké obvody pro Hadamardův test, které jsme získali, můžeme optimalizovat zavedením určitých aproximací a využitím předpokladů o modelovém Hamiltoniánu. Uvažuj například následující obvod pro Hadamardův test:

Předpokládejme, že dokážeme klasicky vypočítat $E_0$, vlastní hodnotu $|0\rangle^N$ pod Hamiltoniánem $H$. To platí, pokud Hamiltonián zachovává symetrii U(1). Ačkoli se to může zdát silným předpokladem, existuje mnoho případů, kdy je bezpečné předpokládat, že existuje vakuový stav (v tomto případě se zobrazí na stav $|0\rangle^N$), který není ovlivněn působením Hamiltoniánu. To platí například pro chemické Hamiltoniány popisující stabilní molekuly (kde je zachováno množství elektronů). Pokud hradlo $\text{Prep} \; \psi$ připravuje požadovaný referenční stav $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, například pro přípravu stavu HF v chemii je $\text{Prep} \; \psi$ součinem jednoqubitových NOT, takže řízené $\text{Prep} \; \psi$ je jen součinem CNOT. Obvod výše pak implementuje následující stav před měřením:

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

kde jsme ve třetím řádku použili klasicky simulovatelný fázový posun $U\ket{0}^N = e^{i\phi}\ket{0}^N$. Střední hodnoty jsou proto získány jako

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

Díky těmto předpokladům jsme schopni vyjádřit střední hodnoty sledovaných operátorů s menším počtem řízených operací. Ve skutečnosti potřebujeme implementovat pouze řízenou přípravu stavu $\text{Prep} \; \psi$, nikoli řízené časové evoluce. Přeformulování výpočtu tímto způsobem nám umožní výrazně snížit hloubku výsledných obvodů.

Rozložení operátoru časové evoluce Trotterovou dekompozicí

Místo přesné implementace operátoru časové evoluce můžeme použít Trotterovu dekompozici pro implementaci jeho aproximace. Opakováním Trotterovy dekompozice určitého řádu několikrát dosáhneme dalšího snížení chyby způsobené aproximací. V následujícím kódu přímo sestavíme Trotterovu implementaci nejefektivnějším způsobem pro interakční graf uvažovaného Hamiltoniánu (pouze interakce nejbližších sousedů). V praxi vkládáme Pauliho rotace $R_{xx}$, $R_{yy}$, $R_{zz}$ s parametrizovaným úhlem $t$, které odpovídají přibližné implementaci $e^{-i (XX + YY + ZZ) t}$. Kvůli rozdílu v definici Pauliho rotací a časové evoluce, kterou se snažíme implementovat, musíme použít parametr $2*dt$, abychom dosáhli časové evoluce $dt$. Navíc pro lichý počet opakování Trotterových kroků obrátíme pořadí operací, což je funkčně ekvivalentní, ale umožňuje syntézovat sousední operace do jediného unitáru $SU(2)$. Tím získáme mnohem mělčí obvod než ten, který by byl získán pomocí obecné funkce 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)

Použití optimalizovaného obvodu pro přípravu stavu

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)

Šablonové obvody pro výpočet maticových prvků $\tilde{S}$ a $\tilde{H}$ pomocí Hadamardova testu

Jediný rozdíl mezi obvody použitými v Hadamardově testu bude fáze v operátoru časové evoluce a měřené observably. Proto můžeme připravit šablonový obvod reprezentující obecný obvod pro Hadamardův test se zástupnými symboly pro hradla závislá na operátoru časové evoluce.

# 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

Kombinací Trotterovy aproximace a neřízených unitárů jsme výrazně snížili hloubku Hadamardova testu.

Krok 3: Spuštění pomocí Qiskit primitiv

Inicializuj Backend a nastav parametry běhového prostředí

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)

Transpilace na QPU

Nejprve vybereme podmnožiny coupling map s „dobře" fungujícími qubity (kde „dobře" je zde poměrně libovolné — chceme hlavně vyloučit velmi špatně fungující qubity) a vytvoříme nový target pro Transpiler.

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

Pak transpilujeme virtuální obvod na nejlepší fyzické rozmístění v tomto novém targetu.

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]

Vytvoření PUBů pro spuštění s Estimatorem

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

Spuštění obvodů

Obvody pro $t=0$ lze vypočítat klasicky.

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)

Spuštění obvodů pro $S$ a $\tilde{H}$ s Estimatorem

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

Krok 4: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu

results = job.result()[0]

Výpočet efektivní matice Hamiltoniánu a překryvové matice

Nejprve vypočítej fázi, kterou stav $\vert 0 \rangle$ nashromáždí při neřízeném časovém vývoji

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

Jakmile máme výsledky spuštění Circuit, můžeme data následně zpracovat a vypočítat maticové prvky $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]$$

A maticové prvky $\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]$$

Nakonec můžeme vyřešit zobecněný problém vlastních čísel pro $\tilde{H}$:

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

a získat odhad energie základního stavu $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

Pro sektor jedné částice můžeme efektivně vypočítat základní stav tohoto sektoru Hamiltoniánu klasicky

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

Příloha: Krylovův podprostor z evolucí v reálném čase

Unitární Krylovův prostor je definován jako

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

pro nějaký časový krok $dt$, který určíme níže. Dočasně předpokládejme, že $r$ je sudé: definujme $d=r/2$. Všimni si, že když do Krylovova prostoru projektujeme Hamiltonián, je to nerozeznatelné od Krylovova prostoru

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

tedy kde jsou všechny časové evoluce posunuty zpět o $d$ časových kroků. Důvodem nerozeznatelnosti je to, že maticové prvky

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

jsou invariantní vůči celkovým posunům evolučního času, protože časové evoluce komutují s Hamiltoniánem. Pro liché $r$ lze použít analýzu pro $r-1$.

Chceme ukázat, že někde v tomto Krylovově prostoru je zaručeně stav s nízkou energií. Učiníme tak pomocí následujícího výsledku odvozeného z Teorému 3.1 v [3]:

Tvrzení 1: existuje funkce $f$ taková, že pro energie $E$ ve spektrálním rozsahu Hamiltoniánu (tedy mezi energií základního stavu a maximální energií)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ pro všechny hodnoty $E$, které leží $\ge\delta$ od $E_0$, tedy je exponenciálně potlačena
3. $f(E)$ je lineární kombinací $e^{ijE\,dt}$ pro $j=-d,-d+1,...,d-1,d$

Níže uvádíme důkaz, který lze bezpečně přeskočit, pokud člověk nechce porozumět celému rigoróznímu argumentu. Teď se zaměříme na důsledky výše uvedeného tvrzení. Z vlastnosti 3 výše vidíme, že výše uvedený posunutý Krylovův prostor obsahuje stav $f(H)|\psi\rangle$. To je náš stav s nízkou energií. Abychom pochopili proč, zapišme $|\psi\rangle$ v energetické vlastní bázi:

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

kde $|E_k\rangle$ je k-tý energetický vlastní stav a $\gamma_k$ je jeho amplituda v počátečním stavu $|\psi\rangle$. Vyjádřen pomocí toho, $f(H)|\psi\rangle$ je dán

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

přičemž využijeme skutečnosti, že $H$ lze nahradit $E_k$, když působí na vlastní stav $|E_k\rangle$. Energetická chyba tohoto stavu je tedy

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

Abychom toto převedli na horní ohraničení, které je snadnější pochopit, nejprve rozdělíme součet v čitateli na členy s $E_k-E_0\le\delta$ a členy s $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}.$$

První člen můžeme shora ohraničit hodnotou $\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,$$

kde první krok plyne z toho, že $E_k-E_0\le\delta$ pro každé $E_k$ v součtu, a druhý krok plyne z toho, že součet v čitateli je podmnožinou součtu ve jmenovateli. Pro druhý člen nejprve ohraničíme jmenovatel zdola hodnotou $|\gamma_0|^2$, protože $f(E_0)^2=1$: celkem dohromady to dává

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

Abychom zjednodušili zbývající část, všimni si, že pro všechna tato $E_k$ víme z definice $f$, že $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Po dalším ohraničení $E_k-E_0<2\|H\|$ shora a ohraničení $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ dostaneme

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

Toto platí pro libovolné $\delta>0$, takže pokud nastavíme $\delta$ rovno naší cílové chybě, pak výše uvedené ohraničení chyby konverguje k ní exponenciálně s Krylovovou dimenzí $2d=r$. Všimni si také, že pokud $\delta<E_1-E_0$, pak člen $\delta$ ve výše uvedeném ohraničení zcela zmizí.

K dokončení argumentu nejprve poznamenáme, že výše uvedené je pouze energetická chyba konkrétního stavu $f(H)|\psi\rangle$, nikoli energetická chyba nejnižšího energetického stavu v Krylovově prostoru. Nicméně z (Rayleigh-Ritzova) variačního principu je energetická chyba nejnižšího energetického stavu v Krylovově prostoru ohraničena shora energetickou chybou libovolného stavu v Krylovově prostoru, takže výše uvedené je také horní ohraničení energetické chyby nejnižšího energetického stavu, tedy výstupu algoritmu Krylovovy kvantové diagonalizace.

Podobnou analýzu jako výše lze provést, přičemž dodatečně zohledníme šum a postup prahování popsaný v notebooku. Tuto analýzu viz [2] a [4].

Příloha: důkaz Tvrzení 1

Následující text je převážně odvozen z [3], Teorém 3.1: nechť $0 < a < b$ a nechť $\Pi^*_d$ je prostor reziduálních polynomů (polynomů, jejichž hodnota v 0 je 1) stupně nejvýše $d$. Řešení

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

je

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

a odpovídající minimální hodnota je

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

Chceme toto převést na funkci, která může být přirozeně vyjádřena pomocí komplexních exponenciál, protože to jsou evoluce v reálném čase, které generují kvantový Krylovův prostor. Za tím účelem je vhodné zavést následující transformaci energií v rámci spektrálního rozsahu Hamiltoniánu na čísla v rozsahu $[0,1]$: definujme

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

kde $dt$ je časový krok takový, že $-\pi < E_0dt < E_\text{max}dt < \pi$. Všimni si, že $g(E_0)=0$ a $g(E)$ roste s tím, jak se $E$ vzdaluje od $E_0$.

Nyní pomocí polynomu $p^*(x)$ s parametry a, b, d nastavenými na $a = g(E_0 + \delta)$, $b = 1$ a d = int(r/2) definujeme funkci:

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

kde $E_0$ je energie základního stavu. Dosazením $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ vidíme, že $f(E)$ je trigonometrický polynom stupně $d$, tedy lineární kombinace $e^{ijE\,dt}$ pro $j=-d,-d+1,...,d-1,d$. Dále z definice $p^*(x)$ výše plyne, že $f(E_0)=p(0)=1$ a pro libovolné $E$ ve spektrálním rozsahu takové, že $\vert E-E_0 \vert > \delta$, platí

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

Reference

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

Průzkum k tutoriálu

Vyplň prosím tento krátký průzkum a poskytni zpětnou vazbu k tomuto tutoriálu. Tvoje podněty nám pomohou zlepšit náš obsah a uživatelský zážitek.

Odkaz na průzkum

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.