グローバーのアルゴリズムの量子回路を組む

こんにちは、LIFULLのshiibassです。
LIFULL Advent Calendar7日目の記事になります。

皆さんはシンギュラリティという言葉をご存知でしょうか。技術レベルは指数関数的に成長していてある時を境に急激に成長を迎え、その点のことをシンギュラリティと呼ぶそうです。
弊社では、そういった未来でどうすればあらゆるLIFEをFULLにできるか、そのために今何をすべきか、などを議論しています。
シンギュラリティを起こす一要素が量子コンピュータで、これによって今までは現実的に計算不可能だった計算も可能になると言われていて、夢がある話ではありますが機械学習エンジニアやデータサイエンティストにとっては技術のキャッチアップで忙しくなるかもしれません。一方でもともと勉強していた人にとっては量子コンピュータのクラウド化が進んでいて将来が待ち遠しいことでしょう。

ということで今回は量子回路を組んでみます。量子コンピュータは無料で利用できるIBM Qを使います。

組む回路はグローバーのアルゴリズムとします。グローバーのアルゴリズムとは、複数の値の中から特定の値だけを抽出するというもので、例えば0,1,2,...,N-1と番号が付けられたデータがあったとして特定の値xを取り出すというものです。for文とif文を使えば実装できるのですが、最悪ケースだとN回のif文判定をする必要があります。一方で量子アルゴリズムを使えばずっと少ない処理数で計算することができるのでデータベースの検索システムに応用可能なのではないかと思います。
本記事ではN=4(2量子ビット)を想定し、$|11⟩$（3番のデータ）を抽出します。記事の中で量子プログラミングで使う演算が出てきますが、量子コンピュータの演算紹介を参考にしてください。

IBM QはGUIで回路を組むこともできるのですが、将来的に巨大な回路を組む時に作業しにくくなるのでQASMファイルの実装、APIでの実装を紹介します。

グローバーのアルゴリズムの概要

はじめにグローバーのアルゴリズムの手順をまとめておきます。
1. 量子重ね合わせ状態を作る
2. 抽出したい値xだけの位相を反転させる
3. 拡散変換でxの出現確率を高め、x以外の出現確率を低くする
4. 2,3を繰り返す
5. 測定する

QASMファイルでの実装

まずQASMファイルで作成する方法を記します。
QASMファイルはプログラミング言語よりも直感的に理解できる言語で問題を記述したファイルです。最終的には以下の図のようになり左側がQASMファイル、右側が生成された回路です。

上から順にコードを読んでいくと、

qreg q[2];
creg c[2];

qreg,cregはそれぞれ入力、出力回路を表します。今回は2量子ビットを定義しています。

1. 量子重ね合わせ状態を作る

h q[0];
h q[1];

hは量子ゲートのHを表していて後に続くビットに対して作用させることを意味しています。GUIで操作する場合は以下の図のように右側からゲートを選んで作用させたいビット番号にドラッグ＆ドロップすればできます。

h q[0];を数式で書くと、

\frac{1}{\sqrt{2}}\left(
\begin{matrix}
1 & 1 \\
1 & -1\\ 
\end{matrix}
\right)\left(
\begin{matrix}
1 \\
0\\ 
\end{matrix}
\right)
= \frac{1}{\sqrt{2}}\left(
\begin{matrix}
1 \\
1 \\ 
\end{matrix}
\right)

となり、重ね合わせの状態が作られます。この状態で出力するとそれぞれのビットで、0.5の確率で0が、0.5の確率で1が返ります（それぞれ成分の上、下の要素の2乗が0,1になる確率）。

2. 抽出したい値xだけの位相を反転させる

h q[1];
cx q[0], q[1];
h q[1];

は制御NOTゲートを拡張した制御Zゲートです。制御Zゲートは|11⟩の位相を反転させるゲートで第2ビットだけの動きを計算すると、

HXH = \frac{1}{\sqrt{2}}\left(
\begin{matrix}
1 & 1 \\
1 & -1\\ 
\end{matrix}
\right)
\left(
\begin{matrix}
0 & 1 \\
1 & 0\\ 
\end{matrix}
\right)
\frac{1}{\sqrt{2}}\left(
\begin{matrix}
1 & 1 \\
1 & -1\\ 
\end{matrix}
\right)
=\left(
\begin{matrix}
1 & 0 \\
0 & -1\\ 
\end{matrix}
\right)

このようになり、第2ビットが1のときに反転することがわかります。第1ビットは制御ビットと呼ばれる判定用のビットなので変化はしません。全体としては$C_Z$の成分は以下の通りです。

C_Z = \left(
\begin{matrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & -1 
\end{matrix}
\right)

これによって|11⟩(3番のデータ)だけ位相反転するので他の値との区別が可能になります。確かめてみると、

C_Z (|11⟩)= \left(
\begin{matrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & -1 
\end{matrix}
\right)
\left(
\begin{matrix}
0\\
0\\
0\\
1 
\end{matrix}
\right)\\
=
\left(
\begin{matrix}
0\\
0\\
0\\
-1 
\end{matrix}
\right)
=-|11⟩

となり、|11⟩の位相が反転しました。ちなみに|00⟩,|01⟩,|10⟩(それぞれ0,1,2番のデータ)は変化しません。前工程では重ね合わせの状態を作ったので、ここでは同時に4パターンすべてを扱えたことになります。今回は|11⟩を探したい値にしましたが、別の値の場合は間に位相を動かすゲートを組み込むことで実現できます。

3. 拡散変換でxの出現確率を高め、x以外の出現確率を低くする

次の処理は以下の通りです。

h q[0];
h q[1];
x q[0];
x q[1];
h q[1];
cx q[0], q[1];
h q[1];
x q[0];
x q[1];
h q[0];
h q[1];

これは拡散変換と呼ばれています。やっていることとしては前工程でずらした位相の違いを利用して|11⟩の出現確率を高め、他の出現確率を低くするというものです。この回路を組み込む前の状態では重ね合わせと制御Zを作用させているので、出力させたとしたら

\frac{1}{2}|00⟩+\frac{1}{2}|01⟩+\frac{1}{2}|10⟩-\frac{1}{2}|11⟩

となっていて、それぞれ等確率となっています。ここに確率振幅(係数)の平均値$m = \frac{1}{4}$($\frac{1}{2}$が3つ,$-\frac{1}{2}$が1つの平均値)から自身の確率振幅を引いて、再度平均値を足すと以下のようになります。

(m - \frac{1}{2} + m)|00⟩+(m-\frac{1}{2}+m)|01⟩+(m-\frac{1}{2}+m)|10⟩+(m-(-\frac{1}{2})+m)|11⟩=|11⟩

このように|11⟩だけが確率1となり、探していた値を抽出することができました。

ここからは難しくなりますが次工程まで飛ばしても概要を理解する上では全く問題ありません。

拡散変換について理解を深めるためにn量子ビットの場合で考えます。n量子ビットの場合には$N=2^n$として拡散変換$D$は以下のように、対角成分が$-1+2/N$でそれ以外が$2/N$の行列で表せます。

D_n = 
\left(
\begin{matrix}
-1+2/N & 2/N & \cdots & 2/N\\
2/N & -1+2/N & \ddots & \vdots\\
\vdots & \ddots & \ddots & 2/N\\
2/N & \cdots & 2/N & -1+2/N 
\end{matrix}
\right)

一つだけ位相が反転した量子重ね合わせの状態に対してこの変換$D$を適用して

\left(
\begin{matrix}
-1+2/N & 2/N & \cdots & 2/N\\
2/N & -1+2/N & \ddots & \vdots\\
\vdots & \ddots & \ddots & 2/N\\
2/N & \cdots & 2/N & -1+2/N 
\end{matrix}
\right)
\frac{1}{\sqrt{N}}\left(
\begin{matrix}
1\\
\vdots\\
1\\
-1\\
1\\
\vdots\\
1 
\end{matrix}
\right)

を計算すると反転していない状態の要素は

(-1+2/N) + (2/N) \times (N-2) + (2/N) \times (-1) = \frac{1}{\sqrt{N}}(1-\frac{4}{N})<\frac{1}{\sqrt{N}}

反転した状態の要素は

(-1+2/N) \times (-1) + (2/N) \times (N-1) = \frac{1}{\sqrt{N}}(3-\frac{4}{N})>\frac{1}{\sqrt{N}}

となり、反転していない状態の出現確率は減り、反転した状態の出現確率が増えました。全体で確率が1のままであること、$N=4(n=2)$のときに2量子ビットのときと同じ結果になることも計算すると確かめられます。

4. 2,3を繰り返す

先ほどの位相反転と拡散変換によって抽出したい状態の確率が上がりましたが、それを繰り返すことで1に近づけることができます。

5. 測定する

measure q[0]->c[0];
measure q[1]->c[1];

は出力回路に結果を渡しています。ここまでを書くと冒頭に記した図にある回路が出来上がります。

実際に動作させてみますと以下のように|11⟩(3番のデータ)の確率が1となっています。すなわち、|11⟩のデータを抽出できたことになります。

APIでの実装

次にAPIを利用した実装を書きます。私はこれが一番使いやすいと思っています。

pip install qiskit

でIBM QのAPIで利用することができます。全体のコードを書くと以下のようになります。qiskitについては以下を参考にしました。
https://www.qiskit.org/documentation/

from qiskit import QuantumProgram

API_TOKEN = ""
URL = 'https://quantumexperience.ng.bluemix.net/api'
BACKEND = "ibmqx_qasm_simulator"

specs = {
    "circuits": [{
        "name": "demo",
        "quantum_registers": [{
            "name": "q",
            "size": 2
        }],
        "classical_registers": [
            {"name": "c",
             "size": 2},
        ]}]
}

qp = QuantumProgram(specs=specs)
qp.set_api(API_TOKEN, URL)
qc = qp.get_circuit("demo")

q = qp.get_quantum_register("q")
c = qp.get_classical_register("c")

qc.h(q[0])
qc.h(q[1])

qc.cz(q[0], q[1])

qc.h(q[0])
qc.h(q[1])
qc.x(q[0])
qc.x(q[1])
qc.cz(q[0], q[1])
qc.x(q[0])
qc.x(q[1])
qc.h(q[0])
qc.h(q[1])

qc.measure(q[0], c[0])
qc.measure(q[1], c[1])

result = qp.execute(["demo"], backend=BACKEND, shots=2048)
print(result.get_ran_qasm("demo"))
print("____________________")
print(result.get_counts("demo"))

全体的にQASMファイルで書いたものとほとんど同じなのですが、Pythonで書けば量子回路のパーツごとにメソッドを切り分けて独自にカスタマイズできる点で便利です（QASMファイルでもそのような記述方法があるのかどうかはわかりません）。
Pythonファイルを実行すると以下のようになります。

OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
creg c[2];
u2(0.0,3.141592653589793) q[1];
u2(0.0,3.141592653589793) q[1];
u2(0.0,3.141592653589793) q[0];
cx q[0],q[1];
u2(0.0,3.141592653589793) q[1];
u2(0.0,3.141592653589793) q[1];
u3(3.141592653589793,0.0,3.141592653589793) q[1];
u2(0.0,3.141592653589793) q[1];
u2(0.0,3.141592653589793) q[0];
u3(3.141592653589793,0.0,3.141592653589793) q[0];
cx q[0],q[1];
u2(0.0,3.141592653589793) q[1];
u3(3.141592653589793,0.0,3.141592653589793) q[1];
u2(0.0,3.141592653589793) q[1];
measure q[1] -> c[1];
u3(3.141592653589793,0.0,3.141592653589793) q[0];
u2(0.0,3.141592653589793) q[0];
measure q[0] -> c[0];

____________________
{'11': 2048}

print("____________________")の上部には量子回路の説明が、下部には出力結果が表示されていて、確かに|11⟩が抽出できています。量子回路の説明でu2やu3が出てきますが、これらは量子ゲートの一般形に$0$や$\pi$などを入力してHやXを作っているからだと思われます。

まとめ

今回は量子回路をQASMファイルとAPIでの2つの方法で組んでみました。2量子ビットの回路でもコードにすると複雑になりがちなので本格的に利用する場合はプログラミングでうまく処理できるように工夫する必要がありそうです。量子プログラミングを業務で活用した事例はまだあまり耳にしないですが、IBM以外でも量子コンピュータのクラウドサービスを提供する予定があることを発表している企業もあり、今後が楽しみな分野です。