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

使用時間の目安：ヘロンで20分 r2 （注：あくまでも目安です。ランタイムは異なるかもしれない)。

背景

このチュートリアルでは、クリロフ量子対角化アルゴリズム（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)

    # 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

ステップ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}

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

Output:

[('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

Output:

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)

Output:

時間発展

我々は、与えられたハミルトニアンによって生成される時間発展演算子を実現することができる。 $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)

Output:

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

Output:

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

Output:

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$ ）だけを実装すればよく、制御された時間発展は必要ない。上記のように計算を組み直すことで、結果として得られる回路の深さを大幅に減らすことができる。

トロッター分解を用いて時間発展演算子を分解する

時間発展演算子を正確に実装する代わりに、トロッター分解を使ってその近似を実装することができる。ある次数のトロッター分解を数回繰り返すことで、近似からもたらされる誤差をさらに減らすことができる。以下では、我々が考えているハミルトニアンの相互作用グラフ（最近傍相互作用のみ）に対して、最も効率的な方法でトロッター実装を直接構築する。実際には、パウリ回転 $R_{xx}$ 、 $R_{yy}$ 、 $R_{zz}$ を、 $e^{-i (XX + YY + ZZ) t}$ の近似実装に対応するパラメトリック角度 $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)

Output:

状態準備に最適化された回路を使用する

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)

Output:

$\tilde{S}$ および $\tilde{H}$ の行列要素をHadamardテストを用いて計算するためのテンプレート回路

ハダマード・テストで使用される回路の違いは、時間発展演算子の位相と測定される観測値だけである。したがって、時間発展演算子に依存するゲートのプレースホルダーを持つ、ハダマード・テストの一般的な回路を表すテンプレート回路を用意することができる。

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

Output:

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

Output:

The optimized circuit has 2Q gates depth:  74

我々は、トロッター近似と非制御ユニタリーの組み合わせにより、ハダマード・テストの深さを大幅に削減した

ステップ3: `Qiskit primitives`を使用して実行する

バックエンドをインスタンス化し、ランタイムパラメータを設定する

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

Output:

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)

Output:

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

Output:

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

Output:

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

Output:

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)

Output:

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

Output:

付録：実時間発展からのクリロフ部分空間

ユニタリー・クリロフ空間は次のように定義される

\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$ が存在する。

$f(E_0)=1$
$|f(E)|\le2\left(1 + \delta\right)^{-d}$ $E_0$ から離れたところにあるのすべての値に対して、つまり指数関数的に抑制される。 $\ge\delta$ $E$
$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,

が固有状態 $|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}.

第一項を $\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_k-E_0\le\delta$ 、2番目のステップは、分子の和が分母の和の部分集合であるため、続く。第2項については、まず分母を $|\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$ のエネルギー誤差に過ぎないことにまず注意する。しかし、(Rayleigh-Ritz) 変分原理により、クリロフ空間内の最低エネルギー状態のエネルギー誤差は、クリロフ空間内の任意の状態のエネルギー誤差によって上界されるため、上記は最低エネルギー状態のエネルギー誤差、つまりクリロフ量子対角化アルゴリズムの出力に対する上界でもある。

上記と同様の分析が、ノイズとノートで説明した閾値処理手順を追加して実施できる。この分析については [2] と [4] を参照のこと。

付録：請求項１の証明

以下はほとんど [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]$ ：defineの範囲内の数に変換する以下の変換を導入するのが便利である

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

ここで $dt$ は $-\pi < E_0dt < E_\text{max}dt < \pi$ のようなタイムステップである。 $E$ が $E_0$ から遠ざかるにつれて、 $g(E_0)=0$ と $g(E)$ が成長することに注意。

ここで、パラメータa, b, dを $a = g(E_0 + \delta)$ , $b = 1$ , d = int( r/2 ) に設定した多項式 $p^*(x)$ を使い、関数を定義する：

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$ は基底状態のエネルギーである。 $f(E)$ が次数 $d$ の三角多項式、つまり $j=-d,-d+1,...,d-1,d$ に対する $e^{ijE\,dt}$ の線形結合であることは、 $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ を挿入することでわかる。さらに、上記の $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. 吉岡、M.アミコ、W.カービーら、"量子プロセッサー上での大規模多体ハミルトニアンの対角化"。 arXiv:2407.14431

[2] イーサン・エッペリー、リン・リン、中務裕司。「量子部分空間対角化の理論」。 SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).

[3] Å. ビョルク「行列計算における数値的手法」。応用数学のテキスト。シュプリンガー・インターナショナル・パブリッシング (2014).

[4] ウィリアム・カービー「誤差を含む量子クリロフ・アルゴリズムの解析」。 Quantum 8, 1457 (2024).

チュートリアル調査

このチュートリアルに関するご意見をお聞かせください。あなたの洞察は、私たちのコンテンツの提供とユーザーエクスペリエンスを向上させるのに役立ちます。

アンケートへのリンク

このページは役に立ちましたか？

バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。

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

背景

要件

セットアップ

ステップ1：古典的な入力を量子問題にマッピングする

クリロフ空間

クリロフ量子対角化アルゴリズム

ハミルトニアンを定義する

アルゴリズムのパラメータを設定する

状態準備

時間発展

ハダマール検定

ステップ2：量子ハードウェア実行に向けた問題の最適化

効率的なハダマール検定

トロッター分解を用いて時間発展演算子を分解する

状態準備に最適化された回路を使用する

S~\tilde{S}S~ および H~\tilde{H}H~ の行列要素をHadamardテストを用いて計算するためのテンプレート回路

ステップ3: Qiskit primitivesを使用して実行する

QPUへのトランスパイル

Estimatorによる実行用のPUBを作成する

回路を実行する

ステップ4：後処理を行い、結果を希望の古典形式で返す

有効ハミルトニアンとオーバーラップ行列を計算する

付録：実時間発展からのクリロフ部分空間

付録：請求項１の証明

参照

チュートリアル調査

$\tilde{S}$ および $\tilde{H}$ の行列要素をHadamardテストを用いて計算するためのテンプレート回路

ステップ3: `Qiskit primitives`を使用して実行する