Groverのアルゴリズム(2019年 Week2)#

Note

このNotebookはQiskit v0.44の仕様に合わせてコードを改変しています。

ウィーク2へようこそ!今週は、量子アルゴリズムの一つであるGrover (グローバー) のアルゴリズムについて学びます。 そして、当チュートリアルの最後には今週の問題(ラーニングチャレンジ演習Ⅱ)がありますので、ぜひチャレンジしてみてください。

古典コンピューターを凌駕する量子コンピューターの数あるアドバンテージの1つに、データベース検索を高速に行えるというのを聞いたことがあるかも知れません。Groverのアルゴリズムはこの能力を実証します。Groverのアルゴリズムは、非構造化データの検索問題に対して二次のオーダーの高速化ができるだけではなく、検索を超えて利用することができます。つまり、その他様々のアルゴリズムの実行時間を二次のオーダーで改善する一般的なテクニック、もしくはサブルーチンとして利用することができます。これは振幅増幅テクニックと呼ばれています。

このページでは、検索問題の詳細を説明し、検索問題用のオラクル(Oracle)を回路で作成した後にGroverのアルゴリズムをQisKitを利用して実装します。

非構造化データの検索#

\(N\)個の大きなアイテムリストがあるとします。その中で、一つだけアタリ\(w\)があるとします。

アタリの箱(マークのついたアイテム)を見つけるためには、古典計算では平均で \(N/2\) 個の箱を探す必要があります。 最悪の場合は、\(N\) 個探す必要があります。ところが、量子コンピューターでは Groverの振幅増幅のテクニックを使って、 おおよそ \(\sqrt N\) ステップでマークされたアイテムを探し出すことができます。 二次の高速化は、大きなリスト内のマークされたアイテムを探すためには実際の所、大きな時間の節約になります。 さらに、このアルゴリズムはリスト自体の内部構造を利用しないので、一般化することができ、多くの古典の問題でも二次の速度向上をもたらしてくれます。

オラクル#

量子コンピューターにこのようなリストをどのように提供すればいいでしょうか?このようなリストをエンコードするためによく用いられる手法としては、マークのついていない\(x\)には\(f(x)=0\)を返し、アタリ\(w\)に対しては\(f(w)=1\)を返すような関数\(f\)を定義する手法が挙げられます。この問題に対して量子コンピューターを使用するには、この関数に対して、重ね合わせ状態にあるアイテムを提供する必要があります。つまり、この関数をオラクルと呼ばれるユニタリー行列としてエンコードする必要があります。そのために、まず初めにアイテムをバイナリ表現にします。すなわち、\(x,w \in \{0,1\}^n\)\(N=2^n\)です。このようにすると量子コンピューター上の量子ビットとして表現することができます。その後、オラクル行列\(U_w\)を任意の標準基底\(|x>\)に対して\(U_w |x> =(-1)^{f(x)}|x>\)となるように定義します。

\(x\)がマークされていないアイテムの場合、オラクルは状態に対して何もしません。ところが、状態\(|w>\)に対してオラクルを適用すると、\(U_w |w> = -|w>\)にマッピングします。幾何学的には、このユニタリー行列は\(N=2^n\)次元のベクトル空間において、マークされたアイテムを原点に対して反転する操作に対応しています。

振幅増幅(Amplitude amplification)#

では、アルゴリズムはどのように動作するのでしょう?リストを調べる前は、私たちはマークされたアイテムがどこにあるのか知りません。従って、私たちの推測は、下記の式で表される均一な重ね合わせ状態での位置特定と大差ありません。 $\( |s> = \frac{1}{\sqrt N}\sum_{x=0}^{N-1} |x> \)$

もしこの時点で標準基底 \(|x> \)でこの重ね合わせ状態を測定した場合、\(\frac{1}{N} = \frac{1}{2^{n}} \) の確率で、標準基底のうちの一つに収束します。予想通り、正しい\(|w>\)を当てる確率は\(\frac{1}{2^{n}}\) です。従って、正しいアイテムを推測するには、平均\(N=2^{n}\)回トライする必要があります。

そこで振幅増幅と呼ばれる処理を加えましょう。この処理により、量子コンピューターが正しいアイテムを見つける確率を大幅に高めることが出来ます。この処理では、マークされたアイテムの振幅を増幅し、その他のアイテムの振幅を小さくします。この結果、最終状態を測定すると、正しいアイテムをほぼ確実に取り出すことができるようになります。

このアルゴリズムには2つの反転という面白い幾何学的解釈があり、2次元平面での回転として表せます。私たちが考慮すべきは、アタリ\(|w>\) と均一な重ね合わせ状態\(|s>\) の2つの特別な状態のみです。この2つのベクトルは、ベクトル空間 \(\mathbb C^{N}\) において、2次元の平面を張ります。\(|w>\) 状態は、\(N^{-1/2}\) の振幅で重ね合わせ状態に入っているため、これら2つのベクトルは完全に直交しているわけではありません。

しかし、\(|s>\) から \(|w>\) を削除し、正規化し直す事で\(|w>\) に直交する追加の状態 \(|s'>\) を導入することができます。

Step 0:#

振幅増幅は均一な重ね合わせ状態 \(|s>\) から開始します。均一な重ね合わせ状態は、\(|s> = H^{\otimes n}|0>^{n}\) により簡単に作成できます。\(t=0\) の時、初期状態は \(|\psi_{0}> = |s>\) です。

左の図は、\(|w>, |s'>\) によって張られる、2次元平面に対応しています。右の図は、\(N=2^n\) の場合の、状態 \(|\psi_{t}>\) の振幅を表す棒グラフです。振幅の平均値は破線で示されています。

Step 1:#

反転のオラクル \(U_{w}\) を状態に適用します \(U_{w}|\psi_{t}> = |\psi_{t'}>\)

幾何学的には、状態 \(|\psi_{t}>\)\(|s'>\) に対して反転させることに対応しています。この変換が意味することは、\(|w>\) 状態の振幅が負の値になるということで、結果として平均振幅が低くなることを意味しています。(訳注:右側のグラフで破線が下がっていることに着目してください)。

Step2:#

次に、\(|s>\) に対する追加の反転 \(U_{s}\) を適用します:\(U_{s} = 2|s> <s| - 1 \) この変換の結果、状態は \(U_{s}|\psi_{t'}>\) となり、変換 \(|\psi_{t+1}> = U_{s}U_{w}|\psi_{t}>\) が完了します。(訳注:右側のグラフでwに対応する振幅が増幅されていることに着目してください)。 2つの反転は常に回転と対応しています。\( U_{s}U_{w}\) による変換は、初期状態 \(|s>\) をアタリ\(|w>\) に近づけるような回転となります。(訳注:step 2の左側の図を参照)。\(U_{s}\) による反転の効果は、振幅の棒グラフにおいて、平均振幅での反転と解釈できます。最初の反転で平均振幅の値が低くなったので、この変換は、負の振幅をもった \(|w>\) をオリジナルの値から大雑把にいって約3倍程度増幅し、他の振幅は小さくします。その後、Step 1 に戻ってこれを繰り返します。アタリ \(w\)に近くなるまで、この処理を何回か繰り返します。

\(t\) 回繰り返した後、状態は \(|\psi_{t}> = (U_{s}U_{w})^{t}|\psi_{0}> \) に変換されます。

回転を何回適用する必要があるでしょうか? おおよそ \(\sqrt N\) 回転で十分なことが分かっています。これは、状態 \(|\psi_{t}>\) の振幅を調べることで明確になります。\(|w>\) の振幅が適用回数と共に線型的(\( \sim tN^{1/2}\))に増えていくことが見てとれます。確率ではなく振幅を扱っているので、ベクトル空間の値には平方根として入ります。そのため、この処理で増幅されるのは、ただの確率ではなく振幅です。

もし解が複数、\(M\)個ある場合、おおよそ\(\sqrt{(N/M)}\) 回転で十分なことが分かっています。

Qiskitでの実装: 2量子ビットのGroverのアルゴリズム#

では、GroverのアルゴリズムをQisKitで実装してみましょう。今回は、2量子ビットを使って、\(\ket{11}\)を見つけるGroverのアルゴリズムを実装します。

最初に環境を準備します。

#initialization
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

# importing Qiskit
from qiskit import IBMQ, BasicAer
from qiskit_ibm_provider import least_busy
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute

# import basic plot tools
from qiskit.tools.visualization import plot_histogram

上のStep1で実施したように求めたい解(今回の場合、\(|11>\))の符号を反転するオラクルを作成します。

def phase_oracle(circuit, register):
    circuit.cz(register[0], register[1])

qr = QuantumRegister(2)
oracleCircuit = QuantumCircuit(qr)
phase_oracle(oracleCircuit, qr)
oracleCircuit.draw(output="mpl")
_images/a29431e00c0828314167f1f91c65903f97fdc7ea4daa6c8bdf4b69d600252db7.png

次にStep 2で実施した平均についての反転(inversion about averate)をする回路を作成します。この回路はDiffusion回路と呼ばれることもあります。

def inversion_about_average(circuit, register):
    """Apply inversion about the average step of Grover's algorithm."""
    circuit.h(register)
    circuit.x(register)
    circuit.h(register[1])
    circuit.cx(register[0], register[1])
    circuit.h(register[1])
    circuit.x(register)
    circuit.h(register)
qAverage = QuantumCircuit(qr)
inversion_about_average(qAverage, qr)
qAverage.draw(output='mpl')
_images/48c8b5cdce4836efdb14568bcfad253f1ae0e35fce7965377ab64f62db51124a.png

回路の初めで均一な重ね合わせ状態を準備し、今まで作成した部品を合わせ、最後に測定を行います。四つの可能性のうち解は一つなので、繰り返し回数は一回であることに注意してください。

qr = QuantumRegister(2)
cr = ClassicalRegister(2)

groverCircuit = QuantumCircuit(qr,cr)
groverCircuit.h(qr)

phase_oracle(groverCircuit, qr)
inversion_about_average(groverCircuit, qr)

groverCircuit.measure(qr,cr)
groverCircuit.draw(output="mpl")
_images/c79778d6a6a71c6a3ac667b8e057a42506e4793b644a794b5a09d013569c31e3.png

シミュレーターでの実験#

上記の回路をシミュレーターで実行してみます。

backend = BasicAer.get_backend('qasm_simulator')
shots = 1024
results = execute(groverCircuit, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
_images/c212f547752e1f4e7cbb030c54400c386cd4afbd1c1ea9e6837aa38b66dc4be9.png

実行結果から分かるようにGroverのアルゴリズムは解を発見しています。

実デバイスでの実験#

以下に示すように、実デバイスで実験を行うこともできます。

Caution

以下のコードは実機を扱うためコメントアウトしています

# Load our saved IBMQ accounts and get the least busy backend device

# IBMQ.load_account()
# provider = IBMQ.get_provider(group='open')
# backend_lb = least_busy(provider.backends(simulator=False, operational=True))
# print("Least busy backend: ", backend_lb)
# run our circuit on the least busy backend. Monitor the execution of the job in the queue
# from qiskit.tools.monitor import job_monitor

# backend = backend_lb
# shots = 1024
# job_exp = execute(groverCircuit, backend=backend, shots=shots)

# job_monitor(job_exp, interval = 2)
# get the results from the computation
# results = job_exp.result()
# answer = results.get_counts(groverCircuit)
# plot_histogram(answer)

Qiskitでの実装: 補助量子ビットを用いた2量子ビットのGroverのアルゴリズム#

同じ\(|11>\)を見つける2量子ビットのGroverのアルゴリズムですが、今回は求めたい解の符号を反転する際に補助量子ビットを用いた回路を実装しましょう。補助量子ビットを用いることで、量子ビットの数が増えた際や問題が難しい際の複雑なオラクルを実装しやすくなります。

最初に環境を準備します。

#initialization
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

# importing Qiskit
from qiskit import IBMQ, BasicAer
from qiskit_ibm_provider import least_busy
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute

# import basic plot tools
from qiskit.tools.visualization import plot_histogram

求めたい解(今回の場合、\(|11>\))の符号を反転するオラクルを作成します。但し、今回は補助量子ビットを用いて\(|11>\)の場合に符号が反転するようにします。この際、正しく符号を反転させるためには事前に補助ビットを\(|1>\)にする必要があることに注意してください。

def phase_oracle(circuit, register,oracle_register):
    circuit.h(oracle_register)
    circuit.ccx(register[0], register[1],oracle_register)
    circuit.h(oracle_register)
    
qr = QuantumRegister(3)
oracleCircuit = QuantumCircuit(qr)
oracleCircuit.x(qr[2])
phase_oracle(oracleCircuit, qr,qr[2])
oracleCircuit.draw(output="mpl")
_images/b1939407867fe0d80309e3f703c06a525de9fce01bfc278fdd5e5eabcd553d11.png

次に、Diffusion回路を準備します。Diffusion回路は補助ビットには作用しないようにしてください。

def inversion_about_average(circuit, register):
    """Apply inversion about the average step of Grover's algorithm."""
    circuit.h(register)
    circuit.x(register)
    circuit.h(register[1])
    circuit.cx(register[0], register[1])
    circuit.h(register[1])
    circuit.x(register)
    circuit.h(register)
qAverage = QuantumCircuit(qr)
inversion_about_average(qAverage, qr[0:2])
qAverage.draw(output='mpl')
_images/1cffb0909c3a7305f561f800e7a8b2341221215ae0204daefe151126e59e0899.png

補助ビットなしで作成した回路と同じように、回路の最初で均一な重ね合わせ状態を作成し(この際、補助ビットにHゲートをかけないように注意してください)、変換を適用し、最後に測定をします。

qr = QuantumRegister(3)
cr = ClassicalRegister(3)

groverCircuit = QuantumCircuit(qr,cr)
groverCircuit.h(qr[0:2])
groverCircuit.x(qr[2])

phase_oracle(groverCircuit, qr,qr[2])
inversion_about_average(groverCircuit, qr[0:2])

groverCircuit.measure(qr,cr)
groverCircuit.draw(output="mpl")
_images/baf6459c2f33b8dfa8094bd67d8e3ccf3fe6418d9108e78b3bea0c0ac19da9f4.png

シミュレーターでの実験#

上記の回路をシミュレーターで実行してみます。

backend = BasicAer.get_backend('qasm_simulator')
shots = 1024
results = execute(groverCircuit, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
_images/601a2b22bb30e2bc4d3f3b99b105329b26cf6f8aa13e203bdced1a74d8671d2d.png

補助ビットなしで作成した回路と同じように\(|11>\)の状態の確率が増幅されていることが確認できます。最上位ビットの1は補助量子ビットのものなので無視して問題ありません。

今週の問題(ラーニングチャレンジ演習 Ⅱ)#

2量子ビット Groverのアルゴリズムの実装

2量子ビットを使って\(\ket{10}\)を見つけるGroverのアルゴリズムを実装してみましょう。
その際、補助量子ビットを使わないものと使うものの二パターン実装してみましょう。 ケット内のより上位の量子ビットがQiskit上の量子レジスタのインデックス数が大きいものに対応しています。

シミュレーターでの実行結果から、最も確率が高い量子ビット列とunrollerを用いてu3とCXに分解した結果を提出してください。

作成した量子回路をUnrollerを用いてu3とCXに分解する方法については、量子回路のコスト計算についてを参照してください。

解答例#

補助量子ビットなし#

#initialization
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

# importing Qiskit
from qiskit import IBMQ, BasicAer
from qiskit_ibm_provider import least_busy
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute

# import basic plot tools
from qiskit.tools.visualization import plot_histogram

\(|10>\)の状態の符号を反転するオラクルを作成します。このオラクルはCZゲートとXゲートを以下のように組み合わせることにより作成出来ます。\(|10>\)の状態の符号のみ反転していることを確認してください。

def phase_oracle(circuit, register):
    circuit.x(register[0])
    circuit.cz(register[0], register[1])
    circuit.x(register[0])

qr = QuantumRegister(2)
oracleCircuit = QuantumCircuit(qr)
phase_oracle(oracleCircuit, qr)
oracleCircuit.draw(output="mpl")
_images/a4f718e4a5bc990398e8fb709f88ead7e9a517d7aacfabffc50fc4f5633ed834.png

解説

オラクルは\(2^{n}\)通りの重ね合わせ状態の中からアタリの状態の符号だけを反転させるものなので、以下のようにXゲートとCZゲートの組み合わせで作成できます。

  1. アタリの状態で0になっているビットについてXゲートを掛けて、アタリの状態が全ビットが1である状態に置き換えられるようにします。

  2. 次にCZゲートを掛けて符号を反転させます。

  3. 先ほどXゲートを掛けたビットについてもう一度Xゲートを掛け、元のアタリの状態に戻します。

def inversion_about_average(circuit, register):
    """Apply inversion about the average step of Grover's algorithm."""
    circuit.h(register)
    circuit.x(register)
    circuit.h(register[1])
    circuit.cx(register[0], register[1])
    circuit.h(register[1])
    circuit.x(register)
    circuit.h(register)
qr = QuantumRegister(2)
cr = ClassicalRegister(2)

groverCircuit = QuantumCircuit(qr,cr)
groverCircuit.h(qr)

phase_oracle(groverCircuit, qr)
inversion_about_average(groverCircuit, qr)

groverCircuit.measure(qr,cr)
groverCircuit.draw(output="mpl")
_images/9b965954758c87c555277942702e8f9013b3936208db10926b1d5973092e600a.png
backend = BasicAer.get_backend('qasm_simulator')
shots = 1024
results = execute(groverCircuit, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
_images/3ba45c09e4e8f302c931c0e01cb9eb3be2afe63fbf991fe8346359d5ec4511cd.png

解説

本文Step1, Step2のように振幅増幅を繰り返して、\(2^{n}\)通りの重ね合わせ状態\(|s>\)\(|w>\)に近付けていきますが、今回はアイテムの総数が4個のため、振幅増幅を1回実施するだけで、理論上は(量子コンピュータのノイズを無視すると)確実にアタリを発見できます。
振幅増幅を行う前の状態は\(|s>=\frac{1}{2}(|00>+|01>+|10>+|11>)\)です。
\(|w>\)に垂直なベクトルを\(|s'>\)とすると、\(|s>=\frac{1}{2}|w>+a|s'>\)と表せます。
\(|s>\)\(|s'>\)のなす角を\(θ\)とすると、\(|s>\)の長さは1なので\(sinθ=\frac{1}{2}\)、したがって\(θ=\frac{π}{6}\)です。
\(|s>\)にオラクル\( U_{w}\)を適用すると(Step1)、\(θ=-\frac{π}{6}\)になります。
さらに\( U_{s}\)を適用すると(Step2)、\(θ=\frac{π}{2}\)となり、\(|s>\)\(|w>\)と一致しアタリを見つけられることが分かります。

2019-w2-amp.png

実際に計算してみると以下のようになり、やはり1回の振幅増幅で\(|s>\)\(|w>\)に一致します。
なお、\( U_{s}|s>=|s>\)および(アイテムの総数が4個の場合)\( U_{s}|x>=|s>-|x>\)を用いています。
\( U_{w}|s>=\frac{1}{2}\sum_{x≠w}^{3} |x> -\frac{1}{2} |w>=|s>-|w>\)
\( U_{s}U_{w}|s>=U_{s}|s>-U_{s}|w>=|s>-(|s>-|w>)=|w>\)

補助量子ビットあり#

#initialization
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

# importing Qiskit
from qiskit import IBMQ, BasicAer
from qiskit_ibm_provider import least_busy
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute

# import basic plot tools
from qiskit.tools.visualization import plot_histogram
def phase_oracle(circuit, register,oracle_register):
    circuit.h(oracle_register)
    circuit.x(register[0])
    circuit.ccx(register[0], register[1],oracle_register)
    circuit.x(register[0])
    circuit.h(oracle_register)

qr = QuantumRegister(3)
oracleCircuit = QuantumCircuit(qr)
oracleCircuit.x(qr[2])
phase_oracle(oracleCircuit, qr,qr[2])
oracleCircuit.draw(output="mpl")
_images/bf352e55c46fcfc0b997dcfb3ba64d4b6070ca54230350035d1a30a6d76e48c9.png

解説

補助量子ビットありの場合は位相反転に用いているゲートが異なりますが、これも補助量子ビットなしの場合に用いているCZゲートと同じ働きをしています。
補助量子ビットありの場合も以下のような手順でオラクルを作成できます。

  1. アタリの状態で0になっているビットについてXゲートを掛けます。

  2. 位相を反転させます。

  3. もう一度Xゲートを掛けます。

def inversion_about_average(circuit, register):
    """Apply inversion about the average step of Grover's algorithm."""
    circuit.h(register)
    circuit.x(register)
    circuit.h(register[1])
    circuit.cx(register[0], register[1])
    circuit.h(register[1])
    circuit.x(register)
    circuit.h(register)
qr = QuantumRegister(3)
cr = ClassicalRegister(3)

groverCircuit = QuantumCircuit(qr,cr)
groverCircuit.h(qr[0:2])
groverCircuit.x(qr[2])

phase_oracle(groverCircuit, qr,qr[2])
inversion_about_average(groverCircuit, qr[0:2])

groverCircuit.measure(qr,cr)
groverCircuit.draw(output="mpl")
_images/d59cad93c99618a1f9de007e983b101774a4957b0d4d7a4795c5787085ec61e5.png
backend = BasicAer.get_backend('qasm_simulator')
shots = 1024
results = execute(groverCircuit, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
_images/c0971fddda79fe736546257e07996b5647dfb287b6e41eddd279c05f4e8d1cac.png

正しく\(|10>\)の状態の確率が増幅されています(最上位の量子ビットは補助量子ビットのため、無視して問題ありません)。

Unrollerを用いた量子回路の分解#

上記解答例をUnrollerを用いてu3とCXに分解します。

from qiskit.transpiler import PassManager
from qiskit.transpiler.passes import Unroller
pass_ = Unroller(['u3', 'cx'])
pm = PassManager(pass_)
new_circuit = pm.run(groverCircuit) 
new_circuit.draw(output='mpl')
_images/c396573ceb284ba6c4d41507ec70dea59ea7ca41a840d66c3bfed94d1cc91ce1.png

解答例の量子回路のコストは、
補助量子ビットなしの場合:\(16+2\times10=36\)
補助量子ビットありの場合:\(26+7\times10=96\)
となります。