Qiskit Runtimeですべてを統合する

まとめ

Victoria Lipinskaが、これまでに学んだ内容の最終的なまとめを行います。

参考文献

上記の動画で参照されている論文です。

• Quantum Chemistry in the Age of Quantum Computing, Cao, et al.
• Quantum computational chemistry, McArdle, et al.

Qiskitパターンを使ったVQE

VQE計算に必要なすべての要素が揃いました。

• ハミルトニアン
• アンサッツ
• 古典的最適化器

あとはこれらをQiskitパターンのフレームワーク上で統合するだけです。

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

前述のとおり、ここでは適切な形式に整えられた対象ハミルトニアンがすでに生成済みであることを前提とします。その方法について疑問がある場合は、ハミルトニアンに関するレッスンを参照してください。以下のコードブロックでは、前のレッスンで説明したコンポーネントをセットアップします。ここではH2をモデルとして選択しています。そのハミルトニアンは書き出すのに十分なほどコンパクトだからです。

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-aer qiskit-ibm-runtime scipy

# General imports
import numpy as np
from qiskit.quantum_info import SparsePauliOp

# Hamiltonian obtained from a previous lesson

H = SparsePauliOp(
    [
        "IIII",
        "IIIZ",
        "IZII",
        "IIZI",
        "ZIII",
        "IZIZ",
        "IIZZ",
        "ZIIZ",
        "IZZI",
        "ZZII",
        "ZIZI",
        "YYYY",
        "XXYY",
        "YYXX",
        "XXXX",
    ],
    coeffs=[
        -0.09820182 + 0.0j,
        -0.1740751 + 0.0j,
        -0.1740751 + 0.0j,
        0.2242933 + 0.0j,
        0.2242933 + 0.0j,
        0.16891402 + 0.0j,
        0.1210099 + 0.0j,
        0.16631441 + 0.0j,
        0.16631441 + 0.0j,
        0.1210099 + 0.0j,
        0.17504456 + 0.0j,
        0.04530451 + 0.0j,
        0.04530451 + 0.0j,
        0.04530451 + 0.0j,
        0.04530451 + 0.0j,
    ],
)

nuclear_repulsion = 0.7199689944489797

最初にefficient_su2回路と最適化器COBYLAを選択します。

# Pre-defined ansatz circuit
from qiskit.circuit.library import efficient_su2

# SciPy minimizer routine
from scipy.optimize import minimize

# Plotting functions

# Random initial state and efficient_su2 ansatz
ansatz = efficient_su2(H.num_qubits, su2_gates=["rx"], entanglement="linear", reps=1)
x0 = 2 * np.pi * np.random.random(ansatz.num_parameters)
print(ansatz.decompose().depth())
ansatz.decompose().draw("mpl")

5

次にコスト関数を構築します。これはハミルトニアンと密接に関係していますが、ハミルトニアンが演算子であるのに対し、コスト関数はEstimatorを使ってその演算子の期待値を返す関数である点が異なります。もちろん、アンサッツと変分パラメーターを使って期待値を計算するため、それらはすべて引数として現れます。以下では、実機やシミュレーター向けにそれぞれ少し異なるバージョンを定義しています。

def cost_func(params, ansatz, H, estimator):
    pub = (ansatz, [H], [params])
    result = estimator.run(pubs=[pub]).result()
    energy = result[0].data.evs[0]
    return energy

# def cost_func_sim(params, ansatz, H, estimator):
#    energy = estimator.run(ansatz, H, parameter_values=params).result().values[0]
#    return energy

ステップ2：量子実行に向けて問題を最適化する

使用するハードウェア上でコードをできる限り効率よく動作させる必要があります。そのため、最適化ステップを開始するバックエンドを選択する必要があります。以下のコードは、利用可能なバックエンドの中で最も待ち行列の少ないものを選択します。

# To run on hardware, select the backend with the fewest number of jobs in the queue
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService(channel="ibm_quantum_platform")
backend = service.least_busy(operational=True, simulator=False)
backend.name

実際のバックエンドで実行するための回路の最適化は、重要かつ奥の深いテーマです。ただし、これはVQEに限った話ではありません。ここでは2つの重要な用語を改めて確認しておきます。

• optimization_level：選択したバックエンドのレイアウトに対して、回路がどの程度最適化されているかを表します。最低の最適化レベルでは、回路をデバイス上で動作させるために最低限必要な処理のみを行います。具体的には、回路の量子ビットをデバイスの量子ビットにマッピングし、2量子ビット演算を可能にするスワップゲートを追加します。最高の最適化レベルははるかに高度で、全体のゲート数を削減するためのさまざまな工夫が施されています。多量子ビットゲートはエラー率が高く、量子ビットは時間とともにデコヒーレンスを起こすため、回路が短いほど良い結果が得られます。
• Dynamical Decoupling（動的デカップリング）：アイドル状態の量子ビットに一連のゲートを適用できます。これにより、環境との不要な相互作用の一部をキャンセルできます。 回路の最適化についての詳細は、リンク先のドキュメントをご参照ください。以下のコードは、qiskit.transpilerのプリセットパスマネージャーを使ってパスマネージャーを生成します。

from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
    ConstrainedReschedule,
)
from qiskit.circuit.library import XGate

target = backend.target
pm = generate_preset_pass_manager(target=target, optimization_level=3)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(target=target),
        ConstrainedReschedule(
            acquire_alignment=target.acquire_alignment,
            pulse_alignment=target.pulse_alignment,
            target=target,
        ),
        PadDynamicalDecoupling(
            target=target,
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)

# Use the pass manager and draw the resulting circuit
ansatz_isa = pm.run(ansatz)
ansatz_isa.draw(output="mpl", idle_wires=False, style="iqp")

同様に、ハミルトニアンにもデバイスのレイアウト特性を適用する必要があります。

hamiltonian_isa = H.apply_layout(ansatz_isa.layout)
hamiltonian_isa

SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXXXII'],
              coeffs=[-0.09820182+0.j, -0.1740751 +0.j, -0.1740751 +0.j,  0.2242933 +0.j,
  0.2242933 +0.j,  0.16891402+0.j,  0.1210099 +0.j,  0.16631441+0.j,
  0.16631441+0.j,  0.1210099 +0.j,  0.17504456+0.j,  0.04530451+0.j,
  0.04530451+0.j,  0.04530451+0.j,  0.04530451+0.j])

ステップ3：Qiskitプリミティブを使って実行する

選択したハードウェア上で実行する前に、シミュレーターを使って大まかなデバッグを行い、場合によってはエラーの推定も行うと良いでしょう。そのため、ここではシミュレーター上でVQEを実行する方法を簡単に紹介します。ただし、古典コンピューター、シミュレーター、GPUのいずれも、高度にエンタングルした127量子ビットの量子コンピューターの全機能を正確にシミュレートすることはできない点に注意が必要です。現在の量子ユーティリティの時代においては、シミュレーターの用途は限られています。

変分回路のパラメーター選択ごとに期待値を計算する必要があること（それが最小化する値だから）はすでにお気づきかと思います。最も効率的な方法は、Qiskitプリミティブ「Estimator」を使うことです。まずはローカルシミュレーターを使い始めます。そのためにはBackendEstimatorというEstimatorのローカル版が必要です。

最適化に使用した実バックエンドを保持したまま、そのデバイスのノイズ特性のモデルをインポートし、選択したローカルシミュレーターと組み合わせて使用できます。ここではaer_simulator_statevectorを使用します。

# We will start by using a local simulator
from qiskit_aer import AerSimulator

# Import an estimator, this time from qiskit (we will import from Runtime for real hardware)
from qiskit.primitives import BackendEstimatorV2

# generate a simulator that mimics the real quantum system
backend_sim = AerSimulator.from_backend(backend)
estimator = BackendEstimatorV2(backend=backend_sim)

いよいよVQEを実装します。選択したハミルトニアン、アンサッツ、古典的最適化器、そして実バックエンドに基づくBackendEstimatorを使ってコスト関数を最小化します。ここでは最大反復回数をかなり小さく設定していますが、これはシミュレーターをデバッグ目的で使っているためです。VQEの最適化ステップは、収束するまでに数百回の反復を要することがよくあります。

res = minimize(
    cost_func,
    x0,
    args=(ansatz_isa, hamiltonian_isa, estimator),
    method="cobyla",
    options={"maxiter": 10, "disp": True},
)

print(getattr(res, "fun") - nuclear_repulsion)
print(res)

Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 10   Least value of F = -0.11556938907226563
The corresponding X is:
[4.11796514 4.52126324 0.69570423 4.12781503 6.55507846 1.80713073
 0.9645473  6.23812214]

-0.8355383835212453
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.11556938907226563
       x: [ 4.118e+00  4.521e+00  6.957e-01  4.128e+00  6.555e+00
            1.807e+00  9.645e-01  6.238e+00]
    nfev: 10
   maxcv: 0.0

このコードは正しく実行されましたが、期待通り収束はしませんでした。次に実機で計算を実行し、その出力について説明します。実際のバックエンドには、Qiskit Runtime Estimatorを使用します。Qiskit Runtimeセッション内で実行し、そのセッションのオプションを指定するのが一般的です。

from qiskit_ibm_runtime import QiskitRuntimeService, Session
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime.options import EstimatorOptions

セッションを使用することで、ジョブがキューで待機するのは最初の1回だけになります。古典的最適化器のその後の反復処理では、キューへの再投入は行われません。セッション内では、レジリエンスレベルと最適化レベルを設定できます。これらのツールは非常に重要なため、VQEにおける各ツールの役割と重要性を簡単に振り返り、詳細を学ぶためのリンクを紹介します。

• Runtimeセッション：VQEは本質的に反復的なアルゴリズムで、古典的最適化器が新しい変分パラメーターを選択するたびに、次の試行で新しいゲートが使用されます。セッションを使わない場合、各試行回路の間に追加のキュー待機が発生する可能性があります。VQE計算をセッション内にまとめることで、ジョブ開始前の初回キュー待機のみで済み、変分ステップの間の追加待機はなくなります。この戦略は前のレッスンの例でも使用されていましたが、ジオメトリを変化させる場合にはさらに重要な役割を果たします。セッションの詳細については、実行モードのドキュメントをご参照ください。
• Estimatorの組み込み最適化：Estimatorには計算を最適化するための組み込みオプションがあります。多くのコンテキスト（Estimatorを含む）では、設定値は0と1に制限されており、0は最適化なし、1（デフォルト）は選択したハードウェアへの一定の回路最適化を意味します。他のコンテキストでは0、1、2、3の設定が可能な場合もあります。各設定で使用される具体的な手法については、ドキュメントをご参照ください。ここでは最適化を0に設定し、skip_transpilation = trueを使用します。これは、上記の最適化セクションでパスマネージャーを使ってすでに回路をトランスパイルしているためです。
• Estimatorの組み込みレジリエンス：最適化と同様に、Estimatorにはエラーに対するレジリエンスの組み込み設定があり、異なるエラー軽減手法に対応しています。レジリエンスレベルの設定については、ドキュメントをご参照ください。

エラー軽減はVQE計算の収束において微妙な役割を果たすことに注目する価値があります。古典的最適化器はエネルギーを最小化するパラメーターを探して、パラメーター空間を探索します。最適パラメーターから大きく外れているときは、エラーがあっても急勾配が古典的最適化器に見えるかもしれません。しかし計算が収束し最適値に近づくにつれて、勾配は小さくなり、エラーによって容易に打ち消されてしまいます。どの程度のエラー軽減を使うべきか、また収束のどのタイミングで適用するか、これらは特定のユースケースに応じて自分で判断する必要があります。

この最初の実機実行では、比較的短時間での実行を優先するため、レジリエンスを0に設定しています。本格的なアプリケーションではエラー軽減を使用することをお勧めします。以下のセルには2種類のオプションがあります。(1) Runtimeセッション用のオプション（ここでは「session_options」と命名）と、(2) 古典的最適化器用のオプション（ここでは単に「options」と呼ぶ）です。

estimator_options = EstimatorOptions(resilience_level=0, default_shots=2000)
with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    res = minimize(
        cost_func,
        x0,
        args=(ansatz_isa, hamiltonian_isa, estimator),
        method="cobyla",
        options={"maxiter": 10, "disp": True},
    )

Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 10   Least value of F = -0.11691688904
The corresponding X is:
[5.11796514 5.52126324 0.69570423 5.12781503 6.55507846 1.80713073
 1.9645473  6.23812214]

ジョブの進捗状況は、IBM Quantum® Platformのワークロードから確認できます。

print(getattr(res, "fun") - nuclear_repulsion)
print(res)

-0.8368858834889796
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.11691688904
       x: [ 5.118e+00  5.521e+00  6.957e-01  5.128e+00  6.555e+00
            1.807e+00  1.965e+00  6.238e+00]
    nfev: 10
   maxcv: 0.0

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

少し立ち止まって、これらの出力の意味を確認しましょう。「fun」の出力は、コスト関数に対して得られた最小値です（最後に計算された値とは限りません）。これは核反発を含む全エネルギーであるため、electron_energyも別に定義しています。

上記のケースでは、関数評価の最大回数に達したというメッセージと、関数評価回数（nfev）が10であることが示されています。これは単に、最適化の収束に関するその他の基準が満たされなかったことを意味し、言い換えると、基底状態エネルギーが見つかったと考える理由はありません。successが「False」であることも同じ意味です。

最後にxがあります。これは変分パラメーターのベクトルです。最小コスト関数（エネルギー期待値）をもたらした計算に使用されたパラメーターです。これら8つの値は、アンサッツ中の可変回転角を持つゲートの8つの回転角に対応しています。

おめでとうございます！IBM Quantum QPU上でVQE計算を実行できました！

次のレッスンでは、ハミルトニアンに変数を含むようにこのワークフローを調整する方法を見ていきます。量子化学の問題においては、分子の形状や結合部位を決定するためにジオメトリを変化させることなどが考えられます。

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

2.1.0
0.40.1