量子汎用化（VQE）を用いたハイゼンベルグ鎖の基底状態エネルギー推定

推定実行時間：Heronプロセッサで37分（注：これはあくまで推定値です。 （実行時間は異なる場合があります。）

学習成果

このチュートリアルを完了すると、以下の内容を理解できるようになります：

• Qiskit を使用してハイゼンベルグスピン鎖を量子ハミルトニアンとしてモデル化する方法
• SPSAオプティマイザーを用いて量子系の基底状態エネルギーを推定する方法
• IBM Quantum のプリミティブとセッションを使用して、 IBM® の量子ハードウェア上で変分ワークフローを実行する方法

前提条件

以下のトピックについて、あらかじめ理解しておくことをお勧めします：

• 量子情報の基礎
• Qiskitパターンの概要
• 変分アルゴリズム設計

背景

ハイゼンベルグスピン鎖は、物性物理学および量子磁性学において最も広く研究されているモデルの一つである。 これは、相互作用する量子スピンの一次元格子について記述したもので、隣接するスピンは交換相互作用を通じて結合している。 外部磁場が加わった等方性ハイゼンベルクモデルのハミルトニアンは、次のように表される：

$H = \sum_{\langle i,j \rangle} \left( J_x X_i X_j + J_y Y_i Y_j + J_z Z_i Z_j \right) + \sum_{i} h_i Z_i,$

ここで、 $X_i$、 $Y_i$、および $Z_i$ は、サイト $i$ に作用するパウリ演算子であり、和 $\langle i,j \rangle$ は最近接ペアに対してとられ、 $J_x = J_y = J_z = 0.5$ は交換結合定数（本チュートリアルでは等方性）であり、 $h_i$ はサイトに依存する外部磁場を表す。 このチュートリアルでは、磁場値は $[-1, 1]$ の範囲からランダムにサンプリングされます。なお、以下の実装において、「最近傍」ペアの集合は、最初の $N$ 個の量子ビット間のハードウェアバックエンド固有の結合によって決定されますが、デバイスのトポロジーによっては、これらが厳密な線形チェーンを形成しない場合があります。

このハミルトニアンの基底状態エネルギーを理解することは、物理学において極めて重要である。 基底状態には、量子相転移、エンタングルメント構造、および磁気秩序に関する情報が含まれている。 従来、スピンの数が増加するにつれて、基底状態の正確なエネルギーを計算することは困難になります。これは、 $N$ 個のスピンを持つ場合、ヒルベルト空間の次元が $2^N$ という指数関数的なスケールで増加するためです。 このため、量子シミュレーションの最適な対象となる。

変分量子固有値ソルバー（VQE）は、ハミルトニアンの基底状態エネルギーを推定するために設計された、量子・古典ハイブリッドアルゴリズムである。 この手法は、量子コンピュータ上でパラメータ化された量子状態 $|\psi(\theta)\rangle$ （アンザッツと呼ばれる）を準備し、その期待値 $\langle \psi(\theta) | H | \psi(\theta) \rangle$ を測定することで機能する。その後、古典的な最適化アルゴリズムが、このエネルギーを最小化するようにパラメータ $\theta$ を反復的に調整する。この際、測定されたエネルギーが常に真の基底状態エネルギーの上限となることを保証する変分原理が活用される。

このチュートリアルでは、Qiskitのcircuitライブラリにあるansatz efficient_su2 を使用します。これは、単一量子ビットの回転ゲートとエンタングルメントゲートを層状に構成するものです。 この最適化には、同時摂動確率近似（SPSA）アルゴリズムが用いられている。このアルゴリズムは、パラメータの数にかかわらず、1回の反復ごとにわずか2回の関数評価のみで勾配を推定するため、ノイズの多い量子ハードウェアに特に適している。

要件

このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：

• Qiskit SDK v2.0 以降、 可視化機能を搭載
• Qiskit Runtime v0.44 またはそれ以降 (pip install qiskit-ibm-runtime)

セットアップ

import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimatorV2
from qiskit.circuit.library import XGate
from qiskit.circuit.library import efficient_su2
from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.scheduling import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2


def visualize_results(results):
    plt.plot(results["cost_history"], lw=2)
    plt.xlabel("Number of function evaluations")
    plt.ylabel("Energy")
    plt.show()

小規模な例

このセクションでは、Qiskitのパターンを小規模な例を用いて段階的に解説し、ワークフローを構築しながら主要な構成要素について説明します。

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

• 入力スピン数
• アウトプットハイゼンベルグ連鎖をモデル化するアンサッツとハミルトニアン

10スピンのハイゼンベルグ鎖をモデル化する近似解とハミルトニアンを構築せよ。 このステップでは、負荷が最も低いバックエンドの結合マップを用いて10スピンのハイゼンベルク・ハミルトニアンを構築し、アンザッツを efficient_su2 準備します。

num_spins = 10
ansatz = efficient_su2(num_qubits=num_spins, reps=2)

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=num_spins, simulator=False
)

coupling = backend.target.build_coupling_map()
reduced_coupling = coupling.reduce(list(range(num_spins)))

edge_list = reduced_coupling.graph.edge_list()
ham_list = []

for edge in edge_list:
    ham_list.append(("ZZ", edge, 0.5))
    ham_list.append(("YY", edge, 0.5))
    ham_list.append(("XX", edge, 0.5))

for qubit in reduced_coupling.physical_qubits:
    ham_list.append(("Z", [qubit], np.random.random() * 2 - 1))

hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)

ansatz.draw("mpl", style="iqp")

Output:

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

• 入力抽象回路、観測可能
• 出力：選択されたQPUに最適化されたターゲット回路と観測値

Qiskitの generate_preset_pass_manager 関数を使用して、選択したQPUに関する回路の最適化ルーチンを自動的に生成します。 プリセットパスマネージャーの中で最高レベルの最適化を提供する optimization_level=3。 また、 ALAPScheduleAnalysis 、 PadDynamicalDecoupling 、デコヒーレンス・エラーを抑制するためのスケジューリング・パスも含んでいる。

target = backend.target
pm = generate_preset_pass_manager(optimization_level=3, target=target)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(durations=target.durations()),
        PadDynamicalDecoupling(
            durations=target.durations(),
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)
isa_ansatz = pm.run(ansatz)
isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)
isa_ansatz.draw("mpl", scale=0.6, style="iqp", fold=-1, idle_wires=False)

Output:

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

• 入力対象回路と観測値
• アウトプット最適化の結果

回路パラメータを最適化することで、系の基底状態の推定エネルギーを最小化する。 最適化中にコスト関数を評価するには、 IBM QuantumEstimator プリミティブを使用します。

ステップ2でバックエンド向けに回路を最適化したため、を設定し skip_transpilation=True 、最適化された回路を渡すことで、Quantum Computeサーバーでのトランスパイル処理を回避できます。 このデモでは、 IBM Quantum プリミティブを使用してQPU上で実行します。 statevector ベース qiskit のプリミティブで実行するには、 IBM Quantum のプリミティブを使用しているコードブロックを、コメントアウトされたブロックに置き換えてください。

このチュートリアルでは、勾配法に基づく最適化アルゴリズムである同時摂動確率近似法（SPSA）を使用します。 次に、その概要を簡単に紹介し、Qiskit v2.0 を使用して SPSA を実装するためのコードを提示します。

SPSAのご紹介

同時摂動確率近似法（SPSA） [1] は、各反復においてわずか2回の関数呼び出しのみを用いて、勾配ベクトル全体を近似する最適化アルゴリズムである。 $f:\mathbb{R}^p\rightarrow \mathbb{R}$ を、最適化対象のパラメータ $p$ を持つコスト関数とし、 $x_i\in \mathbb{R}^p$ を反復計算の $i^{th}$ ステップにおけるパラメータベクトルとする。 勾配を計算するために、サイズ $p$ のランダムベクトル $\Delta_i$ が生成される。ここで、各要素 $\Delta_{ij}$、 $\forall$、 $j\in \{1,2,...,p\}$ は、 $\{-1, 1\}$ から一様分布に従ってサンプリングされる。次に、ランダムベクトル $\Delta_i$ の各要素に小さな値 $c_i$ を乗算し、ランダムな摂動を生成する。 その後、勾配は次のように推定される

$[\nabla f(x_i)]_j \approx \frac{f(x_i + c_i \Delta_i) - f(x_i - c_i \Delta_i)}{2c_i\Delta_{ij}}.$

直感的には、勾配推定の過程でランダムな外乱が加えられるため、ノイズに起因する $f$ の正確な値におけるわずかな偏差は許容され、考慮されることが期待される。 実際、SPSAはノイズに対する耐性が特に高いことで知られており、1回の反復ごとに必要なハードウェア呼び出しはわずか2回です。 したがって、これは変分アルゴリズムを実装する際に、最も好まれる最適化手法の一つである。

このチュートリアルでは、 $i^{th}$ の反復処理におけるハイパーパラメータである $a_i$ および $c_i$ は、次のように計算されます

$a_i = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{c}{(i+1)^\gamma},$

ここで、定数値は $A = 30$、 $\alpha = 0.9$、 $a = 0.3$、 $c = 0.1$、および $\gamma = 0.4$ とする。これらの値は [2] から選択されたものである。 SPSAから良好な性能を引き出すには、ハイパーパラメータを適切に調整する必要があります。

def spsa(
    fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100
):
    nparams = len(x0)
    x = np.copy(x0)

    for i in range(maxiter):
        a_i = a / (A + i + 1) ** alpha
        c_i = c / (i + 1) ** gamma
        delta_i = np.random.choice([-1, 1], nparams)

        # two hardware calls
        eval_1 = fun(x + c_i * delta_i, *args)
        eval_2 = fun(x - c_i * delta_i, *args)

        # compute the gradient and update the parameters
        grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)
        x = x - a_i * grad

    return x

def cost_func(
    params: Sequence,
    ansatz: QuantumCircuit,
    hamiltonian: SparsePauliOp,
    estimator: BaseEstimatorV2,
    cost_history_dict: dict,
) -> float:
    """Ground state energy evaluation."""
    energy = (
        estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs
    )

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = list(params)
    cost_history_dict["cost_history"].append(float(energy[0]))

    print(
        f"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]",
        end="\r",
    )

    return energy


def solve(x0, isa_ansatz, isa_observable, maxiter=150):
    cost_history_dict = {
        "prev_vector": None,
        "iters": 0,
        "cost_history": [],
        "y_min": None,
    }

    # Evaluate the problem using a QPU via `qiskit-ibm-runtime`
    with Session(backend=backend) as session:
        estimator = EstimatorV2(mode=session)
        estimator.skip_transpilation = True
        estimator.options.environment.job_tags = ["TUT_HSVQE"]
        x_opt = spsa(
            cost_func,
            x0=x0,
            args=(isa_ansatz, isa_observable, estimator, cost_history_dict),
            maxiter=maxiter,
        )

        y_min = cost_func(
            x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict
        )

    return y_min, cost_history_dict

np.random.seed(42)
num_params = ansatz.num_parameters
params = 2 * np.pi * np.random.random(num_params)

ここでは、 maxiter = 50. を設定します。 なお、勾配を計算するために各反復ごとに2回の関数呼び出しが必要となるため、関数呼び出しの総数は $2 \times \text{maxiter}$ 回となります。エネルギー推定の精度を高めるために、この値は任意のより大きな maxiter 値に増やすことができます。

maxiter = 50
spsa_min, spsa_history = solve(
    params, isa_ansatz, isa_observable, maxiter=maxiter
)

Output:

Fx Iters. done: 101 [Current cost: -3.03843]

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

• 入力：最適化中の基底状態エネルギーの推定値
• 出力：基底状態の推定エネルギー

print(f"Estimated ground state energy: {spsa_min}")

Output:

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Output:

大規模なハードウェアの例

このチュートリアルには、大規模なハードウェアの例は含まれていません。 量子ビットの数が増えるにつれ、VQEは「 バーレン・プラトー 現象」により大きな課題に直面する。これは、コスト関数の勾配がシステムの規模とともに指数関数的に消失するため、大規模な回路では最適化が事実上不可能になるという現象である。 これにハードウェアノイズが加わるため、VQEをより長いスピン鎖に拡張しても、信頼性の高い再現性のある結果が得られないことになる。 これらの制限を克服する手法については、以下の「次のステップ」のセクションを参照してください。

課題

ハイゼンベルグ連鎖のVQE実装が動作するようになったので、次のことを試してみてください：

1. reps=3アンザッツの深さを変えて試してみてください： （例えば、と efficient_su2``reps を試 reps=1 してみてください）。 アンザッツの深さは、基底状態のエネルギーの推定値や収束速度にどのような影響を与えるのでしょうか？ どの時点で、収益の逓減や不安定さが見られるようになりますか？
2. SPSAのハイパーパラメータを調整する： 学習率のスケジューリングパラメータ（a, c, alpha, gamma, A）を調整し、収束にどのような影響を与えるかを確認する。 ここで使用されているデフォルト設定よりも収束が速い設定を見つけることはできますか？
3. 結合トポロジーの比較： バックエンドのネイティブ結合マップを使用する代わりに、単純な最近傍線形チェーンを構築して、その結果を比較してみてください。 物理ハードウェアの接続性は、トランスパイルされた回路の深さと最終的なエネルギー推定値にどのような影響を与えるのでしょうか？

参照

[1] スポール、J. C. (2002). 確率的最適化のための同時摂動アルゴリズムの実装。 『IEEE Transactions on Aerospace and Electronic Systems』, 34(3), 817-823.

[2] Sahin, M. Emre, 他 (2025). Qiskit Machine Learning ：量子ハードウェアおよび古典シミュレータ上で大規模な量子機械学習タスクを実行するためのオープンソースライブラリ。 arXiv:2505.17756.

次のステップ