格子ハミルトニアンのクリロフ量子対角化

使用時間の目安: Heron r2で20分（注: これは目安です。実行時間は異なる場合があります。）

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

背景

このチュートリアルでは、Qiskitパターンの文脈でクリロフ量子対角化アルゴリズム（KQD）を実装する方法を示します。まずアルゴリズムの背後にある理論を学び、その後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: 古典的な入力を量子問題にマッピングする

クリロフ空間

次数 $r$ のクリロフ空間 $\mathcal{K}^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}$$

$\tilde{H}$ の固有値と固有状態を調べることで、$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)

次に、アルゴリズムのその他のパラメータを設定します。このチュートリアルでは、5次元のクリロフ空間のみを使用するように制限しますが、これはかなり制約があります。

# 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: 量子ハードウェア実行のための問題の最適化

効率的なアダマールテスト

いくつかの近似を導入し、モデルハミルトニアンに関するいくつかの仮定を利用することで、アダマールテストの深い回路を最適化することができます。例えば、以下のアダマールテスト回路を考えてみましょう。

ハミルトニアン $H$ のもとでの $|0\rangle^N$ の固有値 $E_0$ を古典的に計算できると仮定します。 この条件は、ハミルトニアンが 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}$$

ここで、3行目では古典的にシミュレーション可能な位相シフト $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$ のみを実装すればよいのです。上記のように計算を再構成することで、結果として得られる回路の深さを大幅に削減できます。

トロッター分解による時間発展演算子の分解

時間発展演算子を厳密に実装する代わりに、トロッター分解を使用してその近似を実装することができます。ある次数のトロッター分解を複数回繰り返すことで、近似によって導入される誤差をさらに低減できます。以下では、考慮しているハミルトニアンの相互作用グラフ（最近接相互作用のみ）に対して最も効率的な方法で直接トロッター実装を構築します。実際には、パラメータ化された角度 $t$ を持つパウリ回転 $R_{xx}$、$R_{yy}$、$R_{zz}$ を挿入し、これは $e^{-i (XX + YY + ZZ) t}$ の近似的な実装に対応します。パウリ回転と実装しようとしている時間発展の定義の違いにより、$dt$ の時間発展を達成するためにはパラメータ $2*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]

Estimator を使用した実行のための PUB の作成

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

Estimator を使用して $S$ および $\tilde{H}$ の回路を実行します。

# 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] の定理3.1から導かれる以下の結果を用いて示します。

主張1: ハミルトニアンのスペクトル範囲内のエネルギー $E$（すなわち、基底状態エネルギーと最大エネルギーの間）に対して、以下の性質を持つ関数 $f$ が存在します。

1. $f(E_0)=1$
2. $E_0$ から $\ge\delta$ 離れたすべての $E$ の値に対して $|f(E)|\le2\left(1 + \delta\right)^{-d}$、すなわち指数的に抑制されます
3. $f(E)$ は $j=-d,-d+1,...,d-1,d$ に対する $e^{ijE\,dt}$ の線形結合です

証明は以下に示しますが、完全で厳密な議論を理解したい場合を除き、安全にスキップできます。ここでは上記の主張の意味に焦点を当てます。上記の性質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\rangle$ に作用するとき $H$ を $E_k$ に置き換えられるという事実を利用しました。したがって、この状態のエネルギー誤差は次のようになります。

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

第1項は $\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,$$

ここで、第1ステップは和の中のすべての $E_k$ に対して $E_k-E_0\le\delta$ であることから従い、第2ステップは分子の和が分母の和の部分集合であることから従います。第2項については、$f(E_0)^2=1$ であるから分母を $|\gamma_0|^2$ で下から抑えます。すべてをまとめると、次が得られます。

$$\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$ を次数が高々 $d$ の残差多項式（0での値が1である多項式）の空間とします。以下の問題の解は

$$\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$ であり、$E$ が $E_0$ から離れるにつれて $g(E)$ は増加することに注意してください。

次に、多項式 $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$ の三角多項式、すなわち $j=-d,-d+1,...,d-1,d$ に対する $e^{ijE\,dt}$ の線形結合であることがわかります。さらに、上記の $p^*(x)$ の定義から、$f(E_0)=p(0)=1$ であり、$\vert E-E_0 \vert > \delta$ を満たすスペクトル範囲内の任意の $E$ に対して次が成り立ちます。

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

チュートリアルアンケート

このチュートリアルに関するフィードバックをお寄せいただくため、こちらの短いアンケートにご協力ください。皆様のご意見は、コンテンツの改善およびユーザー体験の向上に役立てさせていただきます。