メインコンテンツへスキップ

ユーティリティスケール実験 I

備考

Tamiya Onodera（2024年7月5日）

元の講義の PDFをダウンロード できます。コードスニペットは静的な画像であるため、一部が非推奨になっている場合があります。

この実験を実行するおおよそのQPU時間は45秒です。

1. ユーティリティ論文の概要​

このレッスンでは、2023年6月15日にNature Vol 618で発表された、通称「ユーティリティ論文」に登場するユーティリティスケールの回路を実行します。この論文は、2次元横磁場イジングモデルの時間発展を扱っています。特に、以下のハミルトニアンの時間ダイナミクスを考察します：

H=HZZ+HX=−J∑(i,j)ZiZj+h∑iXiH = H_{ZZ} + H_X = - J \sum_{(i,j)} Z_i Z_j + h \sum_{i} X_i

ここで J>0J > 0 は i<ji < j を満たす最近接スピン間の結合定数、hh はグローバル横磁場です。初期状態からのスピンダイナミクスを、時間発展演算子の一次トロッター分解を用いてシミュレートします：

exp⁡(−iHZZδt)=∏(i,j)exp⁡(iJδtZiZj)=∏(i,j)RZiZj(−2Jδt)exp⁡(−iHXδt)=∏iexp⁡(−ihδtXi)=∏iRXi(2hδt)\begin{aligned} \exp(-i H_{ZZ} \delta t) &= \prod_{(i,j)} \exp (i J \delta t Z_i Z_j) = \prod_{(i,j)} \mathrm{R}_{Z_i Z_j} ( - 2 J \delta t) \\ \exp(-i H_X \delta t) &= \prod_{i} \exp (-i h \delta t X_i ) = \prod_{i} \mathrm{R}_{X_i} ( 2 h \delta t) \end{aligned}

ここで、発展時間 TT は T/δtT / \delta t 個のトロッターステップに離散化され、RZiZj(θJ)\mathrm{R}_{Z_i Z_j}(\theta_J) と RXi(θh)\mathrm{R}_{X_i}(\theta_h) はそれぞれ ZZZZ 回転ゲートと XX 回転ゲートです。

この研究では、127量子ビットのデバイスでヘビーヘックス接続を持つIBM Quantum® Eagle プロセッサ上で実験が行われました。すべての量子ビットに XX 相互作用を適用し、結合マップのすべてのエッジに ZZZZ 相互作用を適用しています。なお、「データ依存性」のため、すべての ZZZZ 相互作用を同時に適用することはできません。そのため、結合マップをカラーリングしてレイヤーにグループ化します。同じレイヤー内のものは同じ色が割り当てられ、並列に適用できます。

また、実験の簡略化のため、θJ=−π/2\theta_J=-\pi /2 の場合に焦点を当てています。

この論文の新しい貢献は、状態ベクトルシミュレーションの限界を超えたスケールで量子回路を構築し、ノイズのある量子コンピュータ上で実行して、信頼性の高い結果を得ることに成功した点です。すなわち、ノイズのある量子コンピュータのユーティリティ（有用性）を実証しました。そのために、確率的誤差増幅（PEA）を用いたゼロノイズ外挿（ZNE）を適用して、ノイズのあるデバイスの誤差を緩和しました。

それ以来、このような実験と回路を「ユーティリティスケール」と呼ぶようになりました。

1.1 目標​

このレッスンの目標は、ユーティリティスケールの回路を構築し、Eagle プロセッサ上で実行することです。PEA は本稿執筆時点ではQiskitの実験的機能であること、また ZNE と PEA を組み合わせるとかなりの時間がかかることから、信頼性の高い結果を得ることはこのノートブックの範囲外とします。

具体的には、論文の Figure 4b に対応する回路を構築・実行し、独自の「未緩和」の点をプロットします。これは ⟨Z62⟩\langle Z_{62} \rangle を観測量とした 127量子ビット × 60レイヤー（20トロッターステップ）の回路です。

image.png

大変そうに見えますか？   心配しないでください。このコースの最後の3つのレッスンがその足がかりを提供します。まず、⟨Z13⟩\langle Z_{13} \rangle を観測量とした 27量子ビット × 6レイヤー（2トロッターステップ）の回路を偽デバイス上で構築・実行する、より小規模な実験を実演します。

以上が概要です。ユーティリティスケールの冒険に出発しましょう！

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

qiskit.__version__
'2.0.2'
#!pip install qiskit_ibm_runtime
#!pip install qiskit_aer
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx

from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit.circuit.library import YGate
from qiskit.quantum_info import SparsePauliOp
from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    fake_provider,
    EstimatorV2 as Estimator,
)
from qiskit_aer import AerSimulator
service = QiskitRuntimeService()

2. 準備​

2.1 RZZ(-π\pi / 2) の構築​

まず、一般に RZZ ゲートは2つの CXCX ゲートを必要とすることを確認します。

from qiskit.circuit.library import RZZGate

θ_h = Parameter("$\\theta_h$")
qc1 = QuantumCircuit(2)
qc1.append(RZZGate(θ_h), [0, 1])
qc1.decompose(reps=1).draw("mpl")

Output of the previous code cell

前述のとおり、この実験では特定の角度 -π\pi / 2 を持つ RZZ ゲートに焦点を当てます。論文に示されているように、これはたった1つの CXCX ゲートで実現できます。

qc2 = QuantumCircuit(2)

qc2.sdg([0, 1])
qc2.append(YGate().power(1 / 2), [1])
qc2.cx(0, 1)
qc2.append(YGate().power(1 / 2).adjoint(), [1])

qc2.draw("mpl")

Output of the previous code cell

後で参照できるよう、この回路をゲートとして定義します。

rzz = qc2.to_gate(label="RZZ")

新しく定義した rzz を試しに使ってみましょう。

qc3 = QuantumCircuit(3)
qc3.append(rzz, [0, 1])
qc3.append(rzz, [0, 2])
display(qc3.draw("mpl"))
# display(qc.decompose(reps=1).draw("mpl"))

Output of the previous code cell

さらに使用する前に、-pi/2 に対する qc1（RZZ ゲート）と新しく定義した rzz または qc2 ゲートが論理的に等価であることを検証しましょう：

from qiskit.quantum_info import Operator

op1 = Operator(qc1.assign_parameters([-np.pi / 2]))
op2 = Operator(qc2)

op1.equiv(op2)
True

2.2 結合マップのカラーリング​

バックエンドの結合マップをどのようにカラーリングするかを学びましょう。これは ZZZZ 相互作用をレイヤーにグループ化するために必要です。

まず、バックエンドの結合マップを可視化しましょう。現在のすべての IBM Quantum デバイスの結合マップはヘビーヘキサゴナルです。

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

backend.coupling_map.draw()

Output of the previous code cell

結合マップをカラーリングするために、グラフや複雑なネットワークを扱うためのPythonパッケージである rustworkx を使用します。このパッケージは複数のカラーリングアルゴリズムを提供していますが、いずれもヒューリスティックであり、最小カラーリングを見つけることは保証されていません。

ただし、ヘビーヘックスグラフは二部グラフであるため、これらのグラフに対して最小カラーリングを見つけられるはずの graph_bipartite_edge_color を使用します。

def color_coupling_map(backend):
    graph = backend.coupling_map.graph
    undirected_graph = graph.to_undirected(multigraph=False)
    edge_color_map = rx.graph_bipartite_edge_color(undirected_graph)
    if edge_color_map is None:
        edge_color_map = rx.graph_greedy_edge_color(undirected_graph)
    # build a map from color to a list of edges
    edge_index_map = undirected_graph.edge_index_map()
    color_edges_map = {color: [] for color in edge_color_map.values()}
    for edge_index, color in edge_color_map.items():
        color_edges_map[color].append(
            (edge_index_map[edge_index][0], edge_index_map[edge_index][1])
        )
    return edge_color_map, color_edges_map

ヘビーヘキサゴナルグラフは3色で塗り分けられるはずです。上記の結合マップでこれを確認しましょう。

edge_color_map, color_edges_map = color_coupling_map(backend)
print(
    f"{backend.name}, {backend.num_qubits}-qubit device, {len(color_edges_map.keys())} colors assigned."
)
ibm_strasbourg, 127-qubit device, 3 colors assigned.

確かにその通りです！

fun として、rustworkx の可視化機能を使って、得られたカラーリングで結合マップを色分けしてみましょう。

color_str_map = {0: "green", 1: "red", 2: "blue"}

undirected_graph = backend.coupling_map.graph.to_undirected(multigraph=False)
for i in undirected_graph.edge_indices():
    undirected_graph.get_edge_data_by_index(i)["color"] = color_str_map[
        edge_color_map[i]
    ]

rx.visualization.graphviz_draw(
    undirected_graph, method="neato", edge_attr_fn=lambda edge: {"color": edge["color"]}
)

Output of the previous code cell

3. 2次元イジングモデルのトロッター化時間発展を解く​

2次元イジングモデルの時間発展に関するユーティリティ論文の回路を構築するルーティンを定義しましょう。このルーティンは、バックエンド、トロッターステップ数を示す整数、バリア挿入を制御する真偽値の3つのパラメータを受け取ります。

def get_utility_circuit(backend, num_steps: int, barrier: bool = False):
    num_qubits = backend.num_qubits
    _, color_edges_map = color_coupling_map(backend)
    θ_h = Parameter("$\\theta_h$")
    qc = QuantumCircuit(num_qubits)

    for i in range(num_steps):
        qc.rx(θ_h, range(num_qubits))

        for _, edge_list in color_edges_map.items():
            for edge in edge_list:
                qc.append(rzz, edge)

        if barrier:
            qc.barrier()
    return qc

構築した回路では、Qubit のマッピングとルーティングをすでに手動で実施していることに注意してください。したがって、後で回路をトランスパイルする際には、Transpiler にQubitのマッピングとルーティングを行わせては（行わせるべきでは）なりません。後ほど示すように、最適化レベルを1、レイアウト方法を "trivial" として呼び出します。

次に、素早い確認のために、構築された回路の情報を取得する簡単なルーティンを定義します。

def get_circuit_info(qc: QuantumCircuit, reps: int = 0):
    qc0 = qc.decompose(reps=reps)
    return (
        f"{qc0.num_qubits} qubits × {qc0.depth(lambda x: x.operation.num_qubits == 2)} layers ({qc0.depth()}-depth)"
        + ",  "
        + f"""Gate breakdown: {", ".join([f"{k.upper()} {v}" for k, v in qc0.count_ops().items()])}"""
    )

これらのルーティンを試してみましょう。27量子ビット × 15レイヤー（5トロッターステップ）の回路が表示されるはずです。偽デバイスには28のエッジがあるため、28×5個のエンタングルメントゲートが存在するはずです。

backend = fake_provider.FakeTorontoV2()
num_steps = 5
qc = get_utility_circuit(backend, num_steps, True)

display(qc.draw(output="mpl", fold=-1))
print(get_circuit_info(qc, reps=0))
print(get_circuit_info(qc, reps=1))

Output of the previous code cell

27 qubits × 15 layers (20-depth),  Gate breakdown: CIRCUIT-165 140, RX 135, BARRIER 5
27 qubits × 15 layers (60-depth),  Gate breakdown: SDG 280, UNITARY 280, CX 140, R 135, BARRIER 5

4. 27量子ビット版の問題を解く​

ここでは、ユーティリティ実験の小規模版を実演します。⟨Z13⟩\langle Z_{13} \rangle を観測量とした 27量子ビット × 6レイヤー（2トロッターステップ）の回路を構築し、AerSimulator と偽デバイスの両方で実行します。

もちろん、Map（マップ）、Optimize（最適化）、Execute（実行）、Post-Process（後処理）からなる4ステップのワークフロー「Qiskit パターン」に従います。具体的には以下のとおりです：

  • 古典的な入力を量子計算にマップする。
  • 量子計算のために回路を最適化する。
  • プリミティブを使用して回路を実行する。
  • 後処理を行い、結果を古典的な形式で返す。

以下では、より小規模な実験の回路を作成するためのMapステップを実施します。次に、AerSimulator用のOptimizeとExecuteのセット、偽デバイス用のOptimizeとExecuteのセットをそれぞれ実施します。最後に、結果をプロットするためのPost-Processステップを実施します。

4.1 ステップ1: Map​

backend = fake_provider.FakeTorontoV2()  # a 27 qubit fake device.
num_steps = 2
qc = get_utility_circuit(backend, num_steps)
obs = SparsePauliOp.from_sparse_list(
    [("Z", [13], 1)], num_qubits=backend.num_qubits
)  # Falcon
angles = [
    0,
    0.1,
    0.2,
    0.3,
    0.4,
    0.5,
    0.6,
    0.7,
    0.8,
    1.0,
    np.pi / 2,
]  # We try 11 angles for theta_h.

4.2 ステップ2と3: Optimize と Execute（シミュレータ）​

backend_sim = AerSimulator()
transpiled_qc_sim = transpile(
    qc, backend_sim, optimization_level=1, layout_method="trivial"
)
transpiled_obs_sim = obs.apply_layout(layout=transpiled_qc_sim.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc_sim, reps=1))
27 qubits × 6 layers (23-depth),  Gate breakdown: SDG 112, UNITARY 112, CX 56, R 54
27 qubits × 6 layers (16-depth),  Gate breakdown: U3 80, CX 56, R 54, U1 32, U 28

あるユーザーが、macOS 14.5 を搭載した 32GB 3LPDDR4X RAM を持つ 2.3 GHz クアッドコア Intel Core i7 プロセッサ搭載の MacBook Pro で次のセルを実行したところ、ウォールタイムで161ミリ秒かかりました。ラップトップによって若干異なります。

%%time
params = [[p] for p in angles]
estimator = Estimator(mode=backend_sim)
pub = (transpiled_qc_sim, transpiled_obs_sim, params)
result_sim = estimator.run([pub]).result()
CPU times: user 231 ms, sys: 186 ms, total: 417 ms
Wall time: 111 ms

4.3 ステップ2と3: Optimize と Execute（偽デバイス）​

backend_fake = fake_provider.FakeTorontoV2()
transpiled_qc_fake = transpile(
    qc, backend_fake, optimization_level=1, layout_method="trivial"
)
transpiled_obs_fake = obs.apply_layout(layout=transpiled_qc_fake.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc_fake, reps=1))
27 qubits × 6 layers (23-depth),  Gate breakdown: SDG 112, UNITARY 112, CX 56, R 54
27 qubits × 6 layers (49-depth),  Gate breakdown: SDG 324, U1 274, H 162, CX 56, U3 14

同じユーザーが上記と同じ環境で次のセルを実行したところ、ウォールタイムで2分19秒かかりました。偽デバイス上で回路を実行するとノイズありシミュレーションが呼び出されるため、厳密なシミュレーションよりもはるかに長い時間がかかります。偽デバイス上では、より大きな回路（27量子ビット × 9レイヤー、3トロッターステップなど）を実行しないことを推奨します。

%%time
params = [[p] for p in angles]
estimator = Estimator(mode=backend_fake)
pub = (transpiled_qc_fake, transpiled_obs_fake, params)
result_fake = estimator.run([pub]).result()
CPU times: user 4min 42s, sys: 9.35 s, total: 4min 51s
Wall time: 38.3 s

4.4 ステップ4: Post-process​

厳密シミュレーションとノイズありシミュレーションの結果をプロットします。FakeToronto におけるノイズの深刻な影響が見て取れます。

plt.plot(angles, result_fake[0].data.evs, "o", label="Fake Device")
plt.plot(angles, result_sim[0].data.evs, "o", label="AerSimulator")
plt.xlabel("$\\mathrm{R_x}$ angle $\\theta_h$")
plt.title("$\\langle Z_{13} \\rangle$")
plt.legend()
plt.show()

Output of the previous code cell

5. 127量子ビット版の問題を解く​

冒頭で述べたユーティリティスケール実験を実行することが目標です。⟨Z62⟩\langle Z_{62} \rangle を観測量とした 127量子ビット × 60レイヤー（20トロッターステップ）の回路を作成・実行します。適切な箇所で27量子ビット版のコードを参考にしながら、自分で挑戦することをお勧めします。ただし、解答はここに示します。

解答：

5.1 ステップ1: Map​

# backend_map = service.backend("ibm_brisbane")
backend_map = service.least_busy(operational=True, simulator=False)

num_steps = 20
qc = get_utility_circuit(backend_map, num_steps)
obs = SparsePauliOp.from_sparse_list(
    [("Z", [62], 1)], num_qubits=backend_map.num_qubits
)  # Eagle
angles = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, np.pi / 2]

5.2 ステップ2と3: Optimize と Execute​

Eagle プロセッサの結合マップには144のエッジがあることに注意してください。

# backend = service.backend("ibm_brisbane")
backend = backend_map

transpiled_qc = transpile(qc, backend, optimization_level=1, layout_method="trivial")
transpiled_obs = obs.apply_layout(layout=transpiled_qc.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc))
156 qubits × 60 layers (221-depth),  Gate breakdown: SDG 7040, UNITARY 7040, CX 3520, R 3120
156 qubits × 60 layers (201-depth),  Gate breakdown: RZ 11933, SX 6240, CZ 3520
params = [[p] for p in angles]
estimator = Estimator(mode=backend)
pub = (transpiled_qc, transpiled_obs, params)
job = estimator.run([pub])

job_id = job.job_id()
print(f"job id={job_id}")
job id=d1479n6qf56g0081sxa0

5.3 Post-process​

ユーティリティ論文の Figure 4b における「緩和済み」の点の値を提供します。これらをあなたの結果と合わせてプロットしてください。

result_paper = [
    1.0171,
    1.0044,
    0.9563,
    0.9602,
    0.8394,
    0.8120,
    0.5466,
    0.4556,
    0.1953,
    0.0141,
    0.0117,
]

# REPLACE WITH YOUR OWN JOB ID
job = service.job(job_id)

plt.plot(angles, job.result()[0].data.evs, "o", label=f"{job.backend().name}")
plt.plot(angles, result_paper, "o", label="Utility Paper")
plt.xlabel("$\\mathrm{R_x}$ angle $\\theta_h$")
plt.title("$\\langle Z_{62} \\rangle$")
plt.legend()
plt.show()

Output of the previous code cell

あなたの結果は Figure 4b の「未緩和」の点に似ていますか？   実験時のデバイスの状態によって、結果は大きく異なる場合があります。結果自体は気にしないでください。確認するのは、コーディングを正しく行ったかどうかです。正しくできていれば、おめでとうございます。あなたはユーティリティ時代のスタートラインに立ちました。

ユーティリティ論文と同様に、世界中の科学者たちがノイズの存在下でも意味のある結果を得るために多大な工夫を凝らしてきました。この集合的な取り組みの最終目標は量子優位性です。それは、量子コンピュータが産業において有用ないくつかの問題を、古典コンピュータよりも速く、より高い精度で、またはより低コストで解けるようになる状態です。これは単一の出来事ではなく、古典的な手法による量子結果の再現に次第に長い時間がかかるようになり、やがてその量子的なリードタイムが決定的な重要性を持つようになる時代となるでしょう。量子優位性について一つ明らかなことがあります：それはユーティリティスケール実験を通じてのみ達成できるということです。もしこのコースが、挑戦とおもしろさに満ちたその探求にあなたを導く一助となれば、私たちにとってこれ以上の喜びはありません。

参考文献​