グローバーのアルゴリズム基礎(3)：3量子ビット以上の場合

前回まででグローバーのアルゴリズムの基本要素である「オラクル回路」「位相反転増幅回路」の理論と実装の対応を解説しました。
前回までは2量子ビットの場合でしたが、一般に3量子ビットの場合の実装方法を解説します。

シリーズ予定：

オラクル回路
位相反転増幅回路
3量子ビット以上の場合 <<今回
まとめと簡単な応用例

対象とする読者

量子プログラミング・量子計算に興味を持って勉強し始めたエンジニア・高校生・大学生
Groverのアルゴリズムについて改めて知りたい方など

前提とする知識

簡単な論理演算(AND, OR, NOT)が分かること
基本的な量子ゲート(H, X, CX, CCX, Z, CZ)が分かること
線形代数(行列・ベクトルの積、直積、あとできれば線形変換)が分かること
ベクトルのブラケット記法( $|0>, |1>$ など、もしくはDirac記法ともいう)が分かること

バージョン情報

Python 3.7.3
Qiskit 0.23.1

Groverのアルゴリズムのおさらい
3量子ビット以上の場合に生じる課題
3量子ビット以上でのGroverのアルゴリズムの実装方法
まとめ

Groverのアルゴリズムのおさらい

オラクル回路(量子オラクル)

オラクル回路は見つけ出したい量子状態をマークして、その状態の位相を反転する回路でした。一般に量子オラクル $U_{\omega}$ はマーキングしたい状態を $| x_i^* \gt$ とし、その個数を $M$ とすると下記のように表します。

$\displaystyle{ U_{\omega} = I - 2 \sum_{i = 1}^{M} | x_i^* \gt \lt x_i^*| }$

2量子ビットの場合で $\sum_{i = 1}^{M} | x_{i}^{*} \gt \lt x_{i}^{*} | = | 11 \gt \lt 11 |$ の場合、この量子オラクルを通過後の状態 $| s' \gt$ は下記になります。

$\begin{aligned} U_{\omega} | s > &= U_{\omega} (| 00 \gt + |01 \gt + | 10 \gt + | 11 \gt) / 2 \\ &= (| 00 \gt + |01 \gt + | 10 \gt - | 11 \gt) / 2 \\ &\equiv | s' \gt \end{aligned}$

位相反転増幅回路

マーキングした状態の確率振幅の大きさを増幅する位相反転増幅は下記です。

$\displaystyle{ U_{s} = 2 | s \gt \lt s | - I }$

2量子ビットの場合で $| 11 \gt$ がマーキングした状態とすると、位相反転増幅後は

$\begin{aligned} U_s | s' \gt &= H^{\otimes 2} (2 | 00 >< 00 | - I ) H^{\otimes 2} \frac{1}{2}(| 00 \gt + |01 \gt + | 10 \gt - | 11 \gt) \\ &= H^{\otimes 2} (2 | 00 >< 00 | - I ) \frac{1}{2}(| 00 \gt + |10 \gt + | 01 \gt - | 11 \gt) \\ &= H^{\otimes 2} \frac{1}{2} (| 00 \gt - | 01 \gt - | 10 \gt + | 11 \gt ) \\ &= |11> \end{aligned}$

ということでマーキングした状態を取り出すことができました。

2量子ビットの場合の量子回路

上記のオラクル回路・位相反転増幅回路を組み合わせたGrover回路の一例は下図です。

f:id:sakumadaisuke32:20210214152132p:plain — 図1. 2量子ビットの場合のGrover回路

3量子ビット以上の場合に生じる課題

結論から説明すると、オラクル回路が現状のCZゲートによる実装方法では3量子ビット以上で困ります。

3量子ビットでは初期状態(オラクル回路への入力 $| s \gt$ )が次の式のようになります。

$\displaystyle{ | 000 \gt + | 001 \gt + | 010 \gt + | 011 \gt + | 100 \gt + | 101 \gt + | 110 \gt + | 111 \gt }$

例えば $| 111 \gt$ がマーキングしたい状態だとして、2量子ビットゲートであるCZゲートだけではこの位相を反転させられません。

なお、位相反転増幅回路は元のまま、MultiControlledToffoliでよいです。

3量子ビット以上でのGroverのアルゴリズムの実装方法

3量子ビット以上の場合の量子回路

結論から、下記の量子回路を実装すればよいです。
(もちろん、マーキングを達成できればこれ以外でもいいと思います。ただ冗長な回路になると思います。)

qr = QuantumRegister(3, 'data')
oracle = QuantumRegister(1, 'oracle')
cr = ClassicalRegister(3, 'output')
grover = QuantumCircuit(qr, oracle, cr)

# 重ね合わせ状態の生成
grover.h(qr[:3])
grover.barrier()

# オラクル
grover.x(oracle)
grover.h(oracle)
grover.mct(qr, oracle)
grover.h(oracle)
grover.x(oracle)
grover.barrier()

# 位相反転増幅
grover.h(qr)
grover.x(qr)
grover.barrier(qr)
grover.h(qr[-1])
grover.mct(qr[:-1], qr[-1])
grover.h(qr[-1])
grover.barrier(qr)
grover.x(qr)
grover.h(qr)

# 測定
grover.measure(qr, cr)

grover.draw(output='mpl')

f:id:sakumadaisuke32:20210307185745p:plain — 図2. 3量子ビット(4量子ビット目は補助ビット)のGrover回路

※4つ目の量子ビット('oracle'と名付けたはオラクル回路のための補助ビットで、そのあとの演算では用いられません。

以降、N量子ビットのGroverのアルゴリズムを実装する場合はMCTゲートを下に伸ばしていく形で、最低でも(N + 1)量子ビットあれば実装可能です。

オラクル回路が正しく動作していることの確認

図2の量子回路が正しくGroverのアルゴリズムの実装となっていることを確かめます。

まず、オラクル回路 $U_{\omega}$ の入力状態は下記です。

$\displaystyle{ | s \gt = \frac{1}{2\sqrt2}(| 0000 \gt + | 0010 \gt + | 0100 \gt + | 0110 \gt + | 1000 \gt + | 1010 \gt + | 1100 \gt + | 1110 \gt ) }$

第4量子ビット(一番右)は図2の補助ビット('oracle')です。

そして、マーキングしたい状態が $| 111 \gt$ (対応する入力は $| 1110 \gt$ )の場合、オラクル回路の出力は

$\begin{aligned} | s' \gt &= U_{\omega}| s > \\ &= (I - 2 | 111 \gt \lt 111|)|s> \\ &= \frac{1}{2\sqrt2} (| 0000 \gt + | 0010 \gt + | 0100 \gt + | 0110 \gt + | 1000 \gt + | 1010 \gt + | 1100 \gt - | 1110 \gt ) \end{aligned}$

でマーキング状態だけ位相が反転するようになっています。

シミュレータで動作を確認します。

statevector_simulatorで動作確認

ソース

qr = QuantumRegister(3, 'data')
oracle = QuantumRegister(1, 'oracle')
cr = ClassicalRegister(3, 'output')
grover = QuantumCircuit(qr, oracle, cr)

# 重ね合わせ状態の生成
grover.h(qr[:3])
grover.barrier()

# オラクル
grover.x(oracle)
grover.h(oracle)
grover.mct(qr, oracle)
grover.h(oracle)
grover.x(oracle)
grover.barrier()

backend = Aer.get_backend('statevector_simulator')
job = execute(grover, backend=backend)
result = job.result()

states = result.get_statevector()
state_key = 0
for state in states: 
    print(f'{int(bin(state_key)[2:]):04}: {state}')
    state_key += 1

grover.draw(output='mpl')

f:id:sakumadaisuke32:20210307192821p:plain — 図3. 3量子ビットのオラクル回路

この状態ベクトルの確率振幅のシミュレーション結果(状態: 確率振幅)は下記です。
(左から順に第4量子ビット、・・・、第1量子ビットで並んでいます。)

0000: (0.35355339059327384-1.2989340843532405e-16j)
0001: (0.35355339059327384-1.2989340843532405e-16j)
0010: (0.35355339059327384-1.2989340843532405e-16j)
0011: (0.35355339059327384-1.2989340843532405e-16j)
0100: (0.35355339059327384-1.2989340843532405e-16j)
0101: (0.35355339059327384-1.2989340843532405e-16j)
0110: (0.35355339059327384-1.2989340843532405e-16j)
0111: (-0.35355339059327384+1.7319121124709873e-16j)
1000: (-1.1496735851465458e-17+4.329780281177469e-17j)
1001: (-1.1496735851465458e-17+4.329780281177469e-17j)
1010: (-1.1496735851465458e-17+4.329780281177469e-17j)
1011: (-1.1496735851465458e-17+4.329780281177469e-17j)
1100: (-1.1496735851465458e-17+4.329780281177469e-17j)
1101: (-1.1496735851465458e-17+4.329780281177469e-17j)
1110: (-1.1496735851465458e-17+4.329780281177469e-17j)
1111: (1.1496735851465458e-17-5.166667716518644e-33j)

取りうる状態がすべてシミュレートされるので1000以降は不要なのですが、大きさが限りになくゼロに近い振幅で表示されます。

0111以前を確認するとマーキングしたい0111だけ大きさが負になっており、目的の回路ができています。

位相反転増幅回路が正しく動作していることの確認

位相反転増幅回路の入力状態 $| s' \gt$ は下記です。

$\begin{aligned} | s' \gt &= \frac{1}{2\sqrt2} (| 0000 \gt + | 0010 \gt + | 0100 \gt + | 0110 \gt + | 1000 \gt + | 1010 \gt + | 1100 \gt - | 1110 \gt ) \\ &= \frac{1}{2\sqrt2} (| 000 \gt + | 001 \gt + | 010 \gt + | 011 \gt + | 100 \gt + | 101 \gt + | 110 \gt - | 111 \gt ) \otimes | 0 \gt_4 \end{aligned}$

第4量子ビットには何も演算されないので位相反転増幅した結果は下記となります。

$\begin{aligned} U_s | s' \gt &= \{ H^{\otimes 3} (2 | 000 \gt \lt 000 | - I ) H^{\otimes 3} \frac{1}{2\sqrt2}(| 000 \gt + | 001 \gt + | 010 \gt +\\ & | 011 \gt + | 100 \gt + | 101 \gt + | 110 \gt - | 111 \gt ) \} \otimes | 0 \gt_4 \\ &= \frac{1}{4\sqrt{2}} (| 000 \gt + | 001 \gt + | 010 \gt + | 011 \gt + \\ &| 100 \gt + | 101 \gt + | 110 \gt + 5 | 111 \gt) \otimes |0>_{4} \end{aligned}$

というように、マーキングした $| 111 \gt$ の確率振幅が大きくなっています。
また、その大きさが規格化されていることも分かります。

$\displaystyle{ \frac{1}{32} \times 7 + \frac{25}{32} = 1 }$

statevector_simulatorで動作確認

ソース

qr = QuantumRegister(3, 'data')
oracle = QuantumRegister(1, 'oracle')
cr = ClassicalRegister(3, 'output')
grover = QuantumCircuit(qr, oracle, cr)

# 重ね合わせ状態の生成
grover.h(qr[:3])
grover.barrier()

# オラクル
grover.x(oracle)
grover.h(oracle)
grover.mct(qr, oracle)
grover.h(oracle)
grover.x(oracle)
grover.barrier()

# 位相反転増幅
grover.h(qr)
grover.x(qr)
grover.barrier(qr)
grover.h(qr[-1])
grover.mct(qr[:-1], qr[-1])
grover.h(qr[-1])
grover.barrier(qr)
grover.x(qr)
grover.h(qr)

backend = Aer.get_backend('statevector_simulator')
job = execute(grover, backend=backend)
result = job.result()

states = result.get_statevector()
state_key = 0
for state in states: 
    print(f'{int(bin(state_key)[2:]):04}: {state}')
    state_key += 1

grover.draw(output='mpl')

f:id:sakumadaisuke32:20210314105947p:plain — 図4. 3量子ビットの場合のGrover回路

この状態ベクトルの確率振幅のシミュレーション結果(状態: 確率振幅)は下記です。 (左から順に第4量子ビット、・・・、第1量子ビットで並んでいます。)

0000: (-0.1767766952966371+1.4071785913826776e-16j)
0001: (-0.17677669529663703+9.74200563264931e-17j)
0010: (-0.176776695296637+9.742005632649305e-17j)
0011: (-0.17677669529663703+5.4122253514718377e-17j)
0100: (-0.17677669529663703+1.190689577323804e-16j)
0101: (-0.17677669529663698+3.2473352108831075e-17j)
0110: (-0.17677669529663695+3.247335210883104e-17j)
0111: (-0.8838834764831851+8.551316055325503e-16j)
1000: (5.748367925732717e-18-3.247335210883103e-17j)
1001: (5.748367925732721e-18-3.247335210883103e-17j)
1010: (5.748367925732723e-18-3.247335210883103e-17j)
1011: (5.748367925732726e-18-3.247335210883103e-17j)
1100: (5.74836792573272e-18-3.247335210883103e-17j)
1101: (5.748367925732726e-18-3.247335210883103e-17j)
1110: (5.7483679257327305e-18-3.247335210883103e-17j)
1111: (2.874183962866363e-17-7.577115492060576e-17j)

0111の状態が一番確率振幅の大きさが大きくなっています。