ショアのアルゴリズム
使用量の目安: Eagle r3プロセッサで約3秒(注意: これは目安です。実際の実行時間は異なる場合があります。)
ショアのアルゴリズムは、1994年にピーター・ショアによって開発された、整数を多項式時間で素因数分解する画期的な量子アルゴリズムです。その重要性は、既知のどの古典アルゴリズムよりも指数関数的に高速に大きな整数を素因数分解できる能力にあります。これは、大きな数の素因数分解の困難さに依存するRSAなどの広く使われている暗号システムのセキュリティを脅かすものです。十分に強力な量子コンピュータ上でこの問題を効率的に解くことにより、ショアのアルゴリズムは暗号技術、サイバーセキュリティ、計算数学などの分野に革命をもたらす可能性があり、量子計算の変革的な力を示しています。
このチュートリアルでは、量子コンピュータ上で15を素因数分解することにより、ショアのアルゴリズムを実演することに焦点を当てます。
まず、位数発見問題を定義し、量子位相推定プロトコルから対応する回路を構築します。次に、トランスパイル可能な最短深度の回路を用いて、実際のハードウェア上で位数発見回路を実行します。最後のセクションでは、位数発見問題と整数の素因数分解を結びつけることにより、ショアのアルゴリズムを完成させます。
チュートリアルの最後に、実際のハードウェア上でのショアのアルゴリズムの他の実演について議論します。汎用的な実装と、15や21などの特定の整数の素因数分解に特化した実装の両方 に焦点を当てます。 注意: このチュートリアルは、ショアのアルゴリズムに関する回路の実装と実演に重点を置いています。内容に関する詳細な学習リソースについては、ジョン・ワトラス博士による量子アルゴリズムの基礎コース、および参考文献セクションの論文をご参照ください。
前提条件
このチュートリアルを開始する前に、以下がインストールされていることを確認してください。
- Qiskit SDK v2.0以降、可視化サポート付き
- Qiskit Runtime v0.40以降 (
pip install qiskit-ibm-runtime)
セットアップ
# Added by doQumentation — required packages for this notebook
!pip install -q numpy pandas qiskit qiskit-ibm-runtime
import numpy as np
import pandas as pd
from fractions import Fraction
from math import floor, gcd, log
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit.transpiler import CouplingMap, generate_preset_pass_manager
from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler
ステップ1: 古典的な入力を量子問題にマッピングする
背景
ショアの整数素因数分解アルゴリズムは、位数発見問題として知られる中間的な問題を利用します。このセクションでは、量子位相推定を用いて位数発見問題を解く方法を示します。
位相推定問題
位相推定問題では、量子ビットの量子状態と、量子ビットに作用するユニタリ量子回路が与えられます。が回路の動作を記述するユニタリ行列の固有ベクトルであることが保証されており、目標はに対応する固有値を計算または近似することです。言い換えれば、回路は以下を満たす数の近似値を出力する必要があります。
位相推定回路の目標は、をビットで近似することです。数学的には、に対してとなるを求めたいということです。以下の図は、量子ビットの測定によりをビットで推定する量子回路を示しています。
上記の回路では、上位量子ビットは状態に初期化され、下位量子ビットはの固有ベクトルであることが保証されているに初期化されます。位相推定回路の最初の要素は、対応する制御量子ビットに対して位相キックバックを実行する制御ユニタリ演算です。これらの制御ユニタリは、最下位ビットから最上位ビットまで、制御量子ビットの位置に応じてべき乗されます。はの固有 ベクトルであるため、下位量子ビットの状態はこの演算の影響を受けませんが、固有値の位相情報が上位量子ビットに伝搬します。
制御ユニタリによる位相キックバック演算の後、上位量子ビットのすべての可能な状態は、ユニタリの各固有ベクトルに対して互いに直交正規であることがわかります。したがって、これらの状態は完全に区別可能であり、それらが形成する基底を計算基底に回転させて測定を行うことができます。数学的な解析により、この回転行列は次元ヒルベルト空間における逆量子フーリエ変換(QFT)に対応することが示されます。この直感的な理由は、モジュラべき乗演算子の周期的構造が量子状態にエンコードされており、QFTがこの周期性を周波数領域の測定可能なピークに変換するためです。
ショアのアルゴリズムでQFT回路が使用される理由のより深い理解については、量子アルゴリズムの基礎コースをご参照ください。 これで、位相推定回路を位数発見に使用する準備が整いました。
位数発見問題
位数発見問題を定義するために、いくつかの数論の概念から始めます。まず、任意の正の整数に対して、集合を次のように定義します。 におけるすべての算術演算はを法として行われます。特に、と互いに素なすべての要素は特別であり、を構成します。 要素に対して、を満たす最小の正の整数を、のを法とする位数と定義します。後で見るように、の位数を求めることでを素因数分解することができます。 位相推定回路から位数発見回路を構築するには、2つの考慮事項が必要です。第一に、位数を求めることを可能にするユニタリを定義する必要があり、第二に、位相推定回路の初期状態を準備するた めのの固有ベクトルを定義する必要があります。
位数発見問題を位相推定に結びつけるために、古典的な状態がに対応するシステム上で、固定された要素を乗じる演算を考えます。特に、この乗算演算子を、各に対してと定義します。ケットの右辺の内部でを法とする積を取っていることは暗黙的であることに注意してください。数学的な解析により、はユニタリ演算子であることが示されます。さらに、はの位数を位相推定問題に結びつけることを可能にする固有ベクトルと固有値のペアを持つことがわかります。具体的には、の任意の選択に対して、はの固有ベクトルであり、対応する固有値はです。ここで、です。 観察により、便利な固有ベクトル/固有値のペアは、を持つ状態であ ることがわかります。したがって、固有ベクトルを見つけることができれば、量子回路で位相を推定し、位数の推定値を得ることができます。しかし、これは容易ではなく、代替案を検討する必要があります。
初期状態として計算基底状態を準備した場合に回路がどうなるかを考えてみましょう。これはの固有状態ではありませんが、上で説明した固有状態の一様な重ね合わせです。言い換えれば、以下の関係が成り立ちます。 上記の方程式の意味は、初期状態をに設定した場合、から一様にランダムにを選択し、位相推定回路の固有ベクトルとしてを使用した場合とまったく同じ測定結果が得られるということです。つまり、上位量子ビットの測定は 、が一様にランダムに選ばれた値の近似値を与えます。これにより、複数回の独立した実行の後に高い確度でを学習することができ、これが我々の目標でした。
モジュラべき乗演算子
ここまでで、量子回路においておよびを定義することにより、位相推定問題を位数発見問題に結びつけました。したがって、残る最後の要素は、に対してのモジュラべき乗を効率的に定義する方法を見つけることです。 この計算を行うために、任意のべき乗に対して、の回路を回繰り返すのではなく、を計算してからの回路を使用することでの回路を作成できることがわかります。必要なべき乗は2のべき乗のみであるため、反復二乗法を用いて古典的に効率よく計算できます。
ステップ2: 量子ハードウェア実行のための問題の最適化
、の具体的な例
ここで一旦立ち止まって、具体的な例について議論し、の位数発見回路を構築しましょう。に対する可能な 非自明なはであることに注意してください。この例ではを選びます。演算子とモジュラべき乗演算子を構築します。 の計算基底状態への作用は以下の通りです。 観察により、基底状態がシャッフルされていることがわかるため、置換行列が得られます。この演算はスワップゲートを用いて4量子ビット上で構築できます。以下では、および制御演算を構築します。
def M2mod15():
"""
M2 (mod 15)
"""
b = 2
U = QuantumCircuit(4)
U.swap(2, 3)
U.swap(1, 2)
U.swap(0, 1)
U = U.to_gate()
U.name = f"M_{b}"
return U
# Get the M2 operator
M2 = M2mod15()
# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)
def controlled_M2mod15():
"""
Controlled M2 (mod 15)
"""
b = 2
U = QuantumCircuit(4)
U.swap(2, 3)
U.swap(1, 2)
U.swap(0, 1)
U = U.to_gate()
U.name = f"M_{b}"
c_U = U.control()
return c_U
# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()
# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)
3量子ビット以上に作用するゲートは、さらに2量子ビットゲートに分解されます。
circ.decompose(reps=2).draw(output="mpl", fold=-1)

次に、モジュラべき乗演算子を構築する必要があります。位相推定で十分な精度を得るために、推定測定に8量子ビットを使用します。したがって、各に対してとなるを構築する必要があります。
def a2kmodN(a, k, N):
"""Compute a^{2^k} (mod N) by repeated squaring"""
for _ in range(k):
a = int(np.mod(a**2, N))
return a
k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]
print(b_list)
[2, 4, 1, 1, 1, 1, 1, 1]
の値のリストからわかるように、先ほど構築したに加えて、とも構築する必要があります。は計算基底状態に対して自明に作用するため、単純に恒等演算子であることに注意してください。
の計算基底状態への作用は以下の通りです。
したがって、この置換は以下のスワップ演算で構築できます。
def M4mod15():
"""
M4 (mod 15)
"""
b = 4
U = QuantumCircuit(4)
U.swap(1, 3)
U.swap(0, 2)
U = U.to_gate()
U.name = f"M_{b}"
return U
# Get the M4 operator
M4 = M4mod15()
# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)
def controlled_M4mod15():
"""
Controlled M4 (mod 15)
"""
b = 4
U = QuantumCircuit(4)
U.swap(1, 3)
U.swap(0, 2)
U = U.to_gate()
U.name = f"M_{b}"
c_U = U.control()
return c_U
# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()
# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)
3量子ビット以上に作用するゲートは、さらに2量子ビットゲートに分解されます。
circ.decompose(reps=2).draw(output="mpl", fold=-1)

演算子は、与えられたに対して置換演算であることを確認しました。ここではで4量子ビットしか必要としないため、置換問題が比較的小規模であることから、これらの演算を目視検査によりSWAPゲートで直接合成することができました。一般的には、これはスケーラブルなアプローチではない可能性があります。代わりに、置換行列を明示的に構築し、QiskitのUnitaryGateクラスとトランスパイル手法を使用してこの置換行列を合成する必要があるかもしれません。ただし、これにより回路の深度が大幅に増加する可能性があります。以下に例を示します。
def mod_mult_gate(b, N):
"""
Modular multiplication gate from permutation matrix.
"""
if gcd(b, N) > 1:
print(f"Error: gcd({b},{N}) > 1")
else:
n = floor(log(N - 1, 2)) + 1
U = np.full((2**n, 2**n), 0)
for x in range(N):
U[b * x % N][x] = 1
for x in range(N, 2**n):
U[x][x] = 1
G = UnitaryGate(U)
G.name = f"M_{b}"
return G
# Let's build M2 using the permutation matrix definition
M2_other = mod_mult_gate(2, 15)
# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2_other, inplace=True)
circ = circ.decompose()
# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)
print(f"qubits: {circ.num_qubits}")
print(
f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.decompose().draw(
output="mpl", fold=-1, style="clifford", idle_wires=False
)
qubits: 4
2q-depth: 94
2q-size: 96
Operator counts: OrderedDict({'cx': 45, 'swap': 32, 'u': 24, 'u1': 7, 'u3': 4, 'unitary': 3, 'circuit-335': 1, 'circuit-338': 1, 'circuit-341': 1, 'circuit-344': 1, 'circuit-347': 1, 'circuit-350': 1, 'circuit-353': 1, 'circuit-356': 1, 'circuit-359': 1, 'circuit-362': 1, 'circuit-365': 1, 'circuit-368': 1, 'circuit-371': 1, 'circuit-374': 1, 'circuit-377': 1, 'circuit-380': 1})

手動で実装したゲートのコンパイル済み回路の深度と、これらのカウントを比較してみましょう。
# Get the M2 operator from our manual construction
M2 = M2mod15()
# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ = circ.decompose(reps=3)
# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)
print(f"qubits: {circ.num_qubits}")
print(
f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.draw(
output="mpl", fold=-1, style="clifford", idle_wires=False
)
qubits: 4
2q-depth: 9
2q-size: 9
Operator counts: OrderedDict({'cx': 9})
ご覧の通り、置換行列によるアプローチは、単一のゲートであっても手動実装と比較して大幅に深い回路になりました。したがって、以前の演算の実装を引き続き使用します。 次に、先ほど定義した制御モジュラべき乗演算子を使用して、完全な位数発見回路を構築する準備が整いました。以下のコードでは、Qiskit回路ライブラリからQFT回路もインポートしています。この回路は各量子ビットにアダマールゲート、一連の制御U1(または位相に応じてZ)ゲート、およびスワップゲートの層を使用します。
# Order finding problem for N = 15 with a = 2
N = 15
a = 2
# Number of qubits
num_target = floor(log(N - 1, 2)) + 1 # for modular exponentiation operators
num_control = 2 * num_target # for enough precision of estimation
# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]
# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)
# Initialize the target register to the state |1>
circuit.x(num_control)
# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
circuit.h(k)
b = b_list[k]
if b == 2:
circuit.compose(
M2mod15().control(), qubits=[qubit] + list(target), inplace=True
)
elif b == 4:
circuit.compose(
M4mod15().control(), qubits=[qubit] + list(target), inplace=True
)
else:
continue # M1 is the identity operator
# Apply the inverse QFT to the control register
circuit.compose(QFT(num_control, inverse=True), qubits=control, inplace=True)
# Measure the control register
circuit.measure(control, output)
circuit.draw("mpl", fold=-1)
が恒等演算子であるため、残りの制御量子ビットからの制御モジュラべき乗演算を省略したことに注意してください。
このチュートリアルの後半では、この回路をibm_marrakeshバックエンドで実行します。これを行うために、この特定のバックエンドに合わせて回路をトランスパイルし、回路の深度とゲート数を報告します。
service = QiskitRuntimeService()
backend = service.backend("ibm_marrakesh")
pm = generate_preset_pass_manager(optimization_level=2, backend=backend)
transpiled_circuit = pm.run(circuit)
print(
f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}"
)
print(
f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}"
)
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(
output="mpl", fold=-1, style="clifford", idle_wires=False
)
2q-depth: 187
2q-size: 260
Operator counts: OrderedDict({'sx': 521, 'rz': 354, 'cz': 260, 'measure': 8, 'x': 4})
