期待値推定のための演算子逆伝播（OBP）#
使用量の目安: Heron r3 プロセッサで 4 分（注意: これはあくまで目安です。実際の実行時間は異なる場合があります。）
学習成果#
このチュートリアルを終えると、次の内容が理解できるようになります。
qiskit-addon-obpを使って、回路実行回数の増加と引き換えに量子回路の深さを削減する方法qiskit-addon-utilsを使って XYZ ハミルトニアンとその時間発展回路を構成する方法
前提知識#
このチュートリアルに進む前に、次のトピックについて理解しておくことをお勧めします。
観測量の期待値を計算するための Estimator プリミティブの使い方
背景#
演算子逆伝播（operator backpropagation）は、量子回路の末尾にある演算を測定対象の観測量へ吸収させる手法であり、観測量の項数が増える代わりに回路の深さを一般に削減します。目標は、観測量が大きくなりすぎない範囲で、回路をできるだけ多く逆伝播させることです。Qiskit ベースの実装は OBP Qiskit アドオンとして提供されています。詳しくは対応する ドキュメント を参照してください。
観測量 \(O = \sum_P c_P P\)（ここで \(P\) はパウリ演算子、\(c_P\) は係数）を測定したい例題の回路を考えます。この回路を 1 つのユニタリー \(U\) として表すと、下図のように \(U = U_C U_Q\) と論理的に分割できます。
演算子逆伝播は、\(O' = U_C^{\dagger}OU_C = \sum_P c_P U_C^{\dagger}PU_C\) として観測量を発展させることで、ユニタリー \(U_C\) を観測量へ吸収します。言い換えると、計算の一部が、観測量を \(O\) から \(O'\) へ発展させることによって古典的に実行されます。これにより元の問題は、ユニタリーが \(U_Q\) である、より深さの小さい新しい回路について観測量 \(O'\) を測定する問題として再定式化できます。
ユニタリー \(U_C\) は複数のスライスとして \(U_C = U_S U_{S-1}...U_2U_1\) と表されます。スライスの定義方法は複数あります。たとえば上の例題回路では、\(R_{zz}\) の各層と \(R_x\) ゲートの各層をそれぞれ個別のスライスと見なせます。逆伝播では \(O' = \Pi_{s=1}^S \sum_P c_P U_s^{\dagger} P U_s\) を古典的に計算します。各スライス \(U_s\) は \(U_s = exp(\frac{-i\theta_s P_s}{2})\) と表せます。ここで \(P_s\) は \(n\) 量子ビットのパウリ演算子、\(\theta_s\) はスカラーです。次のことは容易に確認できます。
上の例で \(\{P,P_s\} = 0\) の場合、期待値を計算するには 1 つではなく 2 つの量子回路を実行する必要があります。したがって逆伝播は観測量の項数を増やす可能性があり、その結果として回路の実行回数も増えます。演算子が大きくなりすぎるのを防ぎながら回路のより深くまで逆伝播させる方法の 1 つは、係数の小さい項を演算子に追加せずに切り捨てることです。たとえば上の例では、\(\theta_s\) が十分小さければ \(P_sP\) を含む項を切り捨てるという選択ができます。項を切り捨てると実行すべき量子回路の数は減りますが、その代わり、切り捨てた項の係数の大きさに比例した誤差が最終的な期待値の計算に生じます。
要件#
このチュートリアルを始める前に、次のものがインストールされていることを確認してください。
Qiskit SDK v2.0 以降（可視化 サポート付き）
Qiskit Runtime v0.22 以降（
pip install qiskit-ibm-runtime）OBP Qiskit アドオン 0.3 以降（
pip install qiskit-addon-obp）Qiskit アドオンユーティリティ 0.3 以降（
pip install qiskit-addon-utils）
セットアップ#
import numpy as np
import matplotlib.pyplot as plt
from qiskit.primitives import StatevectorEstimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler import CouplingMap
from qiskit.synthesis import LieTrotter
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
from qiskit_addon_utils.problem_generators import (
generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth, combine_slices
from qiskit_addon_obp.utils.simplify import OperatorBudget
from qiskit_addon_obp import backpropagate
from qiskit_addon_obp.utils.truncating import setup_budget
from rustworkx.visualization import graphviz_draw
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import EstimatorV2, EstimatorOptions
小規模なシミュレーターでの例#
このチュートリアルでは、OBP Qiskit アドオン を使ってハイゼンベルクスピン鎖の量子ダイナミクスをシミュレートする Qiskit パターン を実装します。なお、ノイズのないシミュレーターでは、逆伝播を行った場合と行わない場合で得られる期待値は同じになります。
ステップ 1: 古典的な入力を量子問題へマッピングする#
量子ハイゼンベルクモデルの時間発展を量子実験へマッピングする#
まず、qiskit-addon-utils の generate_xyz_hamiltonian 関数を使い、与えられた接続グラフ上でハイゼンベルク型のハミルトニアンを生成します。このグラフは rustworkx.PyGraph または CouplingMap のいずれでもかまいません。以下では、10 量子ビットの直線鎖の CouplingMap を使います。
num_qubits = 10
layout = [(i - 1, i) for i in range(1, num_qubits)]
# CouplingMap オブジェクトを生成する
coupling_map = CouplingMap(layout)
graphviz_draw(coupling_map.graph, method="circo")
次に、ハイゼンベルク XYZ ハミルトニアンを表すパウリ演算子を生成します。
ここで \(G(V,E)\) はカップリングマップのグラフです。このチュートリアルでは \(J_x, J_y, J_z\) をそれぞれ \(\frac{\pi}{8}, \frac{\pi}{4}, \frac{\pi}{2}\)、\(h_x, h_y, h_z\) をそれぞれ \(\frac{\pi}{3}, \frac{\pi}{6}, \frac{\pi}{9}\) としています。
# ハイゼンベルク XYZ モデルを表す量子ビット演算子を取得する
hamiltonian = generate_xyz_hamiltonian(
coupling_map,
coupling_constants=(np.pi / 8, np.pi / 4, np.pi / 2),
ext_magnetic_field=(np.pi / 3, np.pi / 6, np.pi / 9),
)
print(hamiltonian)
SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIIY', 'IIIIIIIIIZ', 'IIIIIIIIXI', 'IIIIIIIIYI', 'IIIIIIIIZI', 'IIIIIIIXII', 'IIIIIIIYII', 'IIIIIIIZII', 'IIIIIIXIII', 'IIIIIIYIII', 'IIIIIIZIII', 'IIIIIXIIII', 'IIIIIYIIII', 'IIIIIZIIII', 'IIIIXIIIII', 'IIIIYIIIII', 'IIIIZIIIII', 'IIIXIIIIII', 'IIIYIIIIII', 'IIIZIIIIII', 'IIXIIIIIII', 'IIYIIIIIII', 'IIZIIIIIII', 'IXIIIIIIII', 'IYIIIIIIII', 'IZIIIIIIII', 'XIIIIIIIII', 'YIIIIIIIII', 'ZIIIIIIIII'],
coeffs=[0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j,
0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j,
1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j,
0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j,
0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j,
1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j,
0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 1.04719755+0.j,
0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j,
1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j,
0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j,
1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j,
0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
0.34906585+0.j])
この量子ビット演算子から、その時間発展を表す量子回路を生成できます。ここでは generate_time_evolution_circuit を Lie-Trotter 分解とともに使い、時間発展回路を構成しています。
circuit = generate_time_evolution_circuit(
hamiltonian,
time=0.2,
synthesis=LieTrotter(reps=2),
)
circuit.draw("mpl", style="iqp", fold=-1)
ステップ 2: 量子ハードウェア実行に向けて問題を最適化する#
逆伝播する回路スライスを作成する#
backpropagate 関数は、一度に回路スライス全体を逆伝播します。したがってスライスの分け方の選択は、与えられた問題に対して逆伝播がどれだけうまく機能するかに影響し得ます。ここでは slice_by_depth 関数を使い、同じ種類のゲートをスライスにまとめます。
回路のスライス分割についてより詳しくは、qiskit-addon-utils パッケージの ハウツーガイド をご覧ください。
slices = slice_by_depth(circuit, max_slice_depth=1)
print(f"回路を {len(slices)} 個のスライスに分割しました。")
回路を 18 個のスライスに分割しました。
逆伝播中に演算子が大きくなり得る範囲を制限する#
逆伝播の途中で、演算子の項数は一般に急速に \(2^L\) に近づきます（\(L\) はスライス数）。演算子の 2 つの項が量子ビットごとに可換でない場合、それらに対応する期待値を得るには別々の回路が必要になります。たとえば 2 量子ビットの観測量 \(O = 0.1 XX + 0.3 IZ - 0.5 IX\) を考えると、\([XX,IX] = 0\) なので、これら 2 つの項の期待値は単一の基底での測定で計算できます。しかし \(IZ\) は他の 2 つの項と反可換なので、\(IZ\) の期待値を計算するには別の基底での測定が必要です。言い換えると、\(\langle O \rangle\) を計算するには 1 つではなく 2 つの回路が必要になります。演算子の項数が増えるにつれて、必要な回路実行回数も増える可能性があります。
演算子のサイズは、backpropagate 関数の operator_budget キーワード引数を指定することで制限できます。この引数は OperatorBudget のインスタンスを受け取ります。
割り当てる追加リソースの量（回路実行回数、ひいては必要となる QPU 時間）を制御するため、逆伝播後の観測量が持てる、量子ビットごとに可換なパウリ群の最大数を制限します。ここでは、演算子内の量子ビットごとに可換なパウリ群の数が 8 を超えた時点で逆伝播を停止するよう指定します。
op_budget = OperatorBudget(max_qwc_groups=8)
回路からスライスを逆伝播する#
まず観測量を \(M_Z = \frac{1}{N} \sum_{i=1}^N \langle Z_i \rangle\)（\(N\) は量子ビット数）と指定します。観測量の項を、量子ビットごとに可換な 8 個以下のパウリ群にまとめられなくなるまで、時間発展回路からスライスを逆伝播します。
observable = SparsePauliOp.from_sparse_list(
[("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
num_qubits=num_qubits,
)
observable
SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
coeffs=[0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j,
0.1+0.j, 0.1+0.j])
以下では 6 つのスライスを逆伝播し、項が 8 ではなく 6 個の群にまとめられたことがわかります。これは、もう 1 つスライスを逆伝播するとパウリ群の数が 8 を超えてしまうことを意味します。返されたメタデータを調べれば、実際にそうであることを確認できます。また、この部分の回路変換は厳密であることにも注意してください。つまり、新しい観測量 \(O'\) の項は 1 つも切り捨てられていません。逆伝播後の回路と逆伝播後の演算子は、元の回路と演算子とまったく同じ結果を与えます。
# スライスを観測量へ逆伝播する
bp_obs, remaining_slices, metadata = backpropagate(
observable, slices, operator_budget=op_budget
)
# 逆伝播後に残ったスライスを再結合する
bp_circuit = combine_slices(remaining_slices)
print(f"{metadata.num_backpropagated_slices} 個のスライスを逆伝播しました。")
print(
f"新しい観測量は {len(bp_obs.paulis)} 個の項を持ち、"
f"{len(bp_obs.group_commuting(qubit_wise=True))} 個の群にまとめられます。"
)
print(
f"なお、もう 1 つスライスを逆伝播すると "
f"{metadata.backpropagation_history[-1].num_paulis[0]} 個の項が "
f"{metadata.backpropagation_history[-1].num_qwc_groups} 個の群にわたって生じます。"
)
print("逆伝播後に残った回路は次のようになります:")
bp_circuit.draw("mpl", fold=-1, scale=0.6)
6 個のスライスを逆伝播しました。
新しい観測量は 60 個の項を持ち、6 個の群にまとめられます。
なお、もう 1 つスライスを逆伝播すると 114 個の項が 12 個の群にわたって生じます。
逆伝播後に残った回路は次のようになります:
シミュレーター上の小規模な例では、切り捨てを使いません。これは、ノイズがない場合には逆伝播の有無にかかわらず同じ結果になり、切り捨ては近似が加わることで結果を悪化させるだけだからです。
回路を基底ゲート集合へトランスパイルする#
ここでは、元の回路と逆伝播後の回路の両方を、バックエンドの基底ゲートへトランスパイルします。小規模なインスタンスではシミュレーター上で実行するため、実際のバックエンドに合わせてトランスパイルする必要はありません。
service = QiskitRuntimeService()
backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=133
)
print(backend)
<IBMBackend('ibm_kingston')>
pm_basis = generate_preset_pass_manager(
optimization_level=3
)
isa_circuit = pm_basis.run(circuit)
isa_bp_circuit = pm_basis.run(bp_circuit)
ステップ 3: Qiskit プリミティブを使って実行する#
まず、元の回路と逆伝播後の回路に対応する 2 つの Primitive Unified Bloc（PUB）を作成します。次に、期待値を得るために理想的な Estimator でこれらの PUB を実行します。
pubs = [(isa_circuit, observable), (isa_bp_circuit, bp_obs)]
rng = np.random.default_rng()
estimator = StatevectorEstimator(seed=rng)
job = estimator.run(pubs)
ステップ 4: 後処理を行い、望ましい古典的形式で結果を返す#
ここでは、元の回路と逆伝播後の回路の期待値を取得します。
primitive_result = job.result()
circuit_expval = primitive_result[0].data.evs.item()
bp_circuit_expval = primitive_result[1].data.evs.item()
methods = [
"No backpropagation",
"Backpropagation",
]
values = [circuit_expval, bp_circuit_expval]
ax = plt.gca()
plt.bar(methods, values, color="#a56eff", width=0.4, edgecolor="#8a3ffc")
ax.set_ylim([0.6, 0.92])
ax.set_ylabel(r"$M_Z$", fontsize=12)
Text(0, 0.5, '$M_Z$')
予想どおり、2 つの期待値は一致します。ノイズのない状態ベクトルシミュレーター上で実行しているため、逆伝播は回路と観測量の組に対する厳密な変換であり、元のワークフローと逆伝播後のワークフローは同じ \(M_Z\) の値を返すはずです。逆伝播の利点が現れるのはノイズのあるハードウェア上だけです。そこでは、より短い逆伝播後の回路の方が蓄積する誤差が小さくなります。これを以下の大規模なハードウェアの例で示します。
大規模なハードウェアでの例#
実験を開発する際には、可視化やシミュレーションを容易にするために小さな回路から始めるのが有用です。ここでは、\(J\) と \(h\) のパラメーターを同じ値、観測量も同じ \(M_Z\) としたまま、50 量子ビットのハイゼンベルクハミルトニアンについて、4 トロッターステップでの演算子逆伝播を見ていきます。この規模では理想的な期待値を力任せの方法で計算することはできないため、テンソルネットワークを使い、理想的な期待値が \(\simeq 0.89\) であることを求めています。
この大規模な例では、逆伝播に加えて、切り捨てを伴う逆伝播も導入します。実効的な回路の深さを減らすためには、理想的にはできるだけ多く逆伝播したいところです。しかしそうすると、更新後の観測量に非可換な項が多数生じ、量子側のオーバーヘッドが増えることがよくあります。そこで、切り捨てと呼ばれる手法を用いて、係数の小さい観測量の項を除去できます。切り捨ては更新後の観測量の項数を減らすことでより多くの逆伝播を可能にしますが、その一方で近似も持ち込みます。したがって、近似誤差が、より深い逆伝播によって得られるノイズ低減を上回ってしまわないよう、切り捨てを一定の範囲に制限する必要があります。
切り捨ての量を制限するため、setup_budget 関数を使って、スライスごとの誤差バジェットと、逆伝播した回路全体にわたる総誤差バジェットの両方を割り当てます。これにより、スライスごとにも回路全体としても切り捨てが制御されます。バジェットの割り当て方法については、こちらのガイド も参照してください。
num_qubits = 50
layout = [(i - 1, i) for i in range(1, num_qubits)]
# CouplingMap オブジェクトを生成する
coupling_map = CouplingMap(layout)
hamiltonian = generate_xyz_hamiltonian(
coupling_map,
coupling_constants=(np.pi / 8, np.pi / 4, np.pi / 2),
ext_magnetic_field=(np.pi / 3, np.pi / 6, np.pi / 9),
)
# ハミルトニアンの時間発展回路を生成する
circuit = generate_time_evolution_circuit(
hamiltonian,
time=0.2,
synthesis=LieTrotter(reps=4),
)
# 測定する観測量を定義する
observable = SparsePauliOp.from_sparse_list(
[("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
num_qubits,
)
slices = slice_by_depth(circuit, max_slice_depth=1)
# 逆伝播後の観測量に許容する量子ビットごとに可換な（qwc）群の最大数と、
# 切り捨て誤差のバジェットを定義する
op_budget = OperatorBudget(max_qwc_groups=15)
truncation_error_budget = setup_budget(
max_error_total=0.03, max_error_per_slice=0.005
)
# まず切り捨てなしで逆伝播する
bp_obs, remaining_slices, metadata = backpropagate(
observable, slices, operator_budget=op_budget
)
bp_circuit = combine_slices(remaining_slices)
# 次に、同じ演算子バジェットと定義した切り捨て誤差バジェットを使い、
# 切り捨てを伴う逆伝播を行う
bp_obs_trunc, remaining_slices_trunc, metadata = backpropagate(
observable,
slices,
operator_budget=op_budget,
truncation_error_budget=truncation_error_budget,
)
bp_circuit_trunc = combine_slices(
remaining_slices_trunc, include_barriers=False
)
# 続いて、元の回路と 2 つの逆伝播後の回路をトランスパイルし、
# 対応する観測量にレイアウトを適用する
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)
isa_circuit = pm.run(circuit)
isa_bp_circuit = pm.run(bp_circuit)
isa_bp_circuit_trunc = pm.run(bp_circuit_trunc)
isa_observable = observable.apply_layout(isa_circuit.layout)
isa_bp_observable = bp_obs.apply_layout(isa_bp_circuit.layout)
isa_bp_observable_trunc = bp_obs_trunc.apply_layout(
isa_bp_circuit_trunc.layout
)
# 逆伝播によってどれだけ深さが削減されたかを見るため、
# トランスパイル後の各回路の 2 量子ビット深さを比較する
print(
f"逆伝播なしの 2 量子ビット深さ: "
f"{isa_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
f"逆伝播ありの 2 量子ビット深さ: "
f"{isa_bp_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
f"逆伝播および切り捨てありの 2 量子ビット深さ: "
f"{isa_bp_circuit_trunc.depth(lambda x: x.operation.num_qubits == 2)}"
)
pubs = [
(isa_circuit, isa_observable),
(isa_bp_circuit, isa_bp_observable),
(isa_bp_circuit_trunc, isa_bp_observable_trunc),
]
# 次に、ZNE と測定誤差の低減を有効にしたハードウェア用の Estimator
# プリミティブを生成し、3 つの回路と観測量について計算する
options = EstimatorOptions()
options.default_precision = 0.01
options.resilience_level = 2
options.resilience.zne.noise_factors = [1, 1.2, 1.4]
options.resilience.zne.extrapolator = ["linear"]
estimator = EstimatorV2(mode=backend, options=options)
estimator.options.environment.job_tags = ["TUT_OBP"]
job = estimator.run(pubs)
# 結果と標準偏差を取得する
result_no_bp = job.result()[0].data.evs.item()
result_bp = job.result()[1].data.evs.item()
result_bp_trunc = job.result()[2].data.evs.item()
std_no_bp = job.result()[0].data.stds.item()
std_bp = job.result()[1].data.stds.item()
std_bp_trunc = job.result()[2].data.stds.item()
逆伝播なしの 2 量子ビット深さ: 24
逆伝播ありの 2 量子ビット深さ: 20
逆伝播および切り捨てありの 2 量子ビット深さ: 18
print(f"逆伝播なしの期待値: {result_no_bp}")
print(f"逆伝播後の期待値: {result_bp}")
print(f"切り捨てを伴う逆伝播後の期待値: {result_bp_trunc}")
逆伝播なしの期待値: 0.9337875596354408
逆伝播後の期待値: 0.9107280913377572
切り捨てを伴う逆伝播後の期待値: 0.9198644768080684
# 結果をプロットする
methods = [
"No backpropagation",
"Backpropagation",
"Backpropagation w/ truncation",
]
values = [result_no_bp, result_bp, result_bp_trunc]
error_bars = [std_no_bp, std_bp, std_bp_trunc]
ax = plt.gca()
plt.bar(methods, values, color="#a56eff", width=0.4, edgecolor="#8a3ffc")
plt.errorbar(methods, values, yerr=error_bars, fmt="o", color="r", capsize=5)
plt.axhline(0.89)
ax.set_ylim([0.8, 0.98])
plt.text(0.25, 0.895, "Exact result")
ax.set_ylabel(r"$M_Z$", fontsize=12)
Text(0, 0.5, '$M_Z$')
次のステップ#
この内容に興味を持たれた方には、次の資料もお勧めします。
おすすめ
pauli-prop— Rust で高速化されたパウリ伝播のパッケージ。OBP、古典的な期待値推定、ノイズを含むシミュレーションを扱うチュートリアルが付属しています
© IBM Corp., 2017-2026