$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$
はじめに
前回の記事で、任意の論理関数に対する量子回路が設計できるようになったので、量子オラクルを使った量子アルゴリズムをちょっと試してみたくなりました。今回、その代表選手というか基本の基である「Groverのアルゴリズム」を取り上げて、自作の量子計算シミュレータqlazyで、アルゴリズムの動作を確認します。参考にさせていただいたのは、以下の書籍・記事です。
- 中山茂「量子アルゴリズム」
- Groverの検索アルゴリズムとQuantum Amplitude Amplification/Estimation
- Groverアルゴリズムについて(@kyamaz)
- グローバーのアルゴリズムの量子回路を組む(@shiibass)
- Grover(グローバー)のアルゴリズム(@YuichiroMinato)
問題設定
Groverのアルゴリズムは、データベース検索とかファイル検索のアルゴリズムと呼ばれたりしますが、現在使われているようなデータベースの検索を効率化するものと思ってはいけません。あるクエリ$x$に対して{$0,1$}の値をとる$f(x)$を返すシステムがあったとき、$1$を返す$x$を見つける量子アルゴリズムです。例えば、
\ket{x} \otimes {\ket{y}} \rightarrow \ket{x} \otimes {\ket{y \oplus f(x)}} \tag{0}
のような入出力をする量子オラクルにおいて、なるべく少ない手数で$f(x)=1$となる$x$を見つけるにはどうすれば良いか、というのが基本的な問題設定です。
アルゴリズムの説明
量子振幅増幅
このような問題を考える際の常套手段は、あらゆる可能な解からなる重ね合わせ状態を用意しておいて、最終的に解に相当する状態の確率振幅だけが大きくなるような状態変化を、アルゴリズムとして実現させるというものです。具体的には、
\ket{\psi} = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \ket{k} \tag{1}
という量子状態の均等重ね合わせの中から、$f(x)=1$となる量子状態$\ket{x}$だけを突出させるアルゴリズムを考えれば良いです。そういう意味で、Groverのアルゴリズムは、特定の量子状態の振幅のみを増幅させる「量子振幅増幅」という形に、いまでは一般化されていて、もともとのGroverはその特別なケースとみなされているようです。
というわけで、一般的な議論からはじめます。ある正の整数(2進数)$x$に対して、$0$か$1$のどちらかを取る関数$f(x)$を考えます。$f(x)=0$を満たす$\ket{x}$によって張られる部分空間と、$f(x)=1$を満たす$\ket{x}$によって張られる部分空間は互いに直交していますので、任意の状態$\ket{\psi}$は、それぞれの部分空間(の要素)$\ket{\psi_{0}},\ket{\psi_{1}}$を用いて、
\ket{\psi} = \cos \theta \ket{\psi_{0}} + e^{i\phi} \sin \theta \ket{\psi_{1}} \tag{2}
と書けます。ここで、$\cos^2 \theta$は、量子状態$\ket{\psi}$を測定したときに$f(x)=0$となる$\ket{x}$を得る確率、$\sin^2 \theta$は、$f(x)=1$となる$\ket{x}$を得る確率です。これから、$f(x)=1$となる$\ket{x}$を得る確率をどうやって大きくしていけるかを考えていくので、最初なるべく小さい$\sin^2 \theta$をイメージしておきましょう。以下のように各部分空間を1次元に単純化して図示します(参考記事より)。
$\ket{\psi}$なるべく$\ket{\psi_{1}}$に近づけるために、以下のステップを実施します。
- [STEP.0] $\ket{\psi}$を準備します。
- [STEP.1] $\ket{\psi}$を$\ket{\psi_{0}}$軸に関して反転します($S_{f} \ket{\psi}$)。
- [STEP.2] $S_{f} \ket{\psi}$を元の状態$\ket{\psi}$に関して反転します($-S_{\psi} S_{f} \ket{\psi}$)。
これで、図からも明らかなように、元の$\ket{\psi}$は、$\ket{\psi_{1}}$軸に少し近づきます。これを繰り返すことで、さらに$\ket{\psi_{1}}$軸に近づけることができます。というわけで、
- [STEP.3] STEP.1とSTEP.2を適切な回数繰り返します。
です。量子振幅増幅の基本的な考え方は以上です。
では、$S_{f},S_{\psi}$はどんな変換なのかを考えてみます。まず、$S_{f}$は、
\begin{align}
S_f \ket{\psi_0} &= \ket{\psi_0} \\
S_f \ket{\psi_1} &= -\ket{\psi_1} \tag{3}
\end{align}
を満たす変換なので、
S_f = \ket{\psi_0} \bra{\psi_0} - \ket{\psi_1} \bra{\psi_1} \tag{4}
と書けそうです。また、$S_{\psi}$は、
-S_{\psi} = 2\ket{\psi} \bra{\psi} - I \tag{5}
です。いきなり天下ってしまいましたが、わかりますでしょうか？$\ket{\psi}\bra{\psi}$が$\ket{\psi}$への射影演算子であることに注意しましょう。また、状態は2次元空間上でのベクトルだと思うことにしましょう(実際には多次元ですが)。そうすると、式(5)の右辺は$\ket{\psi} \bra{\psi} + (\ket{\psi} \bra{\psi} - I)$と分解できます。第1項は$\ket{\psi}$への射影で、第2項以降の括弧はその射影と変換前の状態(ベクトル)そのものとの差分ベクトルです。その和なので、$\ket{\psi}$に関する反転演算子になります。
この式(5)は、もう少しわかりやすいゲートに分解できます。$\ket{\psi}$は$\ket{0}$にあるユニタリ演算子$U$を適用して生成するので、以下のように書けます。
\begin{align}
-S_{\psi} &= 2U\ket{0} \bra{0} U^{\dagger} - I \\
&= -U(I-2\ket{0}\bra{0})U^{\dagger} \\
&= -U S_0 U^{\dagger} \tag{6}
\end{align}
ここで、$S_0 = I - 2\ket{0} \bra{0}$とおきましたが、これは$\ket{0}$の符号のみを反転する演算子です。どうでしょう。こちらもわかりますでしょうか？
I = \sum_{k=0}^{N-1} \ket{k} \bra{k}
とおいてみればわかりますね。
この[STEP.1][STEP.2]を繰り返すのですが、何回繰り返せば良いかを、次に考えます。そのため、この[STEP.1][STEP.2]のプロセスを$Q=-S_{\psi} S_f$とおいて$\ket{\psi}$に対して$n$回繰り返した、$Q^n \ket{\psi}$を計算してみます。このままでは見通しが悪いので、$Q$の固有ベクトルで表現し直して計算するのが良いです。上での議論から、$Q$は、以下のように書き直すことができます。
\begin{align}
Q &= -S_{\psi} S_f \\
&= (2 \ket{\psi} \bra{\psi} - I) (\ket{\psi_0} \bra{\psi_0} - \ket{\psi_1} \bra{\psi_1}) \\
&= (2\alpha^2 -1) \ket{\psi_0}\bra{\psi_0} + 2\alpha\beta e^{i\psi}\ket{\psi_1}\bra{\psi_0} - 2\alpha\beta e^{-i\psi} \ket{\psi_0}\bra{\psi_1} - (2\beta^2 -1) \ket{\psi_1}\bra{\psi_1} \tag{7}
\end{align}
$Q$の固有ベクトルと固有値は、天下りですが、以下のようになります($Q$を適用してみれば、確認できますが、計算は省略します)。
固有ベクトル：
\ket{\psi_{\pm}} = \frac{1}{\sqrt{2}} (\ket{\psi_0} \mp \ket{\psi_1}) \tag{8}
固有値：
\lambda_{\pm} = e^{\pm i 2 \theta} \tag{9}
$\ket{\psi}$を固有ベクトルを基底にして表現し直すと、
\begin{align}
\ket{\psi} &= \ket{\psi_{+}}\braket{\psi_{+}}{\psi} + \ket{\psi_{-}}\braket{\psi_{-}}{\psi} \\
&= \frac{1}{\sqrt{2}}(\alpha + i\beta)\ket{\psi_{+}} + \frac{1}{\sqrt{2}}(\alpha - i\beta)\ket{\psi_{-}} \\
&= \frac{1}{\sqrt{2}} (e^{i\theta}\ket{\psi_{+}} + e^{-i\theta}\ket{\psi_{-}}) \tag{10}
\end{align}
となりますので、これに$Q^n$をかけます。
\begin{align}
Q^{n}\ket{\psi} &= \frac{1}{\sqrt{2}}(e^{i\theta} Q^n \ket{\psi_{+}} + e^{-i\theta} Q^n \ket{\psi_{-}}) \\
&= \frac{1}{\sqrt{2}}(e^{i(2n+1)\theta}\ket{\psi_{+}} + e^{-i(2n+1)\theta}\ket{\psi_{-}}) \\
&= \cos(2n+1)\theta \ket{\psi_0} + e^{i\phi} \sin(2n+1)\theta \ket{\psi_1} \tag{11}
\end{align}
ということで、$n$が大きくなるに従い、$\ket{\psi}$が$\ket{\psi_1}$軸に近づいていく様子が数式から見て取れるかと思います。が、無限に大きくすることで漸近的に近づくわけではなく、適切な$n$の値が存在します。式から明らかですが、$(2n+1)\theta$が$\pi/2$になるべく近くなるような$n$が適切な繰り返し回数です。
Groverのアルゴリズム
では、この一般的な議論を具体的な量子回路に落とし込んでみましょう。先程、式(2)のように最初に用意する状態を一般化して記載しましたが、式(1)のように均一重ね合わせ状態を用意しておくのが簡単です。再掲します。
\ket{\psi} = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \ket{k} \tag{1}
$f(x)=1$を満たす$x$の個数を$s$とすると、式(2)は、
\begin{align}
\ket{\psi} &= \cos \theta \ket{\psi_{0}} + e^{i\phi} \sin \theta \ket{\psi_{1}} \\
&= \sqrt{\frac{N-s}{N} }\ket{\psi_{0}} + \sqrt{\frac{s}{N}} e^{i\phi} \ket{\psi_{1}} \tag{12}
\end{align}
と書き直せます。後で使うので覚えておいて、先に進みます。先程のSTEP順に説明していきます。
[STEP.0] 初期状態の準備
量子ビット$n$個($2^n=N$)のレジスタを用意して、各々にアダマールをかければ、均等重ね合わせ状態が実現できます。
[STEP.1] マーキング
振幅増幅したい状態のみマイナスに反転するという意味で、このステップのことを「マーキング」と呼ぶようです。これを具体的にどうするかですが、答えを先に言います。式(0)を実現するような量子オラクルを召喚します。ただし、入力の制御レジスタは$n$量子ビットで、[STEP.0]で用意した均一重ね合わせを入れます。また、標的レジスタの$\ket{y}$には、$\ket{1}$にアダマールをかけたものを入れます(つまり、$(\ket{0}-\ket{1})/\sqrt{2}$)を入れます。
そうすると、どうなるかというと(面倒なので規格化定数は省略します)、
\begin{align}
&\sum_{x=0}^{N-1} \ket{x} \otimes (\ket{0}-\ket{1}) \\
&\rightarrow \sum_{x=0}^{N-1} \ket{x} \otimes (\ket{f(x)-\overline{f(x)}}) \tag{13}
\end{align}
となります。$f(x)-\overline{f(x)}$は、$f(x)$が0か1かによって符号が反転しますので、式(13)を$\ket{x}$に関する線形和とみなすと、その係数の符号によって、$f(x)$が0か1かを識別できる形にすることができるというわけです。
いままでの過程を量子回路で表してみます。
|0> -----H--| |--
|0> -----H--| |--
... |Uf|
|0> -----H--| |--
| |
|0> --X--H--| |--
ここで、Ufと書かれているのが、量子オラクルです。
[STEP.2] 振幅増幅反転
先程の図で示したように、この[STEP.2]の反転によって振幅が少し増幅するので、「振幅増幅反転」と呼ばれているようです。
式(6)で説明したように、$\ket{0}$のときだけ反転する演算子を、$\ket{\psi}$を用意したときに使ったユニタリ演算子でユニタリ変換します。$\ket{\psi}$を用意したときに使ったユニタリ演算子というのは、今の場合、「全ビットにアダマールをかける」です。では、$\ket{0}$のときだけ反転する演算はどうやって実現すれば良いでしょうか。答えは以下です。
--X--*--X--
--X--*--X--
...
--X--*--X--
--X--Z--X--
わかりますでしょうか？全体の$n$量子ビットのレジスタを$n-1$の制御レジスタと$1$ビットの標的レジスタに分けて、Multi-Controlled Zゲートを置いて、両側からXゲートではさみます。こうすると、$\ket{0}$(正確には$\ket{0}\ket{0}\cdots\ket{0}$)以外の場合、何も起きません。$\ket{0}$の場合のみ、制御がはたらいてZゲートがオンになり、かつ、Zへの入力は$\ket{1}$なので、マイナス符号が発生します。そして、Multi-Controlled Zの直後の状態は$-\ket{1}\ket{1}\cdots\ket{1}$となり、その後全体にXゲートをかけると、元の状態にマイナスがついたものになります。
ここまでの過程を量子回路で表すと、以下のようになります。
|0> -----H--| |--H--X--*--X--H--
|0> -----H--| |--H--X--*--X--H--
... |Uf| ...
|0> -----H--| |--H--X--Z--X--H--
| |
|0> --X--H--| |-----------------
[STEP.3] 繰り返し
これを繰り返します。つまり、
|0> -----H--| |--H--X--*--X--H--| |--H--X--*--X--H-- ...
|0> -----H--| |--H--X--*--X--H--| |--H--X--*--X--H-- ...
... |Uf| ... |Uf| ...
|0> -----H--| |--H--X--Z--X--H--| |--H--X--Z--X--H-- ...
| | | |
|0> --X--H--| |-----------------| |----------------- ...
です。このとき、$(2n+1)\theta$が$\pi/2$になるべく近くなるような$n$が適切な繰り返し回数でした。式(12)より、
\sin \theta = \sqrt{\frac{s}{N}} \tag{14}
としましたので、適切な繰り返し回数は、
\begin{align}
n &\sim \frac{\pi}{4\theta} - \frac{1}{2} \\
&= \frac{\pi}{4(\sin^{-1} \sqrt{\frac{s}{N}})} - \frac{1}{2} \tag{15}
\end{align}
となるような整数$n$です。もう少し言うと、右辺の小数を四捨五入すると考えれば、0.5を足して切り捨てすれば良いので、
\begin{align}
n &= \lfloor \frac{\pi}{4\theta} \rfloor \\
&= \lfloor \frac{\pi}{4(\sin^{-1} \sqrt{\frac{s}{N}})} \rfloor \tag{16}
\end{align}
となります。
シミュレータで確認
それでは、シミュレータでアルゴリズムの動作を確認します。今回、$f(x)$の$x$は10ビットを想定してやってみます。つまり、$x$は$0$から$1023$までの値を取ります。で、$f(x)=1$となる$x$は(何でも良かったのですが適当な双子素数にして)、881と883の2つとしてみました。なので、最終的に881と883の状態が大きくなって、測定に引っかかれば成功というわけです。
実装
全体のPythonコードは以下です。
【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。
import math
from qlazypy import QState
# custom gates
def hadamard(self,id):
for i in range(len(id)):
self.h(id[i])
return self
def flip(self,id):
for i in range(len(id)):
self.x(id[i])
return self
def multi_cz(self,id_in,id_out):
bitnum = len(id_in)
psi = 1.0/(2**(bitnum-1)) # unit=pi(radian)
gray_pre = 0
for gray in gray_code(bitnum):
if gray == 0:
continue
msb = len(str(bin(gray)))-3
chb = len(str(bin(gray^gray_pre)))-3
if gray != 1:
if chb == msb:
chb -= 1
self.cx(id_in[chb],id_in[msb])
self.cp(id_in[msb],id_out[0],phase=psi)
psi = -psi
gray_pre = gray
return self
def multi_cx(self,id_in,id_out):
self.h(id_out[0])
self.multi_cz(id_in,id_out)
self.h(id_out[0])
return self
def mpmct(self,id_in,id_out,x):
bitnum = len(id_in)
for i in range(bitnum):
bit = (x>>i)%2
if bit == 0:
self.x(id_in[i])
self.multi_cx(id_in,id_out)
for i in range(bitnum):
bit = (x>>i)%2
if bit == 0:
self.x(id_in[i])
return self
def q_oracle(self,id_in,id_out,func):
num = 2**len(id_in)
for x in range(num):
y = func(x)
if y == 1:
self.mpmct(id_in,id_out,x)
return self
# functions
def gray_code(n):
for k in range(2**n):
yield k^(k>>1)
def func(x):
if x == 881 or x == 883:
y = 1
else:
y = 0
return y
def create_register(num):
return [0]*num
def init_register(*args):
idx = 0
for i in range(len(args)):
for j in range(len(args[i])):
args[i][j] = idx
idx += 1
return idx
def result(self,id,solnum):
# measurement
iid = id[::-1]
freq = self.m(id=iid).frq
res = []
for _ in range(solnum):
idx = freq.index(max(freq))
res.append(idx)
freq[idx] = 0
return res
def opt_iter(N,solnum):
return int(0.25*math.pi/math.asin(math.sqrt(solnum/N)))
if __name__ == '__main__':
# add custom gates
QState.hadamard = hadamard
QState.flip = flip
QState.multi_cz = multi_cz
QState.multi_cx = multi_cx
QState.mpmct = mpmct
QState.q_oracle = q_oracle
QState.result = result
# set parameters
bits = 10
N = 2**bits
solnum = 2
# set register
id_in = create_register(bits)
id_out = create_register(1)
qnum = init_register(id_in,id_out)
id_in_c = id_in[:-1]
id_in_t = [id_in[-1]]
# initial state
qs = QState(qnum)
qs.hadamard(id_in)
qs.x(id_out[0]).h(id_out[0])
iter = opt_iter(N,solnum)
print("== algorithm start (total iter:{}) ==".format(iter))
for i in range(iter):
print("- iter no. = ", i)
qs.q_oracle(id_in,id_out,func)
qs.hadamard(id_in)
qs.flip(id_in)
qs.multi_cz(id_in_c,id_in_t)
qs.flip(id_in)
qs.hadamard(id_in)
# print results
res = qs.result(id_in,solnum)
print("== result==")
print("x =", res)
qs.free()
何をやっているか、簡単に説明します。量子オラクルの部分はほぼ前回の記事と同じものですので、アルゴリズムの実質的な処理部分のみの説明になります。プログラムのmain部の、
# initial state
qs = QState(qnum)
qs.hadamard(id_in)
qs.x(id_out[0]).h(id_out[0])
で、初期状態を用意しています。
iter = opt_iter(N,solnum)
print("== algorithm start (total iter:{}) ==".format(iter))
で、繰り返しの回数を式(16)に基づき計算しています(詳細はopt_iter関数の定義を参照)。
for i in range(iter):
print("- iter no. = ", i)
qs.q_oracle(id_in,id_out,func)
qs.hadamard(id_in)
qs.flip(id_in)
qs.multi_cz(id_in_c,id_in_t)
qs.flip(id_in)
qs.hadamard(id_in)
で、$Q=-S_{\psi} S_f$の演算を繰り返しています。hadamardメソッドは引数のレジスタの全量子ビットに対してアダマールを実行します。flipメソッドは引数の全量子ビットに対してビット反転(Xゲート)を実行します。multi_czはMulti-Controlled Zゲートです。というわけで、先程のアルゴリズムの説明通りに量子ゲートをならべているだけです。
最後に、
# print results
res = qs.result(id_in,solnum)
print("== result==")
print("x =", res)
で、測定結果を表示しています。
結果
では、このプログラムを実行した結果を示します。
== algorithm start (total iter:17) ==
- iter no. = 0
- iter no. = 1
- iter no. = 2
- iter no. = 3
- iter no. = 4
- iter no. = 5
- iter no. = 6
- iter no. = 7
- iter no. = 8
- iter no. = 9
- iter no. = 10
- iter no. = 11
- iter no. = 12
- iter no. = 13
- iter no. = 14
- iter no. = 15
- iter no. = 16
== result==
x = [883, 881]
17回繰り返して、無事881と883が得られました。めでたしめでたし。
おわりに
10量子ビットで繰り返し回数17回程度のGroverであれば、シミュレータで難なく計算することができました。実際の量子コンピュータ(例えば、IBM Qとか)では、そうはうまく行かないはずなので、このあたりの感覚も掴むべく、ちょっと実物も試してみたい気がしてきました。qlazyに機能追加してみっかなー、とか(でも、どうやる？)。このあたり、考えるのも自作の良さ＆楽しみですね(労力かけて車輪をつくっているだけとも言えるが、、まー、勉強にはなります)。
以上