Крилівська квантова діагоналізація решіткових гамільтоніанів

Оцінка використання: 20 хвилин на Heron r2 (ПРИМІТКА: Це лише оцінка. Фактичний час виконання може відрізнятися.)

# 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

Передумови

Цей посібник демонструє, як реалізувати алгоритм крилівської квантової діагоналізації (KQD) в контексті патернів Qiskit. Ви спочатку дізнаєтесь про теорію, що лежить в основі алгоритму, а потім побачите демонстрацію його виконання на QPU.

У різних дисциплінах ми зацікавлені у вивченні властивостей основного стану квантових систем. Приклади включають розуміння фундаментальної природи частинок та сил, прогнозування та розуміння поведінки складних матеріалів і розуміння біохімічних взаємодій та реакцій. Через експоненційне зростання простору Гільберта та кореляції, що виникають у заплутаних системах, класичні алгоритми мають труднощі з розв'язанням цієї проблеми для квантових систем зростаючого розміру. На одному кінці спектра знаходиться існуючий підхід, який використовує переваги квантового обладнання, зосереджуючись на варіаційних квантових методах (наприклад, варіаційний квантовий власний розв'язувач). Ці методи стикаються з викликами на сучасних пристроях через високу кількість викликів функцій, необхідних у процесі оптимізації, що додає значні накладні витрати на ресурси після впровадження передових методів пом'якшення помилок, таким чином обмежуючи їх ефективність малими системами. На іншому кінці спектра знаходяться відмовостійкі квантові методи з гарантіями продуктивності (наприклад, оцінка квантової фази), які вимагають глибоких схем, що можуть бути виконані лише на відмовостійкому пристрої. З цих причин ми представляємо тут квантовий алгоритм, заснований на методах підпросторів (як описано в цій оглядовій статті), алгоритм крилівської квантової діагоналізації (KQD). Цей алгоритм добре працює у великому масштабі [1] на існуючому квантовому обладнанні, має подібні гарантії продуктивності як оцінка фази, сумісний з передовими методами пом'якшення помилок і може надати результати, які є класично недоступними.

Вимоги

Перед початком цього посібника переконайтеся, що у Вас встановлено наступне:

• Qiskit SDK v2.0 або новіше, з підтримкою візуалізації
• Qiskit Runtime v0.22 або новіше ( pip install qiskit-ibm-runtime )

Налаштування

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

Крок 1: Відображення класичних вхідних даних на квантову задачу

Простір Крилова

Простір Крилова $\mathcal{K}^r$ порядку $r$ - це простір, охоплений векторами, отриманими шляхом множення вищих степенів матриці $A$, до $r-1$, на опорний вектор $\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\}$$

Якщо матриця $A$ є гамільтоніаном $H$, ми будемо посилатися на відповідний простір як на степеневий простір Крилова $\mathcal{K}_P$. У випадку, коли $A$ є оператором часової еволюції, згенерованим гамільтоніаном $U=e^{-iHt}$, ми будемо посилатися на простір як на унітарний простір Крилова $\mathcal{K}_U$. Степеневий підпростір Крилова, який ми використовуємо класично, не може бути згенерований безпосередньо на квантовому комп'ютері, оскільки $H$ не є унітарним оператором. Натомість ми можемо використовувати оператор часової еволюції $U = e^{-iHt}$, який, як можна показати, дає подібні гарантії збіжності як степеневий метод. Степені $U$ потім стають різними часовими кроками $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\}$$

Дивіться Додаток для детального виведення того, як унітарний простір Крилова дозволяє точно представляти власні стани з низькою енергією.

Алгоритм крилівської квантової діагоналізації

Маючи гамільтоніан $H$, який ми хочемо діагоналізувати, спочатку розглянемо відповідний унітарний простір Крилова $\mathcal{K}_U$. Мета полягає в тому, щоб знайти компактне представлення гамільтоніана в $\mathcal{K}_U$, яке ми будемо називати $\tilde{H}$. Матричні елементи $\tilde{H}$, проекція гамільтоніана в просторі Крилова, можуть бути обчислені шляхом обчислення наступних очікуваних значень

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

Де $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ є векторами унітарного простору Крилова, а $t_n = n dt$ є кратними обраного часового кроку $dt$. На квантовому комп'ютері обчислення кожного матричного елемента може бути виконано за допомогою будь-якого алгоритму, який дозволяє отримати перекриття між квантовими станами. Цей посібник зосереджується на тесті Адамара. Враховуючи, що $\mathcal{K}_U$ має розмірність $r$, гамільтоніан, спроектований в підпростір, матиме розміри $r \times r$. При достатньо малому $r$ (зазвичай $r<<100$ достатньо для отримання збіжності оцінок власних енергій) ми можемо легко діагоналізувати спроектований гамільтоніан $\tilde{H}$. Однак ми не можемо безпосередньо діагоналізувати $\tilde{H}$ через неортогональність векторів простору Крилова. Нам доведеться виміряти їх перекриття і побудувати матрицю $\tilde{S}$

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

Це дозволяє нам розв'язати задачу на власні значення в неортогональному просторі (також називається узагальненою задачею на власні значення)

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

Потім можна отримати оцінки власних значень і власних станів $H$, розглядаючи їх для $\tilde{H}$. Наприклад, оцінка енергії основного стану отримується шляхом взяття найменшого власного значення $c$, а основний стан - з відповідного власного вектора $\vec{c}$. Коефіцієнти в $\vec{c}$ визначають внесок різних векторів, що охоплюють $\mathcal{K}_U$.

На рисунку показано схемне представлення модифікованого тесту Адамара, методу, який використовується для обчислення перекриття між різними квантовими станами. Для кожного матричного елемента $\tilde{H}_{i,j}$ виконується тест Адамара між станом $\vert \psi_i \rangle$, $\vert \psi_j \rangle$. Це виділено на рисунку кольоровою схемою для матричних елементів і відповідних операцій $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$. Таким чином, для обчислення всіх матричних елементів спроектованого гамільтоніана $\tilde{H}$ потрібен набір тестів Адамара для всіх можливих комбінацій векторів простору Крилова. Верхній провід у схемі тесту Адамара - це допоміжний кубіт, який вимірюється в базисі X або Y, його очікуване значення визначає значення перекриття між станами. Нижній провід представляє всі кубіти системного гамільтоніана. Операція $\text{Prep} \; \psi_i$ готує системний кубіт у стані $\vert \psi_i \rangle$, керований станом допоміжного кубіта (аналогічно для $\text{Prep} \; \psi_j$), а операція $P$ представляє розкладання Паулі системного гамільтоніана $H = \sum_i P_i$. Більш детальне виведення операцій, обчислених тестом Адамара, наведено нижче.

Визначення гамільтоніана

Розглянемо гамільтоніан Гайзенберга для $N$ кубітів на лінійному ланцюжку: $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)]

Встановлення параметрів алгоритму

Ми евристично обираємо значення для часового кроку dt (на основі верхніх меж на нормі гамільтоніана). Посилання [2] показало, що достатньо малий часовий крок становить $\pi/\vert \vert H \vert \vert$, і що краще до певної міри недооцінити це значення, ніж переоцінити, оскільки переоцінка може дозволити внескам від високоенергетичних станів зіпсувати навіть оптимальний стан у просторі Крилова. З іншого боку, вибір занадто малого $dt$ призводить до гіршої обумовленості підпростору Крилова, оскільки базисні вектори Крилова менше відрізняються від кроку до кроку.

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

І встановіть інші параметри алгоритму. Для цілей цього посібника ми обмежимося використанням простору Крилова лише з п'ятьма вимірами, що є досить обмежуючим.

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

Підготовка стану

Виберіть опорний стан $\vert \psi \rangle$, який має деяке перекриття з основним станом. Для цього гамільтоніана ми використовуємо стан із збудженням у середньому кубіті $\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)

Часова еволюція

Ми можемо реалізувати оператор часової еволюції, згенерований даним гамільтоніаном: $U=e^{-iHt}$ через наближення Лі-Троттера.

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>

Тест Адамара

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

Де $P$ є одним із членів у розкладанні гамільтоніана $H=\sum P$, а $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ є керованими операціями, які готують $|\psi_i\rangle$, $|\psi_j\rangle$ вектори унітарного простору Крилова, з $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Щоб виміряти $X$, спочатку застосуємо $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*}$$

... потім вимірюємо:

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

З тотожності $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Аналогічно, вимірювання $Y$ дає

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

Схема тесту Адамара може бути глибокою схемою після того, як ми розкладемо її на нативні вентилі (що збільшиться ще більше, якщо ми врахуємо топологію пристрою)

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

Крок 2: Оптимізація задачі для виконання на квантовому обладнанні

Ефективний тест Адамара

Ми можемо оптимізувати глибокі схеми для тесту Адамара, які ми отримали, вводячи деякі наближення та спираючись на деякі припущення щодо модельного гамільтоніана. Наприклад, розглянемо наступну схему для тесту Адамара:

Припустимо, що ми можемо класично обчислити $E_0$, власне значення $|0\rangle^N$ під дією гамільтоніана $H$. Це виконується, коли гамільтоніан зберігає U(1) симетрію. Хоча це може здаватися сильним припущенням, існує багато випадків, коли безпечно припустити, що існує вакуумний стан (у цьому випадку він відображається на стан $|0\rangle^N$), який не змінюється під дією гамільтоніана. Це справедливо, наприклад, для гамільтоніанів хімії, які описують стабільну молекулу (де кількість електронів зберігається). Враховуючи, що вентиль $\text{Prep} \; \psi$ готує бажаний опорний стан $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, наприклад, для підготовки HF стану для хімії $\text{Prep} \; \psi$ буде добутком однокубітних NOT, тому керований-$\text{Prep} \; \psi$ є просто добутком CNOT. Тоді наведена вище схема реалізує наступний стан перед вимірюванням:

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

де ми використали класично обчислюваний зсув фази $U\ket{0}^N = e^{i\phi}\ket{0}^N$ у третьому рядку. Тому очікувані значення отримуються як

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

Використовуючи ці припущення, ми змогли записати очікувані значення операторів, що становлять інтерес, з меншою кількістю керованих операцій. Фактично, нам потрібно реалізувати лише керовану підготовку стану $\text{Prep} \; \psi$, а не керовану часову еволюцію. Переформулювання нашого обчислення, як показано вище, дозволить нам значно зменшити глибину отриманих схем.

Розклад оператора часової еволюції за допомогою розкладу Троттера

Замість точної реалізації оператора часової еволюції ми можемо використовувати розклад Троттера для реалізації його наближення. Повторення кілька разів розкладу Троттера певного порядку дає нам подальше зменшення похибки, введеної наближенням. У наступному ми безпосередньо будуємо реалізацію Троттера найефективнішим способом для графа взаємодії гамільтоніана, який ми розглядаємо (лише взаємодії найближчих сусідів). На практиці ми вставляємо обертання Паулі $R_{xx}$, $R_{yy}$, $R_{zz}$ з параметризованим кутом $t$, які відповідають приблизній реалізації $e^{-i (XX + YY + ZZ) t}$. Враховуючи різницю в визначенні обертань Паулі та часової еволюції, яку ми намагаємося реалізувати, нам доведеться використовувати параметр $2*dt$ для досягнення часової еволюції $dt$. Крім того, ми змінюємо порядок операцій для непарного числа повторень кроків Троттера, що функціонально еквівалентно, але дозволяє синтезувати суміжні операції в один $SU(2)$ унітар. Це дає набагато менш глибоку схему, ніж та, що отримується за допомогою загальної функціональності 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)

Використання оптимізованої схеми для підготовки стану

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)

Шаблонні схеми для обчислення матричних елементів $\tilde{S}$ і $\tilde{H}$ за допомогою тесту Адамара

Єдина відмінність між схемами, що використовуються в тесті Адамара, буде полягати у фазі в операторі часової еволюції та вимірюваних спостережуваних. Тому ми можемо підготувати шаблонну схему, яка представляє загальну схему для тесту Адамара, з заповнювачами для вентилів, які залежать від оператора часової еволюції.

# 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

Ми значно зменшили глибину тесту Адамара за допомогою комбінації наближення Троттера та некерованих унітарів

Крок 3: Виконання з використанням примітивів Qiskit

Створення екземпляра бекенду та встановлення параметрів середовища виконання

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)

Транспіляція на QPU

Спочатку виберемо підмножини карти зв'язків з "добре" працюючими кубітами (де "добре" тут досить довільне, ми переважно хочемо уникнути дійсно погано працюючих кубітів) та створимо нову ціль для транспіляції

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

Потім транспілюємо віртуальну схему до найкращого фізичного розташування в цій новій цілі

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]

Створення PUB для виконання з Estimator

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

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

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

observables_S = [[observable_S_real], [observable_S_imag]]

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

layout = qc_trans.layout

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

observables_H = observable_trans_list

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

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

Виконання схем

Схеми для $t=0$ є класично обчислюваними

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)

Виконання схем для $S$ і $\tilde{H}$ за допомогою Estimator

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

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

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

job = estimator.run([pub])

Крок 4: Пост-обробка та повернення результату в бажаному класичному форматі

results = job.result()[0]

Обчислення ефективного гамільтоніана та матриць перекриття

Спочатку обчисліть фазу, накопичену станом $\vert 0 \rangle$ під час неконтрольованої часової еволюції

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

Отримавши результати виконання схем, ми можемо пост-обробити дані для обчислення матричних елементів $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]$$

І матричні елементи $\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]$$

Нарешті, ми можемо розв'язати узагальнену проблему власних значень для $\tilde{H}$:

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

і отримати оцінку енергії основного стану $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

Для сектора однієї частинки ми можемо ефективно обчислити основний стан цього сектора гамільтоніана класичним способом

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

Додаток: підпростір Крилова з еволюцій у реальному часі

Унітарний простір Крилова визначається як

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

для деякого часового кроку $dt$, який ми визначимо пізніше. Тимчасово припустимо, що $r$ є парним: тоді визначимо $d=r/2$. Зверніть увагу, що коли ми проєктуємо гамільтоніан у простір Крилова вище, він є нерозрізненним від простору Крилова

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

тобто, де всі часові еволюції зміщені назад на $d$ часових кроків. Причина нерозрізненності полягає в тому, що матричні елементи

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

є інваріантними щодо загальних зсувів часу еволюції, оскільки часові еволюції комутують з гамільтоніаном. Для непарного $r$ ми можемо використовувати аналіз для $r-1$.

Ми хочемо показати, що десь у цьому просторі Крилова гарантовано існує стан з низькою енергією. Ми робимо це за допомогою наступного результату, який виведено з теореми 3.1 у [3]:

Твердження 1: існує функція $f$ така, що для енергій $E$ у спектральному діапазоні гамільтоніана (тобто між енергією основного стану та максимальною енергією)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ для всіх значень $E$, які лежать на відстані $\ge\delta$ від $E_0$, тобто вона експоненційно придушена
3. $f(E)$ є лінійною комбінацією $e^{ijE\,dt}$ для $j=-d,-d+1,...,d-1,d$

Ми наводимо доведення нижче, але його можна сміливо пропустити, якщо Ви не хочете зрозуміти повний, строгий аргумент. Зараз ми зосередимося на наслідках наведеного вище твердження. З властивості 3 вище ми можемо побачити, що зміщений простір Крилова вище містить стан $f(H)|\psi\rangle$. Це наш стан з низькою енергією. Щоб зрозуміти чому, запишемо $|\psi\rangle$ в базисі власних енергій:

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

де $|E_k\rangle$ є k-им власним станом енергії, а $\gamma_k$ є його амплітудою в початковому стані $|\psi\rangle$. Виражений через це, $f(H)|\psi\rangle$ задається як

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

використовуючи той факт, що ми можемо замінити $H$ на $E_k$, коли він діє на власний стан $|E_k\rangle$. Тому похибка енергії цього стану становить

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

Щоб перетворити це на верхню межу, яку легше зрозуміти, ми спочатку розділимо суму в чисельнику на доданки з $E_k-E_0\le\delta$ та доданки з $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}.$$

Ми можемо обмежити зверху перший доданок значенням $\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,$$

де перший крок випливає з того, що $E_k-E_0\le\delta$ для кожного $E_k$ у сумі, а другий крок випливає з того, що сума в чисельнику є підмножиною суми в знаменнику. Для другого доданка спочатку обмежимо знаменник знизу значенням $|\gamma_0|^2$, оскільки $f(E_0)^2=1$: додавши все разом, це дає

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

Щоб спростити те, що залишилося, зверніть увагу, що для всіх цих $E_k$, за визначенням $f$, ми знаємо, що $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Додатково обмежуючи зверху $E_k-E_0<2\|H\|$ і обмежуючи зверху $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$, отримуємо

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

Це справедливо для будь-якого $\delta>0$, тому якщо ми встановимо $\delta$ рівним нашій цільовій похибці, то верхня межа похибки вище збігається до неї експоненційно з розмірністю Крилова $2d=r$. Також зверніть увагу, що якщо $\delta<E_1-E_0$, то доданок $\delta$ фактично повністю зникає у наведеній вище межі.

Щоб завершити аргумент, ми спочатку зазначимо, що наведене вище є лише похибкою енергії конкретного стану $f(H)|\psi\rangle$, а не похибкою енергії стану з найнижчою енергією в просторі Крилова. Однак, згідно з варіаційним принципом (Релея-Рітца), похибка енергії стану з найнижчою енергією в просторі Крилова обмежена зверху похибкою енергії будь-якого стану в просторі Крилова, тому наведене вище також є верхньою межею похибки енергії стану з найнижчою енергією, тобто виходу алгоритму квантової діагоналізації Крилова.

Подібний аналіз, як наведений вище, може бути проведений з додатковим урахуванням шуму та процедури порогового фільтрування, обговорюваної в зошиті. Див. [2] та [4] для цього аналізу.

Додаток: доведення твердження 1

Наступне в основному виведено з [3], теореми 3.1: нехай $0 < a < b$ і нехай $\Pi^*_d$ є простором залишкових поліномів (поліномів, значення яких у 0 дорівнює 1) степеня не більше $d$. Розв'язком

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

є

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

і відповідне мінімальне значення дорівнює

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

Ми хочемо перетворити це у функцію, яку можна природно виразити через комплексні експоненти, оскільки це еволюції в реальному часі, які генерують квантовий простір Крилова. Для цього зручно ввести наступне перетворення енергій у спектральному діапазоні гамільтоніана на числа в діапазоні $[0,1]$: визначимо

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

де $dt$ є часовим кроком таким, що $-\pi < E_0dt < E_\text{max}dt < \pi$. Зверніть увагу, що $g(E_0)=0$ і $g(E)$ зростає, коли $E$ віддаляється від $E_0$.

Тепер, використовуючи поліном $p^*(x)$ з параметрами a, b, d, встановленими як $a = g(E_0 + \delta)$, $b = 1$ і d = int(r/2), ми визначаємо функцію:

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

де $E_0$ є енергією основного стану. Ми можемо побачити, підставивши $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$, що $f(E)$ є тригонометричним поліномом степеня $d$, тобто лінійною комбінацією $e^{ijE\,dt}$ для $j=-d,-d+1,...,d-1,d$. Крім того, з визначення $p^*(x)$ вище ми маємо, що $f(E_0)=p(0)=1$ і для будь-якого $E$ у спектральному діапазоні такого, що $\vert E-E_0 \vert > \delta$, ми маємо

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

Посилання

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

Опитування про підручник

Будь ласка, пройдіть це коротке опитування, щоб надати відгук про цей підручник. Ваші ідеї допоможуть нам покращити наш контент і користувацький досвід.

Link to survey