Skip to main content
IBM Quantum Platform

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

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

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

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

0. クリロフ空間

次数のクリロフ空間 Kr\mathcal{K}^r rr は、行列 AA の高次のべき乗 r−1r-1 までを参照ベクトル ∣v⟩\vert v \rangle と掛け合わせることによって得られるベクトルによってスパンされる空間である。

Kr={∣v⟩,A∣v⟩,A2∣v⟩,...,Ar−1∣v⟩}\mathcal{K}^r = \left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}

行列 AA がハミルトニアン HH の場合、対応する空間はべき乗クリロフ空間 KP\mathcal{K}_P と呼ばれる。 AA がハミルトニアン U=e−iH(dt)U=e^{-iH(dt)} によって生成される時間発展演算子である場合、その空間はユニタリー・クリロフ空間 KU\mathcal{K}_U と呼ばれる。パワークリロフ部分空間は、 HH がユニタリー演算子ではないため、量子コンピュータ上で直接生成することはできません。 その代わりに、時間発展演算子 U=e−iH(dt)U = e^{-iH(dt)}、パワークリロフ空間と同様の収束保証が得られることを示すことができる。 UU のべき乗は、異なる時間ステップ Uk=e−iH(kdt)U^k = e^{-iH(k dt)} となり、 k=0,1,2,...,(r−1)k = 0, 1, 2, ..., (r-1).

KUr={∣ψ⟩,U∣ψ⟩,U2∣ψ⟩,...,Ur−1∣ψ⟩}\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=22L = 22 サイトを持つ反強磁性XX-Z spin-1/2 鎖のハミルトニアンを周期境界条件付きで考える：

H=∑i,jNJxy(XiXj+YiYj)+ZiZj H = \sum_{i, j}^{N} J_{xy} (X_{i} X_{j} + Y_{i} Y_{j}) + Z_{i} Z_{j} 
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つの主要な材料が必要だ：

  1. クリロフ次元( rr )と時間ステップ( dtdt )の選択。
  2. ターゲット（地上）状態と多項式にオーバーラップする初期（参照）状態（ベクトル ∣v⟩\vert v \rangle 上記）。ターゲット状態はスパースである。 この多項式オーバーラップの要件は、量子位相推定アルゴリズムと同じである。
  3. 時間発展演算子 Uk=e−iH(k∗dt)U^{k}=e^{-iH(k * dt)} ( k=0,1,2,...,r−1k = 0, 1, 2, ..., r-1 ).

rr （および、 dtdt ）の選択された値に対して、 rr 量子回路を作成し、それらからサンプルを採取する。 各量子回路は、参照状態の量子回路表現と、 kk 値の時間発展演算子を結合することで作成される。

クリロフ次元が大きいほど、推定エネルギーの収束性が向上する。 このレッスンでは、収束の傾向を説明するために、ディメンションを 55 に設定する。

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

クリロフ次元と時間ステップの他に、時間発展のためのトロッターステップ数を設定する必要がある。 ステップ数が少なすぎるとトロッター化の誤差が大きくなり、ステップ数が多すぎると回路が深くなる。 このレッスンでは、トロッターのステップ数を 66 に設定する。

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

次に、基底状態と重なる参照状態（ ∣ψ⟩\vert \psi \rangle ）を選ぶ必要がある。 このハミルトニアンでは、 1s と 0s ∣...101...010...101⟩\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の合成モジュールには、このような技法が数多く用意されています。 合成には 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:

Output of the previous code cell 
circuits[2].decompose().draw("mpl", fold=-1)

Output:

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)

次に、計算されたエネルギーをクリロフ次元の関数としてプロットし、厳密エネルギーと比較する。 正確なエネルギーは、古典的なブルートフォース法を用いて別途計算される。 推定される基底状態のエネルギーは、クリロフ空間の次元が大きくなるにつれて収束していくことがわかる。 55 のクリロフ次元は制限的であるが、それでも結果は印象的な収束を示しており、クリロフ次元が大きくなるにつれて改善することが期待される [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:

Output of the previous code cell

理解度チェック

下の問題を読んで、答えを考え、三角形をクリックして解答を明らかにしよう。

  • 回答:

    クリロフ次元を上げる。 一般的には、シュート数を増やすこともできるが、上記の計算ではすでにかなり多くなっている。

  • 回答:

    他にも有効な答えがあるかもしれないが、完全な答えには以下を含めること：

    (a) SKQDにはSQDにはない収束保証がある。 SQDでは、計算基底における基底状態のサポートと非常によく重なるようなアナザッツを推測するか、アナザッツのファミリーをサンプリングするために変分要素を計算に導入する必要があります。

    (b)SKQDは、ハダマルドテストによる行列要素の計算を回避できるため、QPUの所要時間が大幅に短縮される。

5. 概要

  • クリロフ基底状態のサンプリングによる基底状態のエネルギー推定は、スピン系、凝縮系問題、格子ゲージ理論などの格子モデルに非常に適している。 VQEやヒューリスティック・アンサッツに基づくSQD（例えば、前のレッスンの化学問題）のように、変分アサッツで多くのパラメータを最適化する必要がないため、このアプローチはVQEよりもはるかに優れたスケールを持つ。
    • 回路の深さを低く抑えるには、フォールト・トレラント・ハードウェアに適した格子問題に取り組むのが賢明である。
  • SKQDでは、VQEのような量子測定の問題は発生しない。 推定すべきパウリ作用素のグループはない。
  • SKQDはノイズの多いサンプルに対しても頑健である。なぜなら、問題固有のポストセレクションルーチン（例えば、問題固有のパターンに合致しないビット列を除外する）を使用したり、従来の対角化のオーバーヘッド（すなわち、より大きな部分空間で対角化を行う）を許容したりすることで、ノイズの影響を効果的に除去できるからである。

参照

[1] ジェフェリー・ユーほか "サンプルベースクリロフ対角化のための量子中心アルゴリズム" (2025). arxiv:quant-ph/2501.09702.

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

[2] N. 波多野とM. 鈴木, "高次の指数積公式の発見" (2005). arXiv:math-ph/0506007.

[4] D. ベリー、G.アホカス、R.クリーブ、B. Sanders, "Efficient quantum algorithms for simulating sparse Hamiltonians" (2006). arXiv:quant-ph/0508139.

このページは役に立ちましたか？
バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。