Skip to main content
IBM Quantum Platform


title: "クイック・スタート" description: "pauli-prop Qiskit アドオンパッケージのクイックスタートガイド"

クイック・スタート

このガイドでは、 1D スピン鎖上の10量子ビットのキック付きイジングモデルの時間的ダイナミクスを、パッケージ pauli-prop を用いて古典的にシミュレーションします。


パウリ伝播のための入力の準備

検討対象となるハミルトニアンは以下の通りである:

H=−J∑⟨i,j⟩ZiZj+h∑iXiH = -J\sum\limits_{\langle i,j \rangle} Z_iZ_j + h\sum\limits_iX_i

ここで、 J>0J>0 は最近接スピン間の結合を表し、 i<ji<j、 hh は全横磁場である。 時間発展演算子の1次トロッター分解は、 2020 回のトロッターステップにわたる量子回路 UU として実装される。 結合定数 JJ は J=−π2J=-\frac{\pi}{2} に、 hh は π6\frac{\pi}{6} にそれぞれ固定される。 ZZZZ の相互作用は、クリフォードゲートを用いて実装される( CXCX, SdgSdg, Y\sqrt{Y} )。

トロッター化された時間発展を量子回路として実装し、x軸周りの非クリフォード回転には π6\frac{\pi}{6} を用いる。 これらの角度がクリフォード角から遠ざかるほど(例えば、 θ=nπ2,n∈Z\theta=n\frac{\pi}{2}, n \in \mathbb{Z} のように)、パウリ伝播法を用いたシミュレーションはより困難になります。

観測量の選択については、単一サイトにおける平均磁化 1N∑i=1N⟨zi⟩\frac{1}{N} \sum_{i=1}^{N} \langle z_i \rangle を考慮する。ここで、 NN はスピンの数である。

import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler import CouplingMap

num_qubits = 10
coupling_map = CouplingMap.from_line(num_qubits, bidirectional=False)

# Num Trotter steps
num_steps = 20
theta_rx = np.pi / 6

# Average single-site magnetization
observable = (
    SparsePauliOp(
        [
            "I" * iq + "Z" + "I" * (num_qubits - iq - 1)
            for iq in range(num_qubits)
        ]
    )
    / num_qubits
)

# Create the Trotter circuit
num_qubits = 10
num_steps = 20
theta_rx = np.pi / 6
circuit = QuantumCircuit(num_qubits)
edges = CouplingMap.from_line(num_qubits, bidirectional=False).get_edges()
for _ in range(num_steps):
    circuit.rx(theta_rx, [i for i in range(num_qubits)])
    for edge in edges:
        circuit.sdg(edge)
        circuit.ry(np.pi / 2, edge[1])
        circuit.cx(edge[0], edge[1])
        circuit.ry(-np.pi / 2, edge[1])
circuit.draw("mpl", fold=-1)

Output:

Output of the previous code cell

パウリ伝播を用いて、系の時間発展をシミュレートする

回路( UU )と観測可能量( OO )が準備できたら、以下の数ステップでシステムを簡単にシミュレーションできます:

  • UU を、Clifford 部分 CC と非 Clifford 部分 PP に分割し、 U=PCU=PC が成り立つようにする。この際、以下の式を用いる。 evolve_through_cliffords
  • OO を PP を通じて展開すると、新しい演算子 O′O^\prime が得られます。その方法は以下の通りです。 pauli_prop.propagate_through_circuit
  • Qiskitに組み込まれているクリフォード進化のサポート機能を使用して、回路のクリフォード部分において O′O^\prime を進化させる
  • O′O^\prime における、完全対角パウリ項(すべての量子ビット Z で I または のいずれかを含むパウリ項)に関連する係数を合計することで、期待値を ⟨0∣O′∣0⟩≈⟨0∣U†OU∣0⟩\langle0|O^\prime|0\rangle \approx \langle0|U^\dagger OU|0\rangle と近似する。 なお、これは近似であることに注意してください。これは、 O′O^\prime の項を、回路の非クリフォード部分を通じて伝播させる際に切り捨てたためです。
import time

from pauli_prop import evolve_through_cliffords, propagate_through_circuit

cliff, non_cliff = evolve_through_cliffords(circuit)

max_terms_list = [10**i for i in range(8)]
approx_evs = []
durations = []
for max_terms in max_terms_list:
    st = time.perf_counter()
    evolved_obs = propagate_through_circuit(
        observable, non_cliff, max_terms=max_terms, atol=1e-12, frame="h"
    )[0]
    evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
    durations.append(time.perf_counter() - st)
    approx_evs.append(
        float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())
    )

より大規模な計算を行うにつれて、期待値による近似の精度は高まります。 この例では、 410≈1064^{10}\approx10^6 付近でパウリ空間全体が飽和しており、これが最後の2つの点の間で曲線が平坦になることに表れています。

以下のプロットは単調収束を示していますが、パウリ伝播のシミュレーションは一般的に単調に収束するわけではありません。 この種のプロットでは、「ぎくしゃくした」動きが見られることは珍しくありません。

import matplotlib.pyplot as plt
from qiskit_aer import AerSimulator

sim_circ = circuit.copy()
sim_circ.save_statevector()
backend = AerSimulator(method="statevector")
psi = backend.run(sim_circ).result().data()["statevector"]
exact_ev = psi.expectation_value(observable)

ax1 = plt.gca()
ax1.plot(max_terms_list, approx_evs, marker="o", label="Approximate")
ax1.axhline(exact_ev, linestyle="--", color="green", label="Exact")
ax1.set_xscale("log")
ax1.set_xlabel("# terms kept")
ax1.set_ylabel(r"$\frac{1}{N} \sum_{i=1}^{N} \langle z_i \rangle$")

ax2 = ax1.twinx()
ax2.plot(
    max_terms_list, durations, marker=".", label="Runtime", color="orange"
)
ax2.set_ylabel("Runtime (s)", color="orange")
ax2.set_yscale("log")

handles1, labels1 = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles1 + handles2, labels1 + labels2, loc="lower right")

plt.title(f"Simulating {num_steps}-step 1D Ising Model")

Output:

Text(0.5, 1.0, 'Simulating 20-step 1D Ising Model')
Output of the previous code cell
このページは役に立ちましたか?
バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。