サンプルベースのクリロフ量子対角化（SKQD）

このサンプルベースのクリロフ量子対角化（SKQD）のレッスンでは、前のレッスンで説明した手法を組み合わせます。Qiskit パターンフレームワークを活用した単一の例から構成されています：

ステップ 1: 問題を量子回路と演算子にマッピングする
ステップ 2: ターゲットハードウェア向けに最適化する
ステップ 3: Qiskit プリミティブを使用して実行する
ステップ 4: 後処理を行う

サンプルベースの量子対角化法において重要なステップは、部分空間のための高品質なベクトルを生成することです。前のレッスンでは、化学ハミルトニアンの部分空間ベクトルを生成するために LUCJ アンサッツを使用しました。このレッスンでは、レッスン 2 で議論したように、量子クリロフ状態[1]を使用します。まず、時間発展演算子を使って量子コンピュータ上でクリロフ空間を生成する方法を復習します。次に、そこからサンプリングを行います。サンプリングされた部分空間にシステムハミルトニアンを射影し、対角化することで基底状態エネルギーを推定します。このアルゴリズムは、レッスン 2 で説明した仮定のもと、基底状態へ確実かつ効率的に収束することが証明されています。

0. クリロフ空間

次数 $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^{-iH(dt)}$ の場合、その空間はユニタリクリロフ空間 $\mathcal{K}_U$ と呼ばれます。冪クリロフ部分空間は、 $H$ がユニタリ演算子でないため、量子コンピュータ上で直接生成することができません。代わりに、時間発展演算子 $U = e^{-iH(dt)}$ を使用することができ、これは冪クリロフ空間と同様の収束保証を与えることが示されています。 $U$ の累乗は、 $k = 0, 1, 2, ..., (r-1)$ として、異なる時間ステップ $U^k = e^{-iH(k dt)}$ になります。

\mathcal{K}_U^r = \left\{ \vert \psi \rangle, U \vert \psi \rangle, U^2 \vert \psi \rangle, ..., U^{r-1} \vert \psi \rangle \right\}

1. 問題を量子回路と演算子にマッピングする

このレッスンでは、周期境界条件を持つ $L = 22$ サイトのスピン-1/2 反強磁性 XX-Z 鎖のハミルトニアンを考えます：

H = \sum_{i, j}^{N} J_{xy} (X_{i} X_{j} + Y_{i} Y_{j}) + Z_{i} Z_{j}

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-sqd qiskit-addon-utils qiskit-ibm-runtime

from qiskit.transpiler import CouplingMap
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

num_spins = 22
coupling_map = CouplingMap.from_ring(num_spins)
H_op = generate_xyz_hamiltonian(coupling_map, coupling_constants=(0.3, 0.3, 1.0))

クリロフ空間を構築するには、3 つの主要な要素が必要です：

クリロフ次元（ $r$ ）と時間ステップ（ $dt$ ）の選択。
目標（基底）状態と多項式オーバーラップを持つ初期（参照）状態（上記のベクトル $\vert v \rangle$ ）。ただし、目標状態はスパースである必要があります。この多項式オーバーラップの要件は、量子位相推定アルゴリズムと同じです。
時間発展演算子 $U^{k}=e^{-iH(k * dt)}$ （ $k = 0, 1, 2, ..., r-1$ ）。

選択した $r$ の値（および $dt$ ）に対して、 $r$ 個の別々の量子回路を作成し、そこからサンプリングします。各量子回路は、参照状態の量子回路表現と $k$ の値に対する時間発展演算子を結合して作成されます。

クリロフ次元を大きくすると、推定エネルギーの収束が改善されます。このレッスンでは、収束の傾向を示すために次元を $5$ に設定します。

参考文献[2]は、KQD に対して十分に小さい時間ステップが $\pi / \vert \vert H \vert \vert$ であり、この値を過大評価するよりも過小評価する方が望ましいことを示しました。一方、 $dt$ を小さくしすぎると、クリロフ部分空間のコンディショニングが悪化します。これはクリロフ基底ベクトルが時間ステップごとに変化しにくくなるためです。さらに、この $dt$ の選択は SKQD の収束に対して証明可能に十分ですが、このサンプリングベースの文脈では、実際の最適な $dt$ の選択は現在も研究が進められているテーマです。このレッスンでは $dt = 0.15$ を設定します。

クリロフ次元と時間ステップに加えて、時間発展のトロッター分解のステップ数を設定する必要があります。ステップ数が少なすぎると大きなトロッター化誤差が生じ、多すぎると回路が深くなります。このレッスンでは、トロッター分解のステップ数を $6$ に設定します。

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

次に、基底状態とある程度のオーバーラップを持つ参照状態 $\vert \psi \rangle$ を選ぶ必要があります。このハミルトニアンには、1 と 0 が交互に並ぶ ネール 状態 $\vert ...101...010...101 \rangle$ を参照状態として使用します。

# Prep `Neel` state as the reference state for evolution
from qiskit import QuantumCircuit

qc_state_prep = QuantumCircuit(num_spins)
for i in range(num_spins):
    if i % 2 == 0:
        qc_state_prep.x(i)

最後に、時間発展演算子を量子回路にマッピングする必要があります。これはレッスン 2 で行いましたが、ここでは Qiskit のメソッド、特に synthesis（合成）と呼ばれるメソッドを活用します。数学的演算子を量子ゲートによる量子回路に合成するための様々な手法があります。そのような手法の多くは Qiskit synthesis モジュールで利用できます。合成には LieTrotter アプローチを使用します[3] [4]。

from qiskit.circuit import QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

evol_gate = PauliEvolutionGate(
    H_op, time=(dt / num_trotter_steps), synthesis=LieTrotter(reps=num_trotter_steps)
)  # `U` operator

qr = QuantumRegister(num_spins)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

circuits = []
for rep in range(krylov_dim):
    circ = qc_state_prep.copy()

    # Repeating the `U` operator to implement U^0, U^1, U^2, and so on, for power Krylov space
    for _ in range(rep):
        circ.compose(other=qc_evol, inplace=True)

    circ.measure_all()
    circuits.append(circ)

circuits[1].decompose().draw("mpl", fold=-1)

Output of the previous code cell

circuits[2].decompose().draw("mpl", fold=-1)

Output of the previous code cell

2. ターゲットハードウェア向けに最適化する

回路を作成したら、ターゲットハードウェア向けに最適化できます。ユーティリティスケールの QPU を選択します。

import warnings

from qiskit import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService

warnings.filterwarnings("ignore")

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

次に、プリセットパスマネージャーを使用して回路をターゲットバックエンド向けにトランスパイルします。

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_circuits = pm.run(circuits=circuits)

3. ターゲットハードウェアで実行する

ハードウェア実行向けに回路を最適化した後、ターゲットハードウェア上で実行し、基底状態エネルギー推定のためのサンプルを収集する準備が整います。

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
job = sampler.run(isa_circuits, shots=100_000)  # Takes approximately 2m 58s of QPU time

counts_all = [job.result()[k].data.meas.get_counts() for k in range(krylov_dim)]

4. 結果を後処理する

次に、増加するクリロフ次元のカウントを累積的に集計します。累積カウントを使用して、クリロフ次元が増加する部分空間を張り、収束の振る舞いを分析します。

from collections import Counter

counts_cumulative = []
for i in range(krylov_dim):
    counter = Counter()
    for d in counts_all[: i + 1]:
        counter.update(d)

    counts = dict(counter)
    counts_cumulative.append(counts)

ハミルトニアンを射影・対角化するために、qiskit-addon-sqd の機能を使用します。このアドオンは、パウリ文字列ベースのハミルトニアンを部分空間に射影し、SciPy を用いて固有値を求める機能を提供します。

from qiskit_addon_sqd.counts import counts_to_arrays
from qiskit_addon_sqd.qubit import solve_qubit

原理的には、部分空間を張る前に、誤ったパターンのビット列をフィルタリングすることができます。例えば、このレッスンの反強磁性ハミルトニアンの基底状態は、通常「アップ」スピンと「ダウン」スピンが等数となります。つまり、ビット列の "1" の数は、システムの総ビット数（スピン数）のちょうど半分であるべきです。次の関数は、"1" の数が指定と異なるビット列をカウントからフィルタリングします。

# Filters out bitstrings that do not have specified number (`num_ones`) of `1` bits.
def postselect_counts(counts, num_ones):
    filtered_counts = {}
    for bitstring, freq in counts.items():
        if bitstring.count("1") == num_ones:
            filtered_counts[bitstring] = freq

    return filtered_counts

正しい数のアップ/ダウン電子を持つビット列を使用して、部分空間を張り、クリロフ次元の増加に伴う固有値を計算します。問題のサイズと利用可能な古典的リソースによっては、部分空間の次元を適切に抑えるためにサブサンプリング（SQD に関するレッスンと同様）を採用する必要があるかもしれません。さらに、レッスン 4 と同様の配置回復の概念を適用することもできます。再構築された固有状態からサイトごとの電子占有率を計算し、その情報を使ってアップ/ダウン電子数が不正確なビット列を補正することができます。これは興味のある読者への演習問題として残しておきます。

import numpy as np

num_batches = 10
rand_seed = 0
scipy_kwargs = {"k": 2, "which": "SA"}

ground_state_energies = []
for idx, counts in enumerate(counts_cumulative):
    counts = postselect_counts(counts, num_ones=num_spins // 2)
    bitstring_matrix, probs = counts_to_arrays(counts=counts)

    eigenvals, eigenstates = solve_qubit(
        bitstring_matrix, H_op, verbose=False, **scipy_kwargs
    )
    gs_en = np.min(eigenvals)
    ground_state_energies.append(gs_en)

次に、クリロフ次元の関数として計算されたエネルギーをプロットし、厳密なエネルギーと比較します。厳密なエネルギーは、古典的なブルートフォース法を用いて別途計算されています。推定された基底状態エネルギーがクリロフ空間の次元の増加とともに収束することがわかります。クリロフ次元 $5$ は制限的ですが、それでも結果は印象的な収束を示しており、より大きなクリロフ次元では改善が期待されます[1]。

import matplotlib.pyplot as plt

exact_gs_en = -23.934184
plt.plot(
    range(1, krylov_dim + 1),
    ground_state_energies,
    color="blue",
    linestyle="-.",
    label="estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [exact_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.ylim([-24, -22.50])
plt.title(
    "Estimating Ground state energy with Sample-based Krylov Quantum Diagonalization"
)
plt.show()

Output of the previous code cell

理解度チェック

以下の質問を読み、答えを考えてから、三角形をクリックして解答を確認してください。

上のプロットの収束を改善するためにできることは何ですか？

解答：

クリロフ次元を増やすことです。一般的には、ショット数を増やすこともできますが、上記の計算ではすでにかなり高い値が設定されています。

SKQD が（a）SQD、および（b）KQD に対して持つ主な利点は何ですか？

解答：

他にも有効な答えがあるかもしれませんが、完全な解答には以下を含める必要があります：

(a) SKQD には SQD にはない収束保証があります。SQD では、基底状態の計算基底サポートとの優れたオーバーラップを持つアンサッツを非常に上手く推測するか、または答えの候補となるアンサッツの族をサンプリングするために変分的なコンポーネントを計算に導入する必要があります。

(b) SKQD は、アダマールテストを通じた行列要素の計算というコストのかかる処理を回避するため、必要な QPU 時間がはるかに少なくて済みます。

5. まとめ

クリロフ基底状態のサンプリングによる基底状態エネルギー推定は、スピン系、物性問題、格子ゲージ理論などの格子モデルに非常に適しています。このアプローチは VQE よりもはるかに良くスケールします。なぜなら、VQE のように変分アンサッツの多くのパラメータを最適化したり、ヒューリスティックなアンサッツベースの SQD（例えば、前のレッスンの化学問題）のような最適化が不要だからです。
- 回路の深さを低く抑えるため、プレフォールトトレラントハードウェアに適した格子問題に取り組むのが賢明です。
SKQD は VQE のような量子測定問題が生じません。推定すべき可換パウリ演算子のグループがないためです。
SKQD はノイズの多いサンプルに対してロバストです。問題固有のポスト選択ルーティン（例えば、問題固有のパターンに従わないビット列をフィルタリングする）を使用したり、古典的な対角化のオーバーヘッドを許容（つまり、より大きな部分空間で対角化する）することで、ノイズの影響を効果的に除去できます。

参考文献

[1] Jeffery Yu et al., "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization" (2025). arxiv:quant-ph/2501.09702.

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

[2] N. Hatano and M. Suzuki, "Finding Exponential Product Formulas of Higher Orders" (2005). arXiv:math-ph/0506007.

[4] D. Berry, G. Ahokas, R. Cleve and B. Sanders, "Efficient quantum algorithms for simulating sparse Hamiltonians" (2006). arXiv:quant-ph/0508139.

0. クリロフ空間​

1. 問題を量子回路と演算子にマッピングする​

2. ターゲットハードウェア向けに最適化する​

3. ターゲットハードウェアで実行する​

4. 結果を後処理する​

理解度チェック​

5. まとめ​

参考文献​

0. クリロフ空間

1. 問題を量子回路と演算子にマッピングする

2. ターゲットハードウェア向けに最適化する

3. ターゲットハードウェアで実行する

4. 結果を後処理する

理解度チェック

5. まとめ

参考文献