量子アルゴリズム：位相推定

備考

Kento Ueda（2024年5月15日）

このノートブックでは、量子フーリエ変換（QFT）と量子位相推定（QPE）の基本概念と実装について説明します。

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

この実験の実行に必要な QPU 時間の目安は約 7 秒です。

1. はじめに

量子フーリエ変換（QFT）

量子フーリエ変換は、古典的な離散フーリエ変換の量子版です。量子状態に適用される線形変換であり、計算基底をフーリエ基底表現にマッピングします。QFT は多くの量子アルゴリズムにおいて重要な役割を担っており、量子状態から周期性の情報を効率的に抽出する手法を提供します。 QFT は、$N$ 量子ビットに対してアダマールゲートや制御位相ゲートなどの量子ゲートを用いて、$O(N^2)$ 回の演算で実装でき、古典的なフーリエ変換に対する指数的な高速化を実現します。

• 応用例：大きな整数の素因数分解や離散対数を求める Shor のアルゴリズムなど、量子アルゴリズムの基盤となっています。

量子位相推定（QPE）

量子位相推定は、ユニタリ演算子の固有ベクトルに対応する位相を推定するための量子アルゴリズムです。このアルゴリズムは、量子状態の抽象的な数学的性質と計算上の応用との橋渡しをします。

• 応用例：ユニタリ行列の固有値の探索や量子システムのシミュレーションなどの問題を解くことができます。

QFT と QPE は、古典コンピュータでは現実的に解けない問題を解く多くの量子アルゴリズムにおいて、本質的な基盤となっています。このノートブックを通じて、これらの技法がどのように実装されるかを理解できるようになります。

2. 量子フーリエ変換（QFT）の基礎

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram, plot_bloch_multivector
from qiskit.quantum_info import Statevector
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import Sampler

from numpy import pi

離散フーリエ変換との類比から、QFT は $N$ 量子ビットの量子状態 $\vert X\rangle = \sum_{j=0}^{N-1} x_j \vert j \rangle$ に作用し、量子状態 $\vert Y\rangle = \sum_{k=0}^{N-1} y_k \vert k \rangle$ にマッピングします。

$$QFT_{N}: \vert j \rangle \mapsto \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}\omega_N^{jk} \vert k \rangle$$

ここで $\omega_N^{jk} = e^{2\pi i \frac{jk}{N}}$ です。

またはユニタリ行列表現として：

$$U_{QFT} = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} \sum_{k=0}^{N-1} \omega_N^{jk} \vert k \rangle \langle j \vert$$

2.1. 直感

量子フーリエ変換（QFT）は、計算基底（Z 基底）とフーリエ基底という 2 つの基底間の変換を行います。ではここでのフーリエ基底とはどのような意味でしょうか？フーリエ変換 $F(\omega)$ は関数 $f(x)$ を周波数 $\omega$ の正弦波関数に畳み込んだものを表すことを思い出してください。平たく言えば、フーリエ変換は $f(x)$ を正弦波（または余弦波）で構成するために各周波数 $\omega$ がどれだけ必要かを示す関数です。QFT がこの文脈でどのような意味を持つかをより直感的に理解するため、4 つの量子ビットを使って整数を 2 進数で符号化した以下のステップ画像を参照してください：

計算基底では、状態 $|0\rangle$ と $|1\rangle$ を使って数を 2 進数で格納します。

異なる量子ビットが変化する頻度に注目してください。一番左の量子ビットは数が 1 増えるごとに反転し、次の量子ビットは 2 増えるごと、3 番目は 4 増えるごとに反転し、以降も同様です。

これらの状態に量子フーリエ変換を適用すると、次のようにマッピングされます

$$|\text{State in Computational Basis}\rangle \quad \xrightarrow[]{\text{QFT}} \quad |\text{State in Fourier Basis}\rangle$$ $$\text{QFT}|x\rangle = |\widetilde{x}\rangle$$

（フーリエ基底における状態を表すのに、チルダ（~）を使うことがよくあります）。

フーリエ基底では、Z 軸周りの異なる回転を使って数を格納します：

フーリエ基底では、Z 軸周りの異なる回転を使って数を格納します：

保存したい数値によって、各量子ビットが Z 軸周りにどの角度で回転するかが決まります。状態 $|\widetilde{0}\rangle$ では、すべての量子ビットが状態 $|{+}\rangle$ にあります。上記の例で見たように、状態 $|\widetilde{5}\rangle$ を4量子ビットにエンコードするために、最も左の量子ビットを $\tfrac{5}{2^n} = \tfrac{5}{16}$ 回転（$\tfrac{5}{16}\times 2\pi$ ラジアン）させました。次の量子ビットはこれの2倍（$\tfrac{10}{16}\times 2\pi$ ラジアン、つまり $10/16$ 回転）回転させ、その次の量子ビットではさらにこの角度が2倍になり、以降も同様に続きます。

改めて、各量子ビットが変化する頻度に注目してください。この場合、最も左の量子ビット（qubit 0）が最も低い周波数を持ち、最も右の量子ビットが最も高い周波数を持ちます。

2.2 例：1 量子ビット QFT

$N = 2^1 = 2$ の場合を考えます。

ユニタリ行列は次のように書けます：

$$\begin{aligned} U_{QFT} & = \frac{1}{\sqrt{2}} \sum_{j=0}^{1} \sum_{k=0}^{1} \omega_2^{jk} \vert k \rangle \langle j \vert \\ & = \frac{1}{\sqrt{2}} (\omega_2^{0} \vert 0 \rangle \langle 0 \vert + \omega_2^{0} \vert 0 \rangle \langle 1 \vert + \omega_2^{0} \vert 1 \rangle \langle 0 \vert + \omega_2^{1} \vert 1 \rangle \langle 1 \vert) \\ & = \frac{1}{\sqrt{2}} (\vert 0 \rangle \langle 0 \vert + \vert 0 \rangle \langle 1 \vert + \vert 1 \rangle \langle 0 \vert - \vert 1 \rangle \langle 1 \vert) \\ & = H \end{aligned}$$

この演算はアダマールゲート（$H$）を適用した結果です。

2.3 QFT の積表現

$N = 2^n$ に対する変換 $QFT_{N}$ を、状態 $\vert x \rangle = \vert x_1\ldots x_n \rangle$ に作用させる形で一般化してみましょう。

$$\begin{aligned} QFT_N\vert x \rangle & = \frac{1}{\sqrt{N}} \sum_{y=0}^{N-1}\omega_N^{xy} \vert y \rangle \\ & = \frac{1}{\sqrt{N}} \sum_{y=0}^{N-1} e^{2 \pi i xy / 2^n} \vert y \rangle ~\text{since}\: \omega_N^{xy} = e^{2\pi i \frac{xy}{N}} \:\text{and}\: N = 2^n \\ & = \frac{1}{\sqrt{N}} \sum_{y=0}^{N-1} e^{2 \pi i \left(\sum_{k=1}^n y_k/2^k\right) x} \vert y_1 \ldots y_n \rangle \:\text{rewriting in fractional binary notation}\: y = y_1\ldots y_n, y/2^n = \sum_{k=1}^n y_k/2^k \\ & = \frac{1}{\sqrt{N}} \sum_{y=0}^{N-1} \prod_{k=1}^n e^{2 \pi i x y_k/2^k } \vert y_1 \ldots y_n \rangle \:\text{after expanding the exponential of a sum to a product of exponentials,} \\ & = \frac{1}{\sqrt{N}} \bigotimes_{k=1}^n \left(\vert0\rangle + e^{2 \pi i x /2^k } \vert1\rangle \right) \:\text{after rearranging the sum and products, and expanding} \sum_{y=0}^{N-1} = \sum_{y_1=0}^{1}\sum_{y_2=0}^{1}\ldots\sum_{y_n=0}^{1} \\ & = \frac{1}{\sqrt{N}} \left(\vert0\rangle + e^{\frac{2\pi i}{2}x} \vert1\rangle\right) \otimes \left(\vert0\rangle + e^{\frac{2\pi i}{2^2}x} \vert1\rangle\right) \otimes \ldots \otimes \left(\vert0\rangle + e^{\frac{2\pi i}{2^{n-1}}x} \vert1\rangle\right) \otimes \left(\vert0\rangle + e^{\frac{2\pi i}{2^n}x} \vert1\rangle\right) \end{aligned}$$

2.4 例：3 量子ビット QFT の回路構成

上記の説明だけでは、QFT 回路をどのように構成するかは明確でないかもしれません。ここでは、3 つの量子ビットの位相がそれぞれ異なる「速度」で変化すると予想される、という点だけを押さえておいてください。これが回路にどのように変換されるかを正確に理解することは、読者への演習問題とします。講義 PDF には複数の図と例が記載されています。また、「量子アルゴリズムの基礎」コースのこのレッスンも参考になります。

ここではシミュレータのみを使って QFT を示すため、量子位相推定に移るまでは Qiskit パターンのフレームワークは使用しません。

# Prepare for 3 qubits circuit
qr = QuantumRegister(3)
cr = ClassicalRegister(3)
qc = QuantumCircuit(qr, cr)

qc.h(2)
qc.cp(pi / 2, 1, 2)  # Controlled-phase gate from qubit 1 to qubit 2
qc.cp(pi / 4, 0, 2)  # Controlled-phase gate from qubit 0 to qubit 2
qc.draw(output="mpl")

qc.h(1)
qc.cp(pi / 2, 0, 1)  # Controlled-phase gate from qubit 0 to qubit 1

qc.draw(output="mpl")

qc.h(0)

qc.draw(output="mpl")

qc.swap(0, 2)

qc.draw(output="mpl")

例として $|5\rangle$ に QFT を適用してみます。

まず、整数 5 の 2 進数表現を確認し、状態 5 を符号化する回路を作成します：

bin(5)

'0b101'

qc = QuantumCircuit(3)

qc.x(0)
qc.x(2)
qc.draw(output="mpl")

Aer シミュレータを使って量子状態を確認します：

statevector = Statevector(qc)
plot_bloch_multivector(statevector)

最後に QFT を追加し、量子ビットの最終状態を確認します：

qc.h(2)
qc.cp(pi / 2, 1, 2)
qc.cp(pi / 4, 0, 2)

qc.h(1)
qc.cp(pi / 2, 0, 1)

qc.h(0)

qc.swap(0, 2)

qc.draw(output="mpl")

statevector = Statevector(qc)
plot_bloch_multivector(statevector)

QFT 関数が正しく動作していることが確認できます。量子ビット 0 は $\tfrac{5}{8}$ 回転、量子ビット 1 は $\tfrac{10}{8}$ 回転（$\tfrac{1}{4}$ 回転に相当）、量子ビット 2 は $\tfrac{20}{8}$ 回転（$\tfrac{1}{2}$ 回転に相当）しています。

2.5 演習：QFT

(1) 4量子ビットの QFT を実装します。

##your code goes here##

(2) $|14\rangle$ に QFT を適用し、Bloch 球を使って状態ベクトルをシミュレートしてプロットします。

##your code goes here##

演習の解答：QFT

（1）

qr = QuantumRegister(4)
cr = ClassicalRegister(4)
qc = QuantumCircuit(qr, cr)

qc.h(3)
qc.cp(pi / 2, 2, 3)
qc.cp(pi / 4, 1, 3)
qc.cp(pi / 8, 0, 3)

qc.h(2)
qc.cp(pi / 2, 1, 2)
qc.cp(pi / 4, 0, 2)

qc.h(1)
qc.cp(pi / 2, 0, 1)

qc.h(0)

qc.swap(0, 3)
qc.swap(1, 2)

qc.draw(output="mpl")

（2）

bin(14)

'0b1110'

qc = QuantumCircuit(4)

qc.x(1)
qc.x(2)
qc.x(3)
qc.draw("mpl")

qc.h(3)
qc.cp(pi / 2, 2, 3)
qc.cp(pi / 4, 1, 3)
qc.cp(pi / 8, 0, 3)

qc.h(2)
qc.cp(pi / 2, 1, 2)
qc.cp(pi / 4, 0, 2)

qc.h(1)
qc.cp(pi / 2, 0, 1)
qc.h(0)

qc.swap(0, 3)
qc.swap(1, 2)

qc.draw(output="mpl")

statevector = Statevector(qc)
plot_bloch_multivector(statevector)

3. 量子位相推定（QPE）の基礎

ユニタリ演算 $U$ が与えられたとき、QPE は $U\vert\psi \rangle =e^{\boldsymbol{2\pi i} \theta }|\psi \rangle$ における $\theta$ を推定します。$U$ はユニタリであるため、その固有値はすべてノルムが 1 です。

3.1 手順

1. 準備

$\vert\psi\rangle$ は一組の量子ビットレジスタに格納されます。追加の $n$ 量子ビットからなるカウンティングレジスタに、値 $2^n\theta$ を記録します：

$$|\psi_0\rangle = \lvert 0 \rangle^{\otimes n} \lvert \psi \rangle$$

2. 重ね合わせ

カウンティングレジスタに $n$ ビットのアダマールゲート演算 $H^{\otimes n}$ を適用します：

$$|\psi_1\rangle = {\frac {1}{2^{\frac {n}{2}}}}\left(|0\rangle +|1\rangle \right)^{\otimes n} \lvert \psi \rangle$$

3. 制御ユニタリ演算

対応する制御ビットが $|1\rangle$ のときにのみ、ターゲットレジスタにユニタリ演算子 $U$ を適用する制御ユニタリ $CU$ を導入する必要があります。$U$ は固有ベクトル $|\psi\rangle$ を持つユニタリ演算子であり、$U|\psi \rangle =e^{\boldsymbol{2\pi i} \theta }|\psi \rangle$ を満たすため、次のようになります：

$$U^{2^{j}}|\psi \rangle =U^{2^{j}-1}U|\psi \rangle =U^{2^{j}-1}e^{2\pi i\theta }|\psi \rangle =\cdots =e^{2\pi i2^{j}\theta }|\psi \rangle$$

3.2 例：T ゲートの QPE

QPE の例として $T$ ゲートを使い、その位相 $\theta$ を推定してみましょう。

$$T|1\rangle = \begin{pmatrix} 1 & 0\\ 0 & e^\frac{i\pi}{4}\\ \end{pmatrix} \begin{pmatrix} 0\\ 1\\ \end{pmatrix} = e^\frac{i\pi}{4}|1\rangle$$

期待される結果は次のとおりです：

$$\theta = \frac{1}{8}$$

ここで

$$T|1\rangle = e^{2i\pi\theta}|1\rangle$$

$X$ ゲートを適用して、$T$ ゲートの固有ベクトル $\vert\psi\rangle = \vert1\rangle$ を初期化します：

qpe = QuantumCircuit(4, 3)
qpe.x(3)
qpe.draw(output="mpl")

次に、カウンティング量子ビットにアダマールゲートを適用します：

for qubit in range(3):
    qpe.h(qubit)
qpe.draw(output="mpl")

制御ユニタリ演算を実行します：

repetitions = 1
for counting_qubit in range(3):
    for i in range(repetitions):
        qpe.cp(pi / 4, counting_qubit, 3)  # This is C-U
    repetitions *= 2
qpe.draw(output="mpl")

逆量子フーリエ変換を適用してカウンティングレジスタの状態を変換し、その後カウンティングレジスタを測定します：

from qiskit.circuit.library import QFT

# Apply inverse QFT
qpe.append(QFT(3, inverse=True), [0, 1, 2])
qpe.draw(output="mpl")

for n in range(3):
    qpe.measure(n, n)
qpe.draw(output="mpl")

Aer シミュレーターを使ってシミュレーションを実行できます：

aer_sim = AerSimulator()
shots = 2048

pm = generate_preset_pass_manager(backend=aer_sim, optimization_level=1)
t_qpe = pm.run(qpe)

sampler = Sampler(mode=aer_sim)
job = sampler.run([t_qpe], shots=shots)
result = job.result()
answer = result[0].data.c.get_counts()

plot_histogram(answer)

結果として 001 が確実に得られており、これは10進数で 1 に相当します。次に、$\theta$ を求めるために結果（1）を $2^n$ で割ります：

$$\theta = \frac{1}{2^3} = \frac{1}{8}$$

ショアのアルゴリズムでは、$\theta$ を未知としてその値を回路から取得することで、数の素因数分解が可能になります。

3.3 演習

3 量子ビットのカウンティングと固有ベクトル用の 1 量子ビットを使って、$\theta = 1/3$ を推定してください。

$P|1\rangle = e^{i\lambda}|1\rangle$。$U|1\rangle = e^{2\pi i \tfrac{1}{3}}|1\rangle$ を実装したいので、$\lambda = \tfrac{2 \pi}{3}$ と設定する必要があります。

##your code goes here##

演習の解答：$\theta = 1/3$

# Create and set up circuit
qpe = QuantumCircuit(4, 3)

# Apply H-Gates to counting qubits:
for qubit in range(3):
    qpe.h(qubit)

# Prepare our eigenstate |psi>:
qpe.x(3)

# Do the controlled-U operations:
angle = 2 * pi / 3
repetitions = 1
for counting_qubit in range(3):
    for i in range(repetitions):
        qpe.cp(angle, counting_qubit, 3)
    repetitions *= 2

# Do the inverse QFT:
qpe.append(QFT(3, inverse=True), [0, 1, 2])

for n in range(3):
    qpe.measure(n, n)

qpe.draw(output="mpl")

aer_sim = AerSimulator()
shots = 4096

pm = generate_preset_pass_manager(backend=aer_sim, optimization_level=1)
t_qpe = pm.run(qpe)

sampler = Sampler(mode=aer_sim)
job = sampler.run([t_qpe], shots=shots)
result = job.result()
answer = result[0].data.c.get_counts()

plot_histogram(answer)

4. IBM Quantum の Sampler プリミティブを使用した実行

実際の量子デバイスを使って QPE を実行し、Qiskit パターンの 4 ステップに従います。

1. 問題を量子ネイティブな形式にマッピングする

2. 回路を最適化する

3. ターゲット回路を実行する

4. 結果を後処理する

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Sampler

# Loading your IBM Quantum® account(s)

service = QiskitRuntimeService()

4.1 ステップ 1：問題を量子回路と演算子にマッピングする

qpe = QuantumCircuit(4, 3)
qpe.x(3)
for qubit in range(3):
    qpe.h(qubit)
repetitions = 1
for counting_qubit in range(3):
    for i in range(repetitions):
        qpe.cp(pi / 4, counting_qubit, 3)  # This is C-U
    repetitions *= 2
qpe.append(QFT(3, inverse=True), [0, 1, 2])
for n in range(3):
    qpe.measure(n, n)
qpe.draw(output="mpl")

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

print(backend)

<IBMBackend('ibm_strasbourg')>

4.2 ステップ 2：ターゲットハードウェア向けに最適化する

# Transpile the circuit into basis gates executable on the hardware
pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
qc_compiled = pm.run(qpe)

qc_compiled.draw("mpl", idle_wires=False)

4.3 ステップ 3：ターゲットハードウェアで実行する

real_sampler = Sampler(mode=backend)
job = real_sampler.run([qc_compiled], shots=1024)
job_id = job.job_id()
print("job id:", job_id)

job id: d13p4zb5z6q00087ag00

job = service.job(job_id)  # Input your job-id between the quotations
job.status()

'DONE'

result_real = job.result()
print(result_real)

PrimitiveResult([SamplerPubResult(data=DataBin(c=BitArray(<shape=(), num_shots=1024, num_bits=3>)), metadata={'circuit_metadata': {}})], metadata={'execution': {'execution_spans': ExecutionSpans([DoubleSliceSpan(<start='2025-06-09 22:39:00', stop='2025-06-09 22:39:00', size=1024>)])}, 'version': 2})

4.4 ステップ 4：結果を後処理する

from qiskit.visualization import plot_histogram

plot_histogram(result_real[0].data.c.get_counts())

# See the version of Qiskit
import qiskit

qiskit.__version__

'2.0.2'