量子位相推定

はじめに

　今回は量子位相推定と呼ばれる手法を解説する。量子フーリエ変換とともに、多くの量子アルゴリズム（暗号解読、量子化学、量子機械学習など）の要素技術として使われる手法である。

量子位相推定とは（概論）

　量子アルゴリズムとは、さまざまな量子回路を組み合わせて構築するものであることを過去のブログで見てきた。量子回路は量子ゲートから構成され、1量子ゲートと2量子ゲートがあることも見た。1量子ゲートは、1量子ビットからなる状態ベクトル $|a\rangle$ （2次元列ベクトル）に $2\times 2$ のユニタリー行列 $U$ （量子アルゴリズムで使われる演算は全てユニタリー性を持たなければならない）を施すものである。

(1) $\begin{equation*} |b\rangle=U|a\rangle \end{equation*}$

図で示すと以下のようになる。左側が入力ビット、右側が出力ビットである。

$U$ の固有値は一般に $e^{i2\pi\phi}$ （ $\phi$ は $0\leq\phi<1$ を満たす実数）と書くことができるので、対応する固有状態を $|\psi\rangle$ とすると

(2) $\begin{equation*} U|\psi\rangle=e^{i2\pi\phi}|\psi\rangle \end{equation*}$

が成り立つ。固有状態 $|\psi\rangle$ が既知のとき、位相 $\phi$ を求めるアルゴリズムが量子位相推定である。

量子位相推定とは（詳細）

　先のブログで量子フーリエ変換 $F$ を解説した。

(3) $\begin{equation*} F|j\rangle=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}e^{i\frac{2\pi kj}{N}}|k\rangle \end{equation*}$

ここで、 $N=2^n$ であり、 $n$ は量子ビット数である。式(3)を眺めると、左辺の状態ベクトル内にある $j$ が右辺の指数関数の肩に移動していることが分かる。つまり、逆フーリエ変換をすれば、指数関数の肩にある値を状態ベクトル内部に移すことができる。

(4) $\begin{equation*} F^{-1}\left(\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}e^{i\frac{2\pi kj}{N}}|k\rangle\right)=|j\rangle \end{equation*}$

この事実を踏まえて式(2)を見直す。いま求めたい値 $\phi$ は指数関数の肩に乗っているので、逆フーリエ変換を利用して状態ベクトル側に下ろす。その後この状態を測定することにより $\phi$ の値を知ることできるだろう。この方針に従って以下の計算を行う。

　まず最初に、 $n$ 個の量子ビットを用意する。これらをまとめて第1レジスタと呼ぶことにする（下図参照）。全ての量子ビットの状態を0に初期化する。第1レジスタとは別に $m$ 個の量子ビットから成る第2レジスタを用意する。第2レジスタ上には固有値を求めたいユニタリー行列 $U$ の固有状態 $|\psi\rangle$ が既に実現されているとする。

第1レジスタと第2レジスタを合わせた全状態を $|\Psi\rangle$ と書くことにすると

(5) $\begin{equation*} |\Psi\rangle=|0\rangle^{\otimes n}\;|\psi\rangle \end{equation*}$

となる。ここで、 $\otimes n$ は同じ状態を $n$ 回繰り返すことを意味する（ $n$ 個の直積）。次に、第1レジスタに属する個々の量子ビットに $H$ ゲート（アダマールゲート）を作用させる（下図参照）。

全状態は以下のようになる。

(6) $\begin{equation*} |\Psi\rangle=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}|k\rangle\;|\psi\rangle \end{equation*}$

ここで、 $N=2^n$ とした。式(6)の第1レジスタの状態がどのような状態なのかを、 $n=3$ の場合に調べてみる。

(7) $\begin{eqnarray*} \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}|k\rangle\ &=& \frac{1}{\sqrt{8}}\left( |0\rangle+|1\rangle+\cdots+|7\rangle \right) \\ &=& \frac{1}{\sqrt{8}}\left( |000\rangle+|001\rangle+\cdots+|111\rangle \right) \\ &=& \frac{1}{\sqrt{8}}\left( |0\rangle|0\rangle|0\rangle+|0\rangle|0\rangle|1\rangle+\cdots+|1\rangle|1\rangle|1\rangle \right) \end{eqnarray*}$

式中の左から右への量子ビットの並び $|q_1q_2q_3\rangle$ は、図における量子ビットの上から下への並びに対応している（下図参照）。

以上が3量子ビットのときの第1レジスタの内容である。この後の説明を一般の $n$ 量子ビットに戻って行うと分かり難くなるので、もうしばらく3量子ビットのまま説明を続ける。

　アダマールゲートの後、第1レジスタの各ビットを制御ビットとして第2レジスタに $U$ の累乗を作用させる（下図参照）。

制御ビットが１の時のみ $U^k$ を作用させると言う意味である。全状態は以下のようになる。

(8) $\begin{equation*} |\Psi\rangle=\frac{1}{\sqrt{8}}\left( |000\rangle\;|\psi\rangle + |001\rangle\;U|\psi\rangle + |010\rangle\;U^2|\psi\rangle + |011\rangle\;U^{1+2}|\psi\rangle +\cdots \right) \end{equation*}$

例えば、式(8)の右辺第2項を見てみる。第1レジスタの $q_3$ ビットは1であるから上の回路図より第2レジスタに $U$ がかかる。また、第4項は $q_2$ ビットと $q_3$ ビットが1である。上の回路図を見ると、 $q_2$ ビットと紐づいている演算は $U^2$ 、 $q_3$ ビットと紐づいている演算は $U^1$ であるからトータルで $U^3$ が第2レジスタに作用する。ここまでの考察で、第1レジスタの2進数を10進数に変換した値が $U$ の肩に乗ることが分かったので

(9) $\begin{equation*} |\Psi\rangle=\frac{1}{\sqrt{8}}\sum_{k=0}^{7}|k\rangle\;U^k|\psi\rangle \end{equation*}$

とまとめることができる。ここまでくれば一般の $n$ ビットの場合に拡張するのは容易である。

(10) $\begin{equation*} |\Psi\rangle=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}|k\rangle\;U^k|\psi\rangle \end{equation*}$

ここで、 $N=2^n$ である。いま

(11) $\begin{equation*} U|\psi\rangle=e^{i2\pi\phi}|\psi\rangle \end{equation*}$

が成り立つのでこれを代入すると

(12) $\begin{eqnarray*} |\Psi\rangle &=& \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}|k\rangle\;U^k|\psi\rangle \\ &=& \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}|k\rangle\;e^{i2\pi\phi k}|\psi\rangle \\ &=& \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}|k\rangle\;e^{i\frac{2\pi k (N\phi)}{N}}|\psi\rangle \end{eqnarray*}$

を得る。式(4)の逆フーリエ変換は以下であった。

(13) $\begin{equation*} F^{-1}\left(\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}e^{i\frac{2\pi kj}{N}}|k\rangle\right)=|j\rangle \end{equation*}$

この式に出てくる $j$ を $N\phi$ と思えば、第1レジスタに逆フーリエ変換をかけた後の状態は次式となる。

(14) $\begin{equation*} F^{-1}|\Psi\rangle=|N\phi\rangle|\psi\rangle \end{equation*}$

このときの量子回路は下図である（3量子ビットのとき）。

最後に第1レジスタを測定し $N$ で割れば $\phi$ を求めることができる。

　上の計算は $N\phi$ が整数となる場合に成り立つ結果である。 $N\phi$ が整数とならない場合でも、その値に近い状態が高い確率で観測される。その誤差は $2\pi$ を $N=2^n$ で割った程度であるから、ビット数 $n$ を増やすほど精度は良くなる。

実装

　ここまでのロジックを、IBMのSDKQiskitを用いて実装してみる（ソースコードはここにある。）。 $U$ として次を考える。

(15) $\begin{equation*} U=\left( \begin{array}{rr} 1 & 0 \\ 0 & e^{i\theta} \\ \end{array} \right) \end{equation*}$

ここで

(16) $\begin{equation*} |1\rangle=\left( \begin{array}{c} 0 \\ 1 \end{array} \right) \end{equation*}$

であることを思い出せば、次式が成り立つことが分かる。

(17) $\begin{equation*} U|1\rangle=e^{i\theta}|1\rangle \end{equation*}$

つまり、 $U$ の固有値は $e^{i\theta}$ 、固有状態 $|\psi\rangle$ は $|1\rangle$ である。式(11)と対応させれば、 $\phi=\theta/(2\pi)$ である。量子回路を作るコードは以下のようになる。

import math
import qiskit


# 第2レジスタのビット数
N_EIGEN_STATE = 1


def execute_IQFT_core(circuit, target, n):
    for control in range(0, target):
        circuit.cp(-math.pi / 2 ** (target - control), control, target)
    circuit.h(target)


def execute_IQFT(circuit, n):
    """
    逆フーリエ変換(IQFT)を実行
    """

    # ビットの並びを反転
    for i in range(math.floor(n / 2)):
        circuit.swap(i, n - (i + 1))

    for i in range(n):
        execute_IQFT_core(circuit, i, n)


def make_circuit(n_encode: int):
    """
    量子位相推定を実行する。

    Parameters
    -----
    n_encode: int
        第2レジスタのビット数

    Returns
    -----
    qc: qiskit.QuantumCircuit
        作成した回路
    theta: qiskit.circuit.Parameter
        後で設定する位相
    """

    n = n_encode + N_EIGEN_STATE
    qc = qiskit.QuantumCircuit(qiskit.QuantumRegister(n), qiskit.ClassicalRegister(n_encode))

    # 第2レジスタをUの固有状態にする。
    qc.x(n_encode)

    # 第1レジスターのそれぞれのビットにアダマールゲートを作用させる。
    for qubit in range(n_encode):
        qc.h(qubit)

    # 第1レジスターの各ビットを制御ビットにして第2レジスタにユニタリー操作をおこなる。
    theta = qiskit.circuit.Parameter('φ')
    r = 1
    for c in range(n_encode):
        for i in range(r):
            qc.cp(theta, control_qubit=c, target_qubit=n_encode)
        r *= 2

    qc.barrier()

    # 逆フーリエ変換
    execute_IQFT(qc, n_encode)

    qc.barrier()

    # 第1レジスタの各ビットを測定する。
    for n in range(n_encode):
        # n番目の量子ビットを測定してn番目の古典ビットに保存する。
        qc.measure(n, n)

    return qc, theta

import math

import qiskit

# 第2レジスタのビット数

N_EIGEN_STATE = 1

def execute_IQFT_core(circuit, target, n):

for control in range(0, target):

circuit.cp(-math.pi / 2 ** (target - control), control, target)

circuit.h(target)

def execute_IQFT(circuit, n):

"""

逆フーリエ変換(IQFT)を実行

"""

# ビットの並びを反転

for i in range(math.floor(n / 2)):

circuit.swap(i, n - (i + 1))

for i in range(n):

execute_IQFT_core(circuit, i, n)

def make_circuit(n_encode: int):

"""

量子位相推定を実行する。

Parameters

-----

n_encode: int

第2レジスタのビット数

Returns

-----

qc: qiskit.QuantumCircuit

作成した回路

theta: qiskit.circuit.Parameter

後で設定する位相

"""

n = n_encode + N_EIGEN_STATE

qc = qiskit.QuantumCircuit(qiskit.QuantumRegister(n), qiskit.ClassicalRegister(n_encode))

# 第2レジスタをUの固有状態にする。

qc.x(n_encode)

# 第1レジスターのそれぞれのビットにアダマールゲートを作用させる。

for qubit in range(n_encode):

qc.h(qubit)

# 第1レジスターの各ビットを制御ビットにして第2レジスタにユニタリー操作をおこなる。

theta = qiskit.circuit.Parameter('φ')

r = 1

for c in range(n_encode):

for i in range(r):

qc.cp(theta, control_qubit=c, target_qubit=n_encode)

r *= 2

qc.barrier()

# 逆フーリエ変換

execute_IQFT(qc, n_encode)

qc.barrier()

# 第1レジスタの各ビットを測定する。

for n in range(n_encode):

# n番目の量子ビットを測定してn番目の古典ビットに保存する。

qc.measure(n, n)

return qc, theta

この回路を実行するコードは以下の通り。

import quantum_phase_estimation as QPE
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt
import numpy as np
import qiskit

# 第1レジスタのビット数
N_ENCODE = 5

if __name__ == "__main__":
        
    # 回路を作る。
    qc, theta = QPE.make_circuit(N_ENCODE)

    # θ=2*pi*φを与える（正解値）
    phase = np.random.rand()

    # 回路にパラメータを設定する。
    qc_parametrized = qc.bind_parameters({theta: phase})
    
    # シミュレータを選択する。
    backend = qiskit.Aer.get_backend('qasm_simulator')

    # 実行回数
    shots = 1024

    # 計算する。
    results = qiskit.execute(qc_parametrized, backend=backend, shots=shots).result()
    
    # 結果をとりだす。
    answer = results.get_counts()

    # 測定された位相を算出する。
    values = list(results.get_counts().values())
    keys = list(results.get_counts().keys())
    idx = np.argmax(list(results.get_counts().values()))
    ans = int(keys[idx], 2)
    phase_estimated = ans / (2 ** N_ENCODE)

    # 正しい位相の値
    true_phase = phase / (2 * np.pi)

    print('True phase: {:.4f}'.format(true_phase))
    print('Estimated phase: {:.4f}'.format(phase_estimated))
    print('Diff: {:.4f}'.format(np.abs(true_phase - phase_estimated)))
    
    # ヒストグラムを描画して保存する。
    plot_histogram(answer, figsize=(20, 7))
    plt.savefig("./histogram.jpg")

import quantum_phase_estimation as QPE

from qiskit.visualization import plot_histogram

import matplotlib.pyplot as plt

import numpy as np

import qiskit

# 第1レジスタのビット数

N_ENCODE = 5

if __name__ == "__main__":

# 回路を作る。

qc, theta = QPE.make_circuit(N_ENCODE)

# θ=2*pi*φを与える（正解値）

phase = np.random.rand()

# 回路にパラメータを設定する。

qc_parametrized = qc.bind_parameters({theta: phase})

# シミュレータを選択する。

backend = qiskit.Aer.get_backend('qasm_simulator')

# 実行回数

shots = 1024

# 計算する。

results = qiskit.execute(qc_parametrized, backend=backend, shots=shots).result()

# 結果をとりだす。

answer = results.get_counts()

# 測定された位相を算出する。

values = list(results.get_counts().values())

keys = list(results.get_counts().keys())

idx = np.argmax(list(results.get_counts().values()))

ans = int(keys[idx], 2)

phase_estimated = ans / (2 ** N_ENCODE)

# 正しい位相の値

true_phase = phase / (2 * np.pi)

print('True phase: {:.4f}'.format(true_phase))

print('Estimated phase: {:.4f}'.format(phase_estimated))

print('Diff: {:.4f}'.format(np.abs(true_phase - phase_estimated)))

# ヒストグラムを描画して保存する。

plot_histogram(answer, figsize=(20, 7))

plt.savefig("./histogram.jpg")

8行目のN_ENCODEは第1レジスタのビット数である。16行目で $\theta(=2\pi\phi)$ の正解値を与える。今回は量子コンピュータの実機ではなくPC上で行うシミュレーターを選択した（22行目）。28行目で計算を実行し、38行目で位相 $\phi$ の推定結果を算出する。45行目で正解値との差分を計算している。実際に計算したときの出力は以下のようになる。

True phase: 0.1178
Estimated phase: 0.1250
Diff: 0.0072

True phase: 0.1178

Estimated phase: 0.1250

Diff: 0.0072

第1レジスタの様子は以下のようになる。

量子ビット数を5としたから横軸の最大値は11111(10進数で31(= $2^5-1$ ))である。縦軸は各ビットの実現確率を表す。このヒストグラムを見ると、状態4(=00100)の実現確率が最も高いことが分かる。従って推定された位相 $\phi$ の値は $0.125(=4/2^5)$ となる。正解値とのずれは0.0072である。

まとめ

　今回は、量子位相推定（QPE）のアルゴリズを解説した。さまざまな分野（量子化学・暗号読解など）での応用が考えられる要素技術である。しかし、QPEを精度よく実行するには、現在の小規模な量子コンピュータ（NISQ）では無理であり、誤り訂正のある大規模な量子コンピュータ（FTQC）の出現を待つ必要がある。その時期は2029年ごろと言われている（記事）。

参考文献

IBM Quantumで学ぶ量子コンピュータ

量子コンピューティング　基本アルゴリズムから量子機械学習まで

量子位相推定アルゴリズム (Quantum Phase Estimation) とは

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

男性育休インタビュー

ユーザーストーリーマッピング

量子位相推定

はじめに

量子位相推定とは（概論）

量子位相推定とは（詳細）

実装

まとめ

参考文献

Kumada Seiya

最近の記事

顔認証（DeepFace vs GPT4o）

時系列解析再訪

πの新しい公式

中心極限定理

DBSCANを用いたクラスタリング

レイアウト解析

量子テレポーテーション

イジングモデルと巡回セールスマン問題

量子フーリエ変換

量子ビットとは何か？

量子もつれ

量子振幅増幅

Google Vision APIでOCR

アーカイブ

カテゴリー