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

Q-CTRLのQiskit Functionsを用いた量子位相推定

使用量の目安: Heron r2プロセッサで約40秒。(注: これは推定値です。実際の実行時間は異なる場合があります。)

背景

量子位相推定(QPE)は、量子コンピューティングにおける基礎的なアルゴリズムであり、Shorのアルゴリズム、量子化学における基底状態エネルギー推定、固有値問題など、多くの重要な応用の基盤となっています。QPEは、ユニタリ演算子の固有状態に関連する位相φ\varphiを推定するもので、以下の関係式に符号化されています。

Uφ=e2πiφφ,U \lvert \varphi \rangle = e^{2\pi i \varphi} \lvert \varphi \rangle,

mm個のカウント量子ビットを用いてϵ=O(1/2m)\epsilon = O(1/2^m)の精度で位相を決定します[1]。これらの量子ビットを重ね合わせ状態に準備し、UUの制御べき乗を適用した後、逆量子フーリエ変換(QFT)を使用して位相を二進符号化された測定結果として抽出します。QPEは、二進分数がφ\varphiを近似するビット列にピークを持つ確率分布を生成します。理想的な場合、最も確率の高い測定結果は位相の二進展開に直接対応し、他の結果の確率はカウント量子ビット数の増加とともに急速に減少します。しかし、深いQPE回路をハードウェア上で実行するには課題があります。多数の量子ビットとエンタングリング操作により、アルゴリズムはデコヒーレンスやゲートエラーに対して非常に敏感になります。その結果、ビット列の分布が広がったりシフトしたりして、真の固有位相が隠されてしまいます。結果として、最も確率の高いビット列がφ\varphiの正しい二進展開に対応しなくなる可能性があります。

このチュートリアルでは、Q-CTRLのFire Opalエラー抑制およびパフォーマンス管理ツールを使用したQPEアルゴリズムの実装を紹介します。これはQiskit Functionとして提供されています(Fire Opalドキュメントを参照)。Fire Opalは、動的デカップリング、量子ビットレイアウトの改善、エラー抑制技術などの高度な最適化を自動的に適用し、より高忠実度の結果を実現します。これらの改善により、ハードウェアのビット列分布がノイズのないシミュレーションで得られる分布に近づき、ノイズの影響下でも正しい固有位相を確実に特定できるようになります。

前提条件

このチュートリアルを始める前に、以下がインストールされていることを確認してください。

  • Qiskit SDK v1.4以降(visualizationサポート付き)
  • Qiskit Runtime v0.40以降(pip install qiskit-ibm-runtime
  • Qiskit Functions Catalog v0.9.0(pip install qiskit-ibm-catalog
  • Fire Opal SDK v9.0.2以降(pip install fire-opal
  • Q-CTRL Visualizer v8.0.2以降(pip install qctrl-visualizer

セットアップ

まず、IBM Quantum APIキーを使用して認証します。次に、以下のようにQiskit Functionを選択します。(このコードは、ローカル環境にアカウントが保存済みであることを前提としています。)

# Added by doQumentation — installs packages not in the Binder environment
%pip install -q qctrlvisualizer
from qiskit import QuantumCircuit

import numpy as np
import matplotlib.pyplot as plt
import qiskit
from qiskit import qasm2
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
import qctrlvisualizer as qv
from qiskit_ibm_catalog import QiskitFunctionsCatalog

plt.style.use(qv.get_qctrl_style())
catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Access Function
perf_mgmt = catalog.load("q-ctrl/performance-management")

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

このチュートリアルでは、既知の単一量子ビットユニタリの固有位相を復元するためにQPEを説明します。推定したい位相を持つユニタリは、ターゲット量子ビットに適用される単一量子ビット位相ゲートです。

U(θ)=(100eiθ)=eiθ1 ⁣1.U(\theta)= \begin{pmatrix} 1 & 0\\[2pt] 0 & e^{i\theta} \end{pmatrix} = e^{i\theta\,|1\rangle\!\langle 1|}.

固有状態ψ=1|\psi\rangle=|1\rangleを準備します。1|1\rangleは固有値eiθe^{i\theta}を持つU(θ)U(\theta)の固有ベクトルであるため、推定すべき固有位相は以下のようになります。

φ=θ2π(mod1)\varphi = \frac{\theta}{2\pi} \pmod{1}

θ=162π\theta=\tfrac{1}{6}\cdot 2\piと設定するため、真の位相はφ=1/6\varphi=1/6となります。QPE回路は、角度θ2k\theta\cdot2^kの制御位相回転を適用することで制御べき乗U2kU^{2^k}を実装し、その後カウントレジスタに逆QFTを適用して測定します。得られるビット列は1/61/6の二進表現の周辺に集中します。

回路はmm個のカウント量子ビット(推定精度を設定するため)と1個のターゲット量子ビットを使用します。まず、QPEの実装に必要な構成要素、すなわち量子フーリエ変換(QFT)とその逆変換、固有位相の10進数と二進分数の間の変換を行うユーティリティ関数、シミュレーション結果とハードウェア結果を比較するために生のカウントを確率に正規化するヘルパー関数を定義します。

def inverse_quantum_fourier_transform(quantum_circuit, number_of_qubits):
"""
Apply an inverse Quantum Fourier Transform the first `number_of_qubits` qubits in the
`quantum_circuit`.
"""
for qubit in range(number_of_qubits // 2):
quantum_circuit.swap(qubit, number_of_qubits - qubit - 1)
for j in range(number_of_qubits):
for m in range(j):
quantum_circuit.cp(-np.pi / float(2 ** (j - m)), m, j)
quantum_circuit.h(j)
return quantum_circuit
def bitstring_count_to_probabilities(data, shot_count):
"""
This function turns an unsorted dictionary of bitstring counts into a sorted dictionary
of probabilities.
"""
# Turn the bitstring counts into probabilities.
probabilities = {
bitstring: bitstring_count / shot_count
for bitstring, bitstring_count in data.items()
}

sorted_probabilities = dict(
sorted(probabilities.items(), key=lambda x: x[1], reverse=True)
)

return sorted_probabilities

ステップ2: 量子ハードウェア実行のための問題の最適化

カウント量子ビットを重ね合わせ状態に準備し、制御位相回転を適用してターゲットの固有位相を符号化し、測定前に逆QFTで仕上げることで、QPE回路を構築します。

def quantum_phase_estimation_benchmark_circuit(
number_of_counting_qubits, phase
):
"""
Create the circuit for quantum phase estimation.

Parameters
----------
number_of_counting_qubits : The number of qubits in the circuit.
phase : The desired phase.

Returns
-------
QuantumCircuit
The quantum phase estimation circuit for `number_of_counting_qubits` qubits.
"""
qc = QuantumCircuit(
number_of_counting_qubits + 1, number_of_counting_qubits
)
target = number_of_counting_qubits

# |1> eigenstate for the single-qubit phase gate
qc.x(target)

# Hadamards on counting register
for q in range(number_of_counting_qubits):
qc.h(q)

# ONE controlled phase per counting qubit: cp(phase * 2**k)
for k in range(number_of_counting_qubits):
qc.cp(phase * (1 << k), k, target)

qc.barrier()

# Inverse QFT on counting register
inverse_quantum_fourier_transform(qc, number_of_counting_qubits)

qc.barrier()
for q in range(number_of_counting_qubits):
qc.measure(q, q)
return qc

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

実験のショット数と量子ビット数を設定し、mm桁の二進数を使用してターゲット位相φ=1/6\varphi = 1/6を符号化します。これらのパラメータを用いて、シミュレーション、デフォルトのハードウェア、およびFire Opal強化バックエンドで実行するQPE回路を構築します。

shot_count = 10000
num_qubits = 35
phase = (1 / 6) * 2 * np.pi
circuits_quantum_phase_estimation = (
quantum_phase_estimation_benchmark_circuit(
number_of_counting_qubits=num_qubits, phase=phase
)
)

MPSシミュレーションの実行

まず、matrix_product_stateシミュレータを使用してリファレンス分布を生成し、カウントを正規化された確率に変換して、後でハードウェア結果と比較できるようにします。

# Run the algorithm on the IBM Aer simulator.
aer_simulator = AerSimulator(method="matrix_product_state")

# Transpile the circuits for the simulator.
transpiled_circuits = qiskit.transpile(
circuits_quantum_phase_estimation, aer_simulator
)
simulated_result = (
aer_simulator.run(transpiled_circuits, shots=shot_count)
.result()
.get_counts()
)
simulated_result_probabilities = []

simulated_result_probabilities.append(
bitstring_count_to_probabilities(
simulated_result,
shot_count=shot_count,
)
)

ハードウェアでの実行

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

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_circuits = pm.run(circuits_quantum_phase_estimation)
# Run the algorithm with IBM default.
sampler = Sampler(backend)

# Run all circuits using Qiskit Runtime.
ibm_default_job = sampler.run([isa_circuits], shots=shot_count)

Fire Opalを使用したハードウェアでの実行

# Run the circuit using the sampler
fire_opal_job = perf_mgmt.run(
primitive="sampler",
pubs=[qasm2.dumps(circuits_quantum_phase_estimation)],
backend_name=backend.name,
options={"default_shots": shot_count},
)

ステップ4: 後処理と所望の古典形式での結果の返却

# Retrieve results.
ibm_default_result = ibm_default_job.result()
ibm_default_probabilities = []

for idx, pub_result in enumerate(ibm_default_result):
ibm_default_probabilities.append(
bitstring_count_to_probabilities(
pub_result.data.c0.get_counts(),
shot_count=shot_count,
)
)
fire_opal_result = fire_opal_job.result()

fire_opal_probabilities = []
for idx, pub_result in enumerate(fire_opal_result):
fire_opal_probabilities.append(
bitstring_count_to_probabilities(
pub_result.data.c0.get_counts(),
shot_count=shot_count,
)
)
data = {
"simulation": simulated_result_probabilities,
"default": ibm_default_probabilities,
"fire_opal": fire_opal_probabilities,
}
def plot_distributions(
data,
number_of_counting_qubits,
top_k=None,
by="prob",
shot_count=None,
):
def nrm(d):
s = sum(d.values())
return {k: (v / s if s else 0.0) for k, v in d.items()}

def as_float(d):
return {k: float(v) for k, v in d.items()}

def to_space(d):
if by == "prob":
return nrm(as_float(d))
else:
if shot_count and 0.99 <= sum(d.values()) <= 1.01:
return {
k: v * float(shot_count) for k, v in as_float(d).items()
}
else:
return as_float(d)

def topk(d, k):
items = sorted(d.items(), key=lambda kv: kv[1], reverse=True)
return items[: (k or len(d))]

phase = "1/6"

sim = to_space(data["simulation"])
dft = to_space(data["default"])
qct = to_space(data["fire_opal"])

correct = max(sim, key=sim.get) if sim else None
print("Correct result:", correct)

sim_items = topk(sim, top_k)
dft_items = topk(dft, top_k)
qct_items = topk(qct, top_k)

sim_keys, y_sim = zip(*sim_items) if sim_items else ([], [])
dft_keys, y_dft = zip(*dft_items) if dft_items else ([], [])
qct_keys, y_qct = zip(*qct_items) if qct_items else ([], [])

fig, axes = plt.subplots(3, 1, layout="constrained")
ylab = "Probabilities"

def panel(ax, keys, ys, title, color):
x = np.arange(len(keys))
bars = ax.bar(x, ys, color=color)
ax.set_title(title)
ax.set_ylabel(ylab)
ax.set_xticks(x)
ax.set_xticklabels(keys, rotation=90)
ax.set_xlabel("Bitstrings")
if correct in keys:
i = keys.index(correct)
bars[i].set_edgecolor("black")
bars[i].set_linewidth(2)
return max(ys, default=0.0)

c_sim, c_dft, c_qct = (
qv.QCTRL_STYLE_COLORS[5],
qv.QCTRL_STYLE_COLORS[1],
qv.QCTRL_STYLE_COLORS[0],
)
m1 = panel(axes[0], list(sim_keys), list(y_sim), "Simulation", c_sim)
m2 = panel(axes[1], list(dft_keys), list(y_dft), "Default", c_dft)
m3 = panel(axes[2], list(qct_keys), list(y_qct), "Q-CTRL", c_qct)

for ax, m in zip(axes, (m1, m2, m3)):
ax.set_ylim(0, 1.05 * (m or 1.0))

for ax in axes:
ax.label_outer()
fig.suptitle(
rf"{number_of_counting_qubits} counting qubits, $2\pi\varphi$={phase}"
)
fig.set_size_inches(20, 10)
plt.show()
experiment_index = 0
phase_index = 0

distributions = {
"simulation": data["simulation"][phase_index],
"default": data["default"][phase_index],
"fire_opal": data["fire_opal"][phase_index],
}

plot_distributions(
distributions, num_qubits, top_k=100, by="prob", shot_count=shot_count
)
Correct result: 00101010101010101010101010101010101

前のコードセルの出力

シミュレーションは正しい固有位相のベースラインを設定します。デフォルトのハードウェア実行では、ノイズが確率を多くの不正確なビット列に分散させるため、結果が不明瞭になります。Q-CTRLパフォーマンス管理を使用すると、分布がより鮮明になり、正しい結果が復元されるため、このスケールでの信頼性の高いQPEが可能になります。

参考文献

[1] Lecture 7: Phase Estimation and Factoring. IBM Quantum Learning - Fundamentals of quantum algorithms. Retrieved October 3, 2025.

チュートリアルアンケート

このチュートリアルに関するフィードバックをお寄せください。皆様のご意見は、コンテンツの改善とユーザーエクスペリエンスの向上に役立てさせていただきます。

アンケートへのリンク