期待値推定のための演算子逆伝播（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$ は係数です。回路を単一のユニタリ $U$ として表し、下図に示すように $U = U_C U_Q$ と論理的に分割できるものとします。

演算子逆伝播は、ユニタリ $U_C$ を $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$ はスカラーです。以下が成り立つことは容易に確認できます。

$$U_s^{\dagger} P U_s = P \qquad \text{if} ~[P,P_s] = 0,$$ $$U_s^{\dagger} P U_s = \qquad cos(\theta_s)P + i sin(\theta_s)P_sP \qquad \text{if} ~\{P,P_s\} = 0$$

上記の例において、$\{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）

セットアップ

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

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)]

# Instantiate a CouplingMap object
coupling_map = CouplingMap(layout)
graphviz_draw(coupling_map.graph, method="circo")

次に、ハイゼンベルグXYZハミルトニアンをモデル化するパウリ演算子を生成します。

$${\hat{\mathcal{H}}_{XYZ} = \sum_{(j,k)\in E} (J_{x} \sigma_j^{x} \sigma_{k}^{x} + J_{y} \sigma_j^{y} \sigma_{k}^{y} + J_{z} \sigma_j^{z} \sigma_{k}^{z}) + \sum_{j\in V} (h_{x} \sigma_j^{x} + h_{y} \sigma_j^{y} + h_{z} \sigma_j^{z}),}$$

ここで $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}$ を使用しています。

# Get a qubit operator describing the Heisenberg XYZ model
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を用いて時間発展回路を構築します。

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"Separated the circuit into {len(slices)} slices.")

Separated the circuit into 18 slices.

逆伝播中の演算子の大きさの制約

逆伝播中、演算子の項数は一般的に急速に $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'$ の項は一切切り捨てられていません。逆伝播された回路と逆伝播された演算子は、元の回路と演算子と正確に同じ結果を与えます。

# Backpropagate slices onto the observable
bp_obs, remaining_slices, metadata = backpropagate(
    observable, slices, operator_budget=op_budget
)
# Recombine the slices remaining after backpropagation
bp_circuit = combine_slices(remaining_slices)

print(f"Backpropagated {metadata.num_backpropagated_slices} slices.")
print(
    f"New observable has {len(bp_obs.paulis)} terms, which can be combined into "
    f"{len(bp_obs.group_commuting(qubit_wise=True))} groups."
)
print(
    f"Note that backpropagating one more slice would result in "
    f"{metadata.backpropagation_history[-1].num_paulis[0]} terms "
    f"across {metadata.backpropagation_history[-1].num_qwc_groups} groups."
)
print("The remaining circuit after backpropagation looks as follows:")
bp_circuit.draw("mpl", fold=-1, scale=0.6)

Backpropagated 6 slices.
New observable has 60 terms, which can be combined into 6 groups.
Note that backpropagating one more slice would result in 114 terms across 12 groups.
The remaining circuit after backpropagation looks as follows:

小規模のシミュレータの例では、切り捨ては使用しません。ノイズがない場合、逆伝播ありと逆伝播なしの回路は同じ結果をもたらすため、切り捨てにより近似が加わると結果が悪化するためです。

回路を基底ゲートセットにトランスパイルする

ここでは、元の回路と逆伝播された回路の両方をバックエンドの基底ゲートにトランスパイルします。小規模な例ではシミュレータで実行するため、実際のバックエンドでトランスパイルする必要はありません。

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, basis_gates=backend.configuration().basis_gates
)
isa_circuit = pm_basis.run(circuit)
isa_bp_circuit = pm_basis.run(bp_circuit)

ステップ3: Qiskitプリミティブを使用した実行

まず、元の回路と逆伝播された回路に対応する2つのPrimitive Unified Blocs（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)]

# Instantiate a CouplingMap object
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),
)

# Generate a time evolution circuit for the Hamiltonian
circuit = generate_time_evolution_circuit(
    hamiltonian,
    time=0.2,
    synthesis=LieTrotter(reps=4),
)

# Define the observable to measure
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)

# Define the maximum number of qwc groups allowed in the
# backpropagated observable,
# and the truncation error budget
op_budget = OperatorBudget(max_qwc_groups=15)
truncation_error_budget = setup_budget(
    max_error_total=0.03, max_error_per_slice=0.005
)

# First backpropagation without truncation
bp_obs, remaining_slices, metadata = backpropagate(
    observable, slices, operator_budget=op_budget
)
bp_circuit = combine_slices(remaining_slices)

# Now backpropagate with truncation, using the same operator budget and
# the defined truncation error budget
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
)

# Now we transpile the original circuit and the two backpropagated circuits,
# and apply the layout to the corresponding observables
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
)

# Compare the 2-qubit depth of each transpiled circuit to see how much
# depth backpropagation saved
print(
    f"2-qubit depth without backpropagation: "
    f"{isa_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
    f"2-qubit depth with backpropagation: "
    f"{isa_bp_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
    f"2-qubit depth with backpropagation and truncation: "
    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),
]

# Now we instantiate the Estimator primitive for the hardware with
# ZNE and measurement error
# mitigation and compute the three circuits and observables
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)

# Retrieve the results and the standard deviations
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-qubit depth without backpropagation: 24
2-qubit depth with backpropagation: 20
2-qubit depth with backpropagation and truncation: 18

print(f"Expectation value without backpropagation: {result_no_bp}")
print(f"Backpropagated expectation value: {result_bp}")
print(f"Backpropagated expectation value with truncation: {result_bp_trunc}")

Expectation value without backpropagation: 0.9543907942381811
Backpropagated expectation value: 0.9445337385406468
Backpropagated expectation value with truncation: 0.934050286970965

# Plot the results
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: OBP、古典的な期待値推定、ノイズシミュレーションをカバーするチュートリアル付きのRustによる高速化パウリ伝播パッケージ