Kwantowa diagonalizacja Krylova sieciowych hamiltonianów

Szacowany czas wykonania: 20 minut na procesorze Heron r2 (UWAGA: To jedynie szacunek. Rzeczywisty czas może się różnić.)

# 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

Podstawy teoretyczne

Ten samouczek pokazuje, jak zaimplementować algorytm kwantowej diagonalizacji Krylova (KQD) w kontekście wzorców Qiskit. Najpierw poznasz teorię stojącą za algorytmem, a następnie zobaczysz demonstrację jego wykonania na QPU.

W wielu dziedzinach nauki interesujemy się własnościami stanu podstawowego układów kwantowych. Przykłady obejmują rozumienie fundamentalnej natury cząstek i sił, przewidywanie i rozumienie zachowania złożonych materiałów oraz rozumienie biochemicznych interakcji i reakcji. Z uwagi na wykładniczy wzrost przestrzeni Hilberta i korelacje pojawiające się w splątanych układach, klasyczne algorytmy mają trudności z rozwiązaniem tego problemu dla kwantowych układów o rosnącym rozmiarze. Na jednym końcu spektrum są istniejące podejścia wykorzystujące sprzęt kwantowy, skupiające się na wariacyjnych metodach kwantowych (na przykład wariacjnym wyznaczaniu wartości własnych). Te techniki napotykają trudności z obecnymi urządzeniami ze względu na dużą liczbę wywołań funkcji wymaganych w procesie optymalizacji, co powoduje duże nakłady zasobów przy zastosowaniu zaawansowanych technik korekcji błędów, ograniczając ich skuteczność do małych układów. Na drugim końcu spektrum znajdują się tolerujące błędy metody kwantowe z gwarancjami wydajności (na przykład kwantowa estymacja fazy), które wymagają głębokich obwodów możliwych do wykonania jedynie na urządzeniu odpornym na błędy. Z tych powodów przedstawiamy tu algorytm kwantowy oparty na metodach podprzestrzeni (opisany w tym artykule przeglądowym) — algorytm kwantowej diagonalizacji Krylova (KQD). Algorytm ten dobrze skaluje się dla dużych układów [1] na istniejącym sprzęcie kwantowym, posiada podobne gwarancje wydajności jak estymacja fazy, jest kompatybilny z zaawansowanymi technikami łagodzenia błędów i może dostarczać wyniki niedostępne dla klasycznych komputerów.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane:

• Qiskit SDK v2.0 lub nowszy, z obsługą wizualizacji
• Qiskit Runtime v0.22 lub nowszy ( pip install qiskit-ibm-runtime )

Konfiguracja

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)

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

    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: Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Przestrzeń Krylova

Przestrzeń Krylova $\mathcal{K}^r$ rzędu $r$ to przestrzeń rozpiętą przez wektory uzyskane przez mnożenie wyższych potęg macierzy $A$, aż do $r-1$, przez wektor referencyjny $\vert v \rangle$.

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

Jeśli macierz $A$ jest hamiltonianem $H$, odpowiadającą przestrzeń nazwiemy potęgową przestrzenią Krylova $\mathcal{K}_P$. W przypadku gdy $A$ jest operatorem ewolucji czasowej generowanym przez hamiltonian $U=e^{-iHt}$, przestrzeń nazwiemy unitarną przestrzenią Krylova $\mathcal{K}_U$. Potęgowej podprzestrzeni Krylova używanej klasycznie nie można bezpośrednio wygenerować na komputerze kwantowym, ponieważ $H$ nie jest operatorem unitarnym. Zamiast tego możemy użyć operatora ewolucji czasowej $U = e^{-iHt}$, który — jak można wykazać — daje podobne gwarancje zbieżności jak metoda potęgowa. Potęgi $U$ stają się wówczas różnymi krokami czasowymi $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\}$$

Szczegółowe wyprowadzenie tego, jak unitarna przestrzeń Krylova pozwala dokładnie reprezentować niskenergetyczne stany własne, znajdziesz w Dodatku.

Algorytm kwantowej diagonalizacji Krylova

Mając hamiltonian $H$, który chcemy zdiagonalizować, najpierw rozważamy odpowiadającą mu unitarną przestrzeń Krylova $\mathcal{K}_U$. Celem jest znalezienie zwartej reprezentacji hamiltonianu w $\mathcal{K}_U$, którą oznaczymy $\tilde{H}$. Elementy macierzowe $\tilde{H}$, czyli rzutowanie hamiltonianu na przestrzeń Krylova, można obliczyć wyznaczając następujące wartości oczekiwane:

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

gdzie $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ są wektorami unitarnej przestrzeni Krylova, a $t_n = n dt$ są wielokrotnościami wybranego kroku czasowego $dt$. Na komputerze kwantowym obliczenie każdego elementu macierzowego można wykonać dowolnym algorytmem pozwalającym uzyskać iloczyn skalarny między stanami kwantowymi. Ten samouczek koncentruje się na teście Hadamarda. Przyjmując, że $\mathcal{K}_U$ ma wymiar $r$, hamiltonian rzutowany na podprzestrzeń będzie miał wymiary $r \times r$. Przy wystarczająco małym $r$ (zazwyczaj $r<<100$ wystarcza do uzyskania zbieżności szacowań energii własnych) możemy wówczas łatwo zdiagonalizować rzutowany hamiltonian $\tilde{H}$. Nie możemy jednak bezpośrednio zdiagonalizować $\tilde{H}$ ze względu na nieortogonalność wektorów przestrzeni Krylova. Musimy zmierzyć ich iloczyny skalarne i skonstruować macierz $\tilde{S}$:

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

Pozwala to rozwiązać zagadnienie własne w przestrzeni nieortogonalnej (zwane też uogólnionym zagadnieniem własnym):

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

Można następnie uzyskać szacunki wartości własnych i stanów własnych $H$ poprzez zbadanie wartości własnych $\tilde{H}$. Na przykład szacunek energii stanu podstawowego otrzymuje się, biorąc najmniejszą wartość własną $c$, a stan podstawowy z odpowiadającego wektora własnego $\vec{c}$. Współczynniki $\vec{c}$ wyznaczają udział różnych wektorów rozpinających $\mathcal{K}_U$.

Rysunek przedstawia reprezentację obwodową zmodyfikowanego testu Hadamarda — metody stosowanej do obliczania iloczynu skalarnego między różnymi stanami kwantowymi. Dla każdego elementu macierzowego $\tilde{H}_{i,j}$ wykonuje się test Hadamarda między stanami $\vert \psi_i \rangle$ i $\vert \psi_j \rangle$. Jest to zaznaczone na rysunku schematem kolorów elementów macierzy i odpowiadających im operacji $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$. Do obliczenia wszystkich elementów macierzowych rzutowanego hamiltonianu $\tilde{H}$ potrzebny jest zestaw testów Hadamarda dla wszystkich możliwych kombinacji wektorów przestrzeni Krylova. Górna linia w obwodzie testu Hadamarda to qubit pomocniczy (ancilla), mierzony w bazie X lub Y; jego wartość oczekiwana wyznacza wartość iloczynu skalarnego między stanami. Dolna linia reprezentuje wszystkie qubity hamiltonianu układu. Operacja $\text{Prep} \; \psi_i$ przygotowuje qubit układu w stanie $\vert \psi_i \rangle$ sterowanym stanem qubita pomocniczego (analogicznie $\text{Prep} \; \psi_j$), a operacja $P$ reprezentuje dekompozycję Pauliego hamiltonianu układu $H = \sum_i P_i$. Bardziej szczegółowe wyprowadzenie operacji obliczanych przez test Hadamarda podane jest poniżej.

Definicja hamiltonianu

Rozważmy hamiltonian Heisenberga dla $N$ qubitów na łańcuchu liniowym: $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)]

Ustawianie parametrów algorytmu

Heurystycznie dobieramy wartość kroku czasowego dt (na podstawie górnych ograniczeń normy hamiltonianu). Praca [2] wykazała, że wystarczająco małym krokiem czasowym jest $\pi/\vert \vert H \vert \vert$, oraz że do pewnego stopnia lepiej niedoszacować tę wartość niż ją przeszacować — przeszacowanie może bowiem dopuścić wkłady stanów wysokoenergetycznych, które zakłócają nawet optymalny stan w przestrzeni Krylova. Z drugiej strony zbyt małe $dt$ prowadzi do gorszego uwarunkowania podprzestrzeni Krylova, ponieważ wektory bazy Krylova różnią się od siebie mniej z kroku na krok.

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

Ustawiamy też pozostałe parametry algorytmu. Na potrzeby tego samouczka ograniczymy się do przestrzeni Krylova o zaledwie pięciu wymiarach, co jest dość restrykcyjne.

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

Przygotowanie stanu

Wybieramy stan referencyjny $\vert \psi \rangle$, który ma niezerowe nakładanie ze stanem podstawowym. Dla tego hamiltonianu używamy stanu z wzbudzeniem na środkowym qubicie $\vert 00..010...00 \rangle$ jako naszego stanu referencyjnego.

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

Ewolucja czasowa

Operator ewolucji czasowej generowany przez dany hamiltonian $U=e^{-iHt}$ możemy zrealizować za pomocą aproksymacji Lie-Trottera.

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 Hadamarda

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

gdzie $P$ jest jednym z wyrazów dekompozycji hamiltonianu $H=\sum P$, a $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ to sterowane operacje przygotowujące wektory $|\psi_i\rangle$, $|\psi_j\rangle$ unitarnej przestrzeni Krylova, przy czym $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Aby zmierzyć $X$, najpierw stosujemy $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*}$$

... a następnie mierzymy:

$$\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 tożsamości $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Analogicznie, pomiar $Y$ daje:

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

Obwód testu Hadamarda może być głębokim obwodem po dekompozycji do bramek natywnych (co jeszcze bardziej wzrośnie przy uwzględnieniu topologii urządzenia).

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: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Wydajny test Hadamarda

Możemy zoptymalizować głębokie obwody dla testu Hadamarda, które uzyskaliśmy, wprowadzając pewne przybliżenia i opierając się na założeniach dotyczących hamiltoniana modelu. Rozważmy na przykład następujący obwód dla testu Hadamarda:

Załóżmy, że możemy klasycznie obliczyć $E_0$, wartość własną $|0\rangle^N$ pod hamiltonianem $H$. Warunek ten jest spełniony, gdy hamiltonian zachowuje symetrię U(1). Choć może się to wydawać silnym założeniem, istnieje wiele przypadków, w których bezpiecznie można przyjąć, że istnieje stan próżni (w tym przypadku odpowiadający stanowi $|0\rangle^N$), na który działanie hamiltonianu nie ma wpływu. Jest tak na przykład w przypadku hamiltonianów chemicznych opisujących stabilną cząsteczkę (gdzie liczba elektronów jest zachowana). Zakładając, że bramka $\text{Prep} \; \psi$ przygotowuje pożądany stan referencyjny $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, na przykład w celu przygotowania stanu HF dla chemii $\text{Prep} \; \psi$ byłby iloczynem jednoQubitowych bramek NOT, a zatem kontrolowane-$\text{Prep} \; \psi$ jest jedynie iloczynem bramek CNOT. Wówczas powyższy obwód implementuje następujący stan przed pomiarem:

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

gdzie w trzecim wierszu użyliśmy klasycznie symulowanego przesunięcia fazy $U\ket{0}^N = e^{i\phi}\ket{0}^N$. W związku z tym wartości oczekiwane oblicza się 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}$$

Dzięki tym założeniom udało nam się zapisać wartości oczekiwane interesujących nas operatorów przy użyciu mniejszej liczby operacji kontrolowanych. W rzeczywistości musimy zaimplementować jedynie kontrolowane przygotowanie stanu $\text{Prep} \; \psi$, a nie kontrolowane ewolucje czasowe. Przeformułowanie naszego obliczenia w powyższy sposób pozwoli nam znacznie zmniejszyć głębokość wynikowych obwodów.

Dekompozycja operatora ewolucji czasowej przy użyciu dekompozycji Trottera

Zamiast implementować operator ewolucji czasowej w sposób dokładny, możemy użyć dekompozycji Trottera, aby uzyskać jego przybliżenie. Wielokrotne powtarzanie dekompozycji Trottera określonego rzędu pozwala nam dodatkowo zmniejszyć błąd wprowadzony przez to przybliżenie. W poniższym kodzie budujemy implementację Trottera w najbardziej efektywny sposób dla grafu interakcji rozpatrywanego hamiltoniana (tylko najbliższe interakcje sąsiedzkie). W praktyce wstawiamy rotacje Pauliego $R_{xx}$, $R_{yy}$, $R_{zz}$ z parametryzowanym kątem $t$, które odpowiadają przybliżonej implementacji $e^{-i (XX + YY + ZZ) t}$. Ze względu na różnicę w definicji rotacji Pauliego i ewolucji czasowej, którą chcemy zaimplementować, musimy użyć parametru $2*dt$, aby uzyskać ewolucję czasową o $dt$. Ponadto odwracamy kolejność operacji dla nieparzystej liczby kroków Trottera, co jest funkcjonalnie równoważne, ale pozwala na syntezę sąsiednich operacji w jedną unitarną $SU(2)$. Daje to znacznie płytszy obwód niż ten uzyskiwany za pomocą ogólnej funkcjonalności 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)

Użycie zoptymalizowanego obwodu do przygotowania stanu

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)

Obwody szablonowe do obliczania elementów macierzowych $\tilde{S}$ i $\tilde{H}$ za pomocą testu Hadamarda

Jedyną różnicą między obwodami używanymi w teście Hadamarda będzie faza w operatorze ewolucji czasowej oraz mierzone obserwable. Możemy więc przygotować obwód szablonowy reprezentujący ogólny obwód dla testu Hadamarda, z miejscami zastępczymi dla bramek zależących od operatora ewolucji czasowej.

# 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

Znacząco zmniejszyliśmy głębokość testu Hadamarda dzięki kombinacji przybliżenia Trottera i niekontrolowanych unitarnych.

Krok 3: Wykonanie przy użyciu prymitywów Qiskit

Utwórz instancję backendu i ustaw parametry środowiska uruchomieniowego

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)

Transpilacja do QPU

Najpierw wybierzmy podzbiory mapy sprzężeń z „dobrymi" kubitami (pojęcie „dobry" jest tutaj dość arbitralne — głównie chcemy unikać kubitów o bardzo słabych parametrach) i utwórzmy nowy cel dla Transpilera

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

Następnie dokonaj transpilacji układu wirtualnego do najlepszego układu fizycznego w tym nowym celu

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]

Tworzenie PUB-ów do wykonania z 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)

Uruchamianie układów

Układy dla $t=0$ można obliczyć klasycznie

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)

Wykonaj układy dla $S$ i $\tilde{H}$ za pomocą Estimatora

# 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: Przetwarzanie końcowe i zwracanie wyniku w żądanym formacie klasycznym

results = job.result()[0]

Obliczanie efektywnego Hamiltonianu i macierzy nakładania

Najpierw oblicz fazę nagromadzoną przez stan $\vert 0 \rangle$ podczas niekontrolowanej ewolucji czasowej

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

Gdy mamy wyniki wykonań układów, możemy przetworzyć dane w celu obliczenia elementów macierzy $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]$$

I elementy macierzy $\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]$$

Wreszcie możemy rozwiązać uogólniony problem wartości własnych dla $\tilde{H}$:

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

i uzyskać oszacowanie energii stanu podstawowego $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

Dla sektora jednocząstkowego możemy efektywnie obliczyć klasycznie stan podstawowy tego sektora Hamiltonianu

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

Dodatek: Podprzestrzeń Krylova z ewolucji w czasie rzeczywistym

Unitarna przestrzeń Krylova jest zdefiniowana jako

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

dla pewnego kroku czasowego $dt$, który wyznaczymy później. Załóżmy tymczasowo, że $r$ jest parzyste: wtedy zdefiniuj $d=r/2$. Zauważmy, że kiedy rzutujemy Hamiltonian na powyższą przestrzeń Krylova, jest ona nierozróżnialna od przestrzeni Krylova

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

to znaczy takiej, gdzie wszystkie ewolucje czasowe są przesunięte wstecz o $d$ kroków czasowych. Nierozróżnialność wynika z tego, że elementy macierzy

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

są niezmiennicze ze względu na globalne przesunięcia czasu ewolucji, ponieważ ewolucje czasowe komutują z Hamiltonianem. Dla nieparzystego $r$ można użyć analizy dla $r-1$.

Chcemy pokazać, że gdzieś w tej przestrzeni Krylova jest zagwarantowany stan niskoenergetyczny. Robimy to za pomocą następującego twierdzenia, które pochodzi z Twierdzenia 3.1 w [3]:

Twierdzenie 1: istnieje funkcja $f$ taka, że dla energii $E$ w zakresie spektralnym Hamiltonianu (to znaczy między energią stanu podstawowego a energią maksymalną)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ dla wszystkich wartości $E$ oddalonych o $\ge\delta$ od $E_0$, czyli jest wykładniczo tłumiona
3. $f(E)$ jest kombinacją liniową $e^{ijE\,dt}$ dla $j=-d,-d+1,...,d-1,d$

Poniżej podajemy dowód, który można bezpiecznie pominąć, o ile nie zależy nam na pełnym, rygorystycznym uzasadnieniu. Skupmy się teraz na konsekwencjach powyższego twierdzenia. Z własności 3 powyżej wynika, że przesunięta przestrzeń Krylova zawiera stan $f(H)|\psi\rangle$. To jest nasz stan niskoenergetyczny. Aby to zrozumieć, zapiszmy $|\psi\rangle$ w bazie stanów własnych energii:

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

gdzie $|E_k\rangle$ to k-ty stan własny energii, a $\gamma_k$ to jego amplituda w stanie początkowym $|\psi\rangle$. W tej reprezentacji $f(H)|\psi\rangle$ jest dane wzorem

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

korzystając z faktu, że $H$ można zastąpić przez $E_k$, gdy działa na stan własny $|E_k\rangle$. Błąd energetyczny tego stanu wynosi zatem

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

Aby zamienić to na łatwiejsze do interpretacji ograniczenie górne, najpierw rozdzielamy sumę w liczniku na wyrazy z $E_k-E_0\le\delta$ i wyrazy z $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}.$$

Możemy ograniczyć pierwszy wyraz z góry przez $\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,$$

gdzie pierwszy krok wynika z $E_k-E_0\le\delta$ dla każdego $E_k$ w sumie, a drugi krok wynika z tego, że suma w liczniku jest podzbiorem sumy w mianowniku. Dla drugiego wyrazu najpierw ograniczamy mianownik z dołu przez $|\gamma_0|^2$, ponieważ $f(E_0)^2=1$: łącząc wszystko razem, otrzymujemy

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

Aby uprościć to, co pozostało, zauważmy, że dla wszystkich tych $E_k$, z definicji $f$ wiemy, że $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Dodatkowo ograniczając $E_k-E_0<2\|H\|$ z góry i ograniczając $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ z góry, otrzymujemy

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

Zachodzi to dla dowolnego $\delta>0$, więc jeśli ustawimy $\delta$ równe naszemu docelowemu błędowi, powyższe ograniczenie błędu dąży do niego wykładniczo wraz z wymiarem przestrzeni Krylova $2d=r$. Należy też zauważyć, że jeśli $\delta<E_1-E_0$, to wyraz $\delta$ całkowicie znika w powyższym ograniczeniu.

Aby zakończyć rozumowanie, najpierw zauważmy, że powyższe dotyczy tylko błędu energetycznego konkretnego stanu $f(H)|\psi\rangle$, a nie błędu energetycznego stanu o najniższej energii w przestrzeni Krylova. Jednak na mocy zasady wariacyjnej (Rayleigha-Ritza) błąd energetyczny stanu o najniższej energii w przestrzeni Krylova jest ograniczony z góry przez błąd energetyczny dowolnego stanu w tej przestrzeni, więc powyższe jest również ograniczeniem górnym błędu energetycznego stanu o najniższej energii, czyli wyjścia algorytmu kwantowej diagonalizacji Krylova.

Podobna analiza jak powyższa może być przeprowadzona z uwzględnieniem szumu i procedury progowania omówionej w notatniku. Patrz [2] i [4] w celu zapoznania się z tą analizą.

Dodatek: dowód Twierdzenia 1

Poniższe wywodzi się głównie z [3], Twierdzenia 3.1: niech $0 < a < b$ i niech $\Pi^*_d$ będzie przestrzenią wielomianów resztkowych (wielomianów, których wartość w 0 wynosi 1) stopnia co najwyżej $d$. Rozwiązanie

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

wynosi

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

a odpowiadająca minimalna wartość to

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

Chcemy przekształcić to w funkcję, która może być wyrażona naturalnie w kategoriach eksponencjałów zespolonych, ponieważ to właśnie rzeczywiste ewolucje czasowe generują kwantową przestrzeń Krylova. W tym celu wygodnie jest wprowadzić następującą transformację energii w zakresie spektralnym Hamiltonianu na liczby z przedziału $[0,1]$: zdefiniujmy

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

gdzie $dt$ jest krokiem czasowym takim, że $-\pi < E_0dt < E_\text{max}dt < \pi$. Zauważmy, że $g(E_0)=0$ i $g(E)$ rośnie, gdy $E$ oddala się od $E_0$.

Teraz używając wielomianu $p^*(x)$ z parametrami a, b, d ustawionymi na $a = g(E_0 + \delta)$, $b = 1$, oraz d = int(r/2), definiujemy funkcję:

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

gdzie $E_0$ to energia stanu podstawowego. Wstawiając $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$, można zobaczyć, że $f(E)$ jest wielomianem trygonometrycznym stopnia $d$, czyli kombinacją liniową $e^{ijE\,dt}$ dla $j=-d,-d+1,...,d-1,d$. Ponadto z definicji $p^*(x)$ powyżej wynika, że $f(E_0)=p(0)=1$ i dla dowolnego $E$ w zakresie spektralnym takiego, że $\vert E-E_0 \vert > \delta$, mamy

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

Literatura

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

Ankieta dotycząca samouczka

Wypełnij tę krótką ankietę, aby przekazać nam opinię o tym samouczku. Twoje uwagi pomogą nam udoskonalić nasze materiały i poprawić doświadczenie użytkowników.

Link do ankiety

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.