科学しよう

量子計算のプログラミングの解説をメインに、データサイエンス・機械学習について勉強したことをご紹介します

MENU

PennyLaneを用いた量子回路最適化

僕が学生時代に量子情報処理の研究に関わっていたのはこちらで自己紹介した通りで、 研究者ではありませんが今も量子情報処理・量子計算には興味があります。
そして今は機械学習エンジニア・データサイエンティストとしてソフトウェア開発に携わっています。
そんな自分が今量子コンピュータの界隈でできたらいいな・機会があれば取り組みたい・自分の経歴を生かすことができると思っていることが、量子機械学習です。

量子機械学習についてはQiitaでKuma氏がまとめてくださっています
主に以下の目的で研究されています。

ただ量子計算も機械学習も両方わからないといけないのでなかなかハードです。
みんなが量子コンピュータ使える頃には使い方をちゃんと知っている人しか使いこなせないという参入障壁が高い状態になっていると困ります。
だから今のうちにやりたい、やりましょう!

ここからは量子機械学習を実装したフレームワークの一つであるPennyLaneのチュートリアルをベースに、要素技術を深掘りして解説していきます。
今回のシリーズは量子回路のパラメータ最適化(回転ゲートの回転角の調整)です。
こちらのチュートリアルです。

対象とする読者

  • 量子計算に興味をがあるエンジニア・高校生・大学生の方
  • 量子機械学習に興味がある機械学習エンジニア・データサイエンティスト など

前提とする知識

バージョン情報

  • Python 3.9.5
  • PennyLane 0.15.1

目次

PennyLaneとは

PannyLaneはカナダのXanado社が開発している量子機械学習を実装した(量子機械学習をするのに必要な・あったら便利な関数やクラスを用意してくれている)ライブラリの一つで、機械学習フレームワークであるPyTorchをベースにしています。

pennylane.ai

PyTorchで作られているので計算結果をPyTorchに連携できます。
ただし、現状では自分のPC内にシミュレータを構築して計算しているので、実機ではありません。

その他に量子機械学習を実装しているフレームワークは2021年6月現在把握してある限り以下です。

あまりないですね、、
中でもこのPennyLaneはチュートリアルが豊富で、ノートブックも用意されているのでとりあえず実行することができます!

コアとなる概念: 量子ノード(Quantum Node, QNode)

PennyLaneを使う上で重要な概念が量子ノード(Quantum Node, QNode)です。
下図のように量子ビットや量子ゲートのパラメータ(入力)、測定結果(出力)、そしてその間の何らかの量子回路からなる計算単位です。
このQNodeを一つの計算グラフとして計算結果を後続の処理に渡したり、組み合わせて量子ニューラルネットワークを構築するのでしょう。
今回のシリーズではQNodeは1つしか使わないので、入力から量子ゲートを経て出力までをセットにした関数と思っていただければ良いです。

f:id:sakumadaisuke32:20210608003041p:plain
図1. 量子ノード概念図

量子回路の計算結果を目的関数として最小化問題を解く

古典の機械学習では解きたい課題を表す目的関数を定義し、様々なパラメータでその目的関数を計算し、目的関数の値が最小もしくは最大となるパラメータの組み合わせを見つけます。
量子機械学習でも同様に、量子回路での計算結果を使って目的関数を定義し、量子回路のパラメータを色々変えて目的関数を計算し、目的関数の値が最小もしくは最大となるパラメータの組み合わせを見つけます。 通常、TensorFlowやPyTorchなどの機械学習フレームワークでは最小とするように計算するようになっています。

今回の解く問題と量子回路

今回解こうとしている問題は「1つの量子ビットの状態を | 0 \gtから | 1 \gtに遷移させるための \mathrm{R_x},  \mathrm{R_Y}ゲートの角度はいくらか」(図2)です。
これは明らかにXゲートもしくはYゲートをかければ良いのですが、あえてこれを探索的に見つけようという練習課題です。

f:id:sakumadaisuke32:20210612084233p:plain
図2. 問題設定を示すBloch球

量子回路は下記です。

from pennylane import numpy as np

dev1 = qml.device("default.qubit", wires=1)

@qml.qnode(dev1)
def circuit(params):
    qml.RX(params[0], wires=0)
    qml.RY(params[1], wires=0)
    return qml.expval(qml.PauliZ(0))

drawer = qml.draw(circuit)
print(drawer([0.54, 0.12]))

circuit関数で量子回路を定義し、QNodeとして扱うためにqnodeデコレータを指定しています。
この回路のprintした結果は下記です。

 0: ──RX(0.54)──RY(0.12)──┤ ⟨Z⟩ 

可視化のために \mathrm{R_x}ゲートと \mathrm{R_Y}ゲートの角度を指定していますが、これをこの後で探索します。
なお、角度の単位はradianです。

今回の誤差関数

PennyLaneでは測定にPauli演算子を使います。
今回はZ軸への射影測定を想定してPauliZ測定します。

pennylane.readthedocs.io

演算子 \hat{\sigma}_zの期待値は状態 | 0 \gtについては1,  | 1 \gtについては-1であることを利用して、下記となります。
上で定義した量子回路そのままです。

def cost(x):
    return circuit(x)

誤差関数を数式で表すと、演算子 \hat{\sigma}_zを任意の状態 | \psi \gtで期待値を求めるということなので、細かい計算は割愛しますが下記になります。


\begin{aligned}
\lt \hat{\sigma}_z \gt &= \lt \psi | \hat{\sigma}_z | \psi \gt \\
&= \lt 0 | \hat{R}_Y^\dagger (\phi) \hat{R}_X^\dagger (\theta) \hat{\sigma}_z R_Y (\theta) R_X (\phi) | 0 \gt \\
&= \cos \theta \cos \phi
\end{aligned}

 ( \theta , \phi) = (0, 0)の時に最大で \lt \hat{\sigma}_z \gt = 1,  ( \theta , \phi) = (0, \pi)もしくは (\theta , \phi) = (\pi, 0)の時に最小で  \lt \hat{\sigma}_z \gt = -1となります。

勾配降下法で最適化

勾配降下法で最適化します。
コードとしては以下です。

# initialise the optimizer
opt = qml.GradientDescentOptimizer(stepsize=0.4)

# set the number of steps
steps = 100
# set the initial parameter values
init_params = np.array([0.011, 0.012])
params = init_params

for i in range(steps):
    # update the circuit parameters
    params = opt.step(cost, params)

    if (i + 1) % 5 == 0:
        print("Cost after step {:5d}: {: .7f}".format(i + 1, cost(params)))

print("Optimized rotation angles: {}".format(params))

やることは古典の機械学習と同じです。
誤差関数の測定値と勾配に基づいて、パラメータ更新をします。

誤差関数の勾配は下記です。


\begin{aligned}
\nabla \lt \hat{\sigma}_z \gt &= \left( \cos \phi \frac{\partial}{\partial \theta} \cos \theta , \cos \theta \frac{\partial}{\partial \phi} \cos \phi \right) \\
&= \left( - \sin \theta \cos \phi  , - \cos \theta \sin \phi \right) \\
\end{aligned}

計算したらわかりますが ( \theta , \phi) = (0, 0)の時に勾配もゼロとなりパラメータ更新されなくなってしまうので、初期値を少しずらしてあるのがミソです。

パラメータ設定と誤差測定を繰り返して数値的に最適なパラメータを探索する

上のコードを実行すると下記のように誤差関数の更新が行われて、最初に期待した通りにY軸周りに \pi回転する場合が解として求まりました。

Cost after step     5:  0.9961778
Cost after step    10:  0.8974944
Cost after step    15:  0.1440490
Cost after step    20: -0.1536720
Cost after step    25: -0.9152496
Cost after step    30: -0.9994046
Cost after step    35: -0.9999964
Cost after step    40: -1.0000000

(中略)

Cost after step   100: -1.0000000
Optimized rotation angles: [7.15266381e-18 3.14159265e+00]

まとめ

問題設定を量子回路で表現するまでできれば、あとは古典の機械学習と同じ流れでできることを確かめました。

  • PennyLaneでは量子回路をQNodeで定義する
  • 問題設定・誤差関数を量子回路で表現する
  • パラメータ最適化は古典機械学習と同じく誤差関数の勾配から勾配降下法で行える


【送料無料】 量子コンピュータによる機械学習 / Maria Schuld 【本】

赤池情報量基準(AIC)を多項式の回帰で実装

こんにちは、さくまっちょです。

今回は本業のデータサイエンス関連で気になっていたこととして赤池情報量基準(AIC)を実際に実験用にデータを作って計算し その意味や威力を確認してみました。

機械学習・データサイエンスやっているとモデルの汎化性能を上げるために複雑なモデルを使わない方が良いという説明が PRMLやはじパタなど各種機械学習系のテキストには必ずと言っていいほど載っています。
数式もちゃんと載っているのですが、いざ自分で計算してみようと思っても実装をしている記事などがなかなか出てこず、仕方がないので数式を調べながら実装してみました。

同じ疑問をお持ちの方のお役に立てば幸いです。

もしよければコメントいただきたいのですが普段みなさんってAICって計算しないのでしょうか?
もしくはいくつかのライブラリや計算クラスではAICやそれに準ずるモデル選択のメソッドが用意されているから自分では実装しないのでしょうか?

すべての実装についてはこちらの僕のノートブックをご覧ください。

対象とする読者

  • 機械学習エンジニア
  • データサイエンティスト
  • 機械学習・データサイエンスに興味のある学生・エンジニア

前提とする知識

  • 回帰分析がわかること
  • 基本的な数学(sin, logや線形代数)

バージョン情報

目次

赤池情報量基準(AIC)とは

Wikipediaによると

  • 統計数理研究所所長の赤池弘次が1971年に考案し1973年に発表した
  • 統計モデルの良さを評価するための指標
  • 「モデルの複雑さと、データとの適合度とのバランスを取る」ために使用される
  • AIC最小のモデルを選択すれば、多くの場合、良いモデルが選択できる

です。 また、詳細な意味や議論については赤池氏の記事をぜひ読んでみてください。
よく比較されるBICとの違いについてもご本人が議論されています。

詳細な議論はともかくとして、AICを使うことによって様々に作ったモデルの中から「これが一番データに合ってそうだ」というモデルを決めるのに役立ちます。

多項式回帰での次数と汎化性能

PRMLやはじパタで説明されているように、僕もまずはデータを生成し多項式で回帰して汎化性能を確かめてみました。

作成したデータがこちらです。

f:id:sakumadaisuke32:20210523171654p:plain
図1. 実験データ. 実線が理想的な正弦波. プロットが実際に訓練データとして使うために正弦波にノイズを与えたデータ.

これに0次から20次までの多項式で回帰した結果がこちらです。

f:id:sakumadaisuke32:20210523172024p:plain
図2. 多項式で回帰した結果. 黒実線が理想的な正弦波. 緑実線が回帰した結果. 赤点が訓練データ. 青点が検証データ.

6次くらいまでは正弦波っぽく描かれていますが、7次以降は赤点(訓練データ)との距離は縮んでそうですが歪みがあるのと青点(検証データ)との距離は離れてそうです。

実際に訓練データと検証データの回帰線とのRMSEを比較したのがこちらです。

f:id:sakumadaisuke32:20210523172436p:plain
図3. 多項式の次数とRMSE

3次までは訓練・検証ともRMSEは下がっていますね。
それ以降も横ばいな感じはありますが、訓練データは10次まで、検証データも9次まではRMSEが改善してそうです。
このグラフから得られる結果とグラフの形状から、なんとなく3次~6次の間でモデルを決定したいです。

ここで「なんとなく」ではなく定量的に次数を決めるための指標がAICです。

データの分布が正規分布を仮定した場合のAIC

AICの定義

Wikipediaによると

\displaystyle{
\mathrm{AIC} = -2 \ln {L} + 2 k
}

 Lは尤度、 kはパラメータの数(今回は多項式の次数 + 1(定数項を含むため)) です。

尤度を計算するのに確率分布を導入しないといけなくて、このあたりの実装が厄介だったのですがこちらの41ページで正規分布の場合は簡単に整理できるという資料がありましたので、この通りに計算してみたのが次です。

データの分布が正規分布を仮定した場合のAIC

\displaystyle{
\mathrm{AIC} = n \ln{2 \pi} + n \ln{\mathrm{MSE}} + n + 2(k + 2)
}

 nはデータの個数、 \mathrm{MSE}は平均自乗誤差です。 実装はこちらになります。

from sklearn.metrics import mean_squared_error as MSE

def AIC(ys_actual, ys_pred, num_feature):
    num_data = len(ys_actual)
    mse = MSE(ys_actual, ys_pred)
    return num_data * np.log(2 * np.pi) + num_data * np.log(mse) + num_data + 2 * (num_feature + 2)

いざ書いてみたらシンプルでした。

AICを使った多項式回帰の次数の決定

上記のAICの式を用いて各次数について計算してみて最適次数を決定しました。

f:id:sakumadaisuke32:20210523173758p:plain
図4. 多項式の次数とAICの関係. 青線が訓練データのAIC. オレンジ線が検証データのAIC

検証データのAICが最小になる3次にするのがよさそうと分かりました。

以上でAICの実装と使い方、そして実際に次数を決められることを示しました。

まとめ

  • AICはモデルの複雑さとデータの適合度を評価する指標
  • AICを使うと複数作成したモデルから最適なものを決めらる


パターン認識と機械学習 上 ベイズ理論に/C.M.ビショップ元田浩【3000円以上送料無料】


はじめてのパターン認識 [ 平井 有三 ]


Pythonで動かして学ぶ!あたらしい機械学習の教科書 第2版 (AI & TECHNOLOGY) [ 伊藤 真 ]

qRAMの作り方

qRAM(Quantum Random Access Memory)を量子回路上で実装する方法を解説します。 こちらこちらの論文をベースにエッセンスを解説し、これがあることで解ける問題の幅が広がることを紹介します。

対象とする読者

  • 量子プログラミングに興味を持ち始めたばかりのエンジニア
  • 量子計算に興味を持ち始めたばかりの高校生・大学生 など

前提とする知識

  • 簡単な論理演算(AND, OR, NOT)が分かること
  • 基本的な量子ゲート(X, CX, CCX)が分かること

バージョン情報

目次

qRAM(Quantum Random Access Memory)とは

定義

\displaystyle{
\sum_j{\psi _j | j \gt_a} \overset{\mathrm{qRAM}}{\longrightarrow} \sum_j{\psi _j | j \gt_a} | D_j \gt_d
}

右辺に注目してもらえればいいのですが、入力の状態 | j \gt_aに対してデータ | D \gt_dエンタングルして1対1で紐づいているという状態を指します。 つまり、通常のプログラミングでの辞書型のデータ構造に相当しており、あるキー(状態)を用いてデータを取り出して以後の処理で用いることができます。 注意としてはqRAMという名称ですがあくまでデータ構造なので古典コンピュータのRAM(メインメモリ)のようにそういうハードウェアがあるのではなく、量子回路の組み合わせで構成します。

もちろん量子なので入力が重ね合わせ状態なら取り出せるデータも重ね合わせで取り出せます。

今回の記事ではこのqRAMのデータ構造の作り方・量子回路の組み方を解説します。

RAM(Random Access Memory)って何?

WikipediaによるRAMの説明を整理すると下記です。

  • コンピュータで使用するメモリ(主記憶装置)の一分類
  • 格納されたデータに任意の順序でアクセスできる(ランダムアクセス)メモリ
  • 電源が落ちれば記憶内容が消えてしまう短期メモリの意で使われていることが専ら

つまり、普段プログラミングするときに見聞きする「メモリ」です。 プログラムそのものやプログラムが用いるデータなどを実行している間だけメモリ上に展開して、必要な時にあるアドレスに格納されている値を読み取って処理を進めます。

古典RAMと量子RAM(qRAM)の違い

ありません、ほぼ同じです。

強いてあげるなら、前述したようにqRAMはqRAMというハードウェアを指すのではなく、メモリ動作をする状態を指します。 ただ、古典RAMもトランジスタをメモリとして扱えるようにパッケージしたものなので、使うビットが古典ビットor量子ビット、回路が古典回路or量子回路かという違いだけかもしれません。

もちろん重ね合わせ状態を使えるのはqRAMだけです。

qRAMがあると何ができる?

古典コンピュータと同じく、そのプログラムにとって必要なデータをプロセッサのすぐ近くに置いておくことができ、アクセスしやすくなります。 つまり、パフォーマンスがよくなります。

量子コンピュータもプログラムを実行しないときはデータをハードディスクのような古典の補助記憶装置にデータを保存していると思います。 そのデータを使おうとするたびにハードディスクにアクセスしていたら待ち時間のちりが積もって長くなります。 これは古典コンピュータでも同じですが、そのような処理の待ち時間を改善します。 また、計算結果で後々使いたいデータもいちいちハードディスクなどに一時保存して使いたいときに参照するのも時間がかかってしまいます。 すぐ近くにあるのがいいです。

また、このデータ構造があるおかげで実行できるアルゴリズムと解ける問題・アプリケーションがあります。 それが以下です。

qRAMが必要なアプリケーション

探索・検索などの計算量が大きい問題で威力が発揮されるようです。 これらはいずれもID-データというデータ構造をしており、答えとして知りたいものはIDですが計算するのはデータとなっています。 qRAMのデータ構造していると都合がよさそうです。

QRAMを実装するときの考え方

定義で紹介したこちらのエンタングルメントが形成されるように、CXゲートMCTゲートを組み合わせればよいです。

QRAMの実装例

例えば、次のような辞書で古典的なデータがあると仮定し、それを量子回路にエンコードします。

ID IDを示す状態 値を示す状態
0 |00> 2 |0010>
1 |01> 4 |0100>
2 |10> 6 |0110>
3 |11> 8 |1000>

この表の組み合わせを満たすようにエンタングルメントを形成すればよいです。

inputs = QuantumRegister(2, 'input')
data = QuantumRegister(4, 'data')
outputs = ClassicalRegister(6, 'output')
circuit = QuantumCircuit(inputs, data, outputs)

ram = {
    '00': '0010', 
    '01': '0100', 
    '10': '0110', 
    '11': '1000', 
}

def qram(circuit, ram):
    for key, val in ram.items(): 
        for idx_input in range(len(inputs)):
            if key[idx_input] == '0':
                circuit.x(inputs[idx_input])
            
        for idx_data in range(len(data))[::-1]:
            if val[idx_data] == '1':
                circuit.mct(inputs, data[len(data) - (idx_data + 1)])

        for idx_input in range(len(inputs)):
            if key[idx_input] == '0':
                circuit.x(inputs[idx_input])
        circuit.barrier()

circuit.h(inputs)
circuit.barrier()

qram(circuit, ram)

circuit.measure(inputs, outputs[:len(inputs)])
circuit.measure(data, outputs[len(inputs):len(inputs) + len(data)])

backend = Aer.get_backend('qasm_simulator')
job = execute(circuit, backend=backend)
result = job.result()

count = result.get_counts()

circuit.draw(output='mpl')

f:id:sakumadaisuke32:20210505180437p:plain
図1. qRAM動作確認回路. barrier()で仕切られた枠毎に入力状態はそれぞれ左から|00>, |01>, |10>, |11>を示します.(0, 1の数式上の表記と回路上の位置は順序逆であることに注意してください.)

ちゃんとエンタングルしているか測定して確認します。

f:id:sakumadaisuke32:20210505180956p:plain
図2. 図1の回路の測定結果

状態を示す0/1の並びの左4桁がデータで右2桁が入力状態です。 上の表の通りにエンタングルしていますね。

QuantumChallenge2020の例

その他のqRAMの応用例としてQuantumChallenge2020のコードを紹介します。

まとめ

  • qRAMはIDとデータという辞書型のデータ構造
  • 辞書を満たすようにエンタングルメントを形成すれば作れる
  • 探索・検索系のアルゴリズムで重要

参考


量子コンピューティング 基本アルゴリズムから量子機械学習まで [ 情報処理学会出版委員会 ]


コンピュータアーキテクチャ技術入門 ーー高速化の追求×消費電力の壁【電子書籍】[ Hisa Ando ]

グローバーのアルゴリズム基礎(4):まとめと簡単な応用例

すみません、忙しくて更新が遅くなりました。
このシリーズ最後はGroverのアルゴリズムの簡単な応用を紹介します。
(自分で考えたトイプロブレムです。)

シリーズ予定:

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

前回までで入力の量子ビットについて、特定の状態について振幅を増幅できるようになりました。
ただこれだけだと入力Nビットの何番目のビットが立った状態を増幅すればいいというだけなので、実行する前から答えがわかっている計算しかできません。

より実用に近づけるには、「入力Nビットについてある条件を満たすビットを知りたいのだけど、予めわからないのでGroverのアルゴリズムで探索・特定したい」という課題にアプローチできるようになる必要があります。
その一例はQuantumChallenge2020の出題の通りですが、これはQRAMとかやや複雑な数理があるので、今回は単純な問題設定にして解説します。

今回の問題設定は「a + b = 0もしくは2の解となるa, bをGroverのアルゴリズムで見つける」です。
(= 1となる解は増幅されないのでやりません。
理由はわからないのですが僕の理解が追いついていないからだと思うので、わかったらまたポストします。
もしくはご存じの方がいたらコメントください、、。)

対象とする読者

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

前提とする知識

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

バージョン情報

目次

問題設定

問題は「a + b = 0もしくは2の解となるa, bをGroverのアルゴリズムで見つける」です。
この場合の答えはa = 0, b = 0a = 1, b = 1で人間にとってはすぐわかりますが、これをGroverのアルゴリズムで解きます。

なお、a + b = 1の解はa = 0, b = 1もしくはa = 1, b = 0でこちらも人間にとってはすぐわかりますが、Gorverのアルゴリズムで増幅されないので今回は対象外とします。
理由はわかってないのですが、複数の解が含まれる場合は増幅されないのかもしれません。
数理的に理解ができたら改めて解説しようと思います。

ラクル回路の構築

ラクル回路の演算子の確認

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

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

例えばa + b = 0を満たす入力はa = 0, b = 0なので状態としては | 00 \gtを、 またa + b = 2を満たす入力はa = 1, b = 1なので状態としては | 11 \gtをマーキング(位相を反転)するような回路を構成できれば良いということです。

条件を満たす量子状態の振幅を反転

a + b = 0を満たす解

これを満たす入力を検出する回路は下記です。
最終的にMCTゲートで検出しています。

共通関数

def XOR(qc, control1, control2, target):
    qc.cx(control2, target)
    qc.cx(control1, target)
    return qc

def AND(qc, control1, control2, target):
    qc.ccx(control1, control2, target)
    return qc
    
def adder(qc, input1, input2, carry, summ):
    """半加算器
    """
    qc = XOR(qc, input1, input2, summ)
    qc = AND(qc, input1, input2, carry)
    return qc
inputs = QuantumRegister(2, 'input')
calcs = QuantumRegister(2, 'calculation')
oracle = QuantumRegister(1, 'oracle')
c = ClassicalRegister(3)
grover = QuantumCircuit(inputs, calcs, oracle, c)

grover.h(inputs)
grover.barrier()

grover = adder(grover, inputs[0], inputs[1], calcs[1], calcs[0]) # calcs[1]がcarry calcs[0]がsum
grover.x(calcs)
grover.mct(calcs, oracle)
grover.barrier()

grover.measure(inputs, c[:2])
grover.measure(oracle, c[2])

f:id:sakumadaisuke32:20210411170136p:plain
図1. a + b = 0 の解を検出する場合

この測定結果は下記です。

f:id:sakumadaisuke32:20210411170308p:plain
図2. 図1の回路の測定結果

この時点ではすべて等確率ですが、MCToracleビットが反転しているのは入力が00の場合(図2の一番右)のみです。

a + b = 2を満たす解

同様にこれを満たす回路は下記です。

inputs = QuantumRegister(2, 'input')
calcs = QuantumRegister(2, 'calculation')
oracle = QuantumRegister(1, 'oracle')
c = ClassicalRegister(3)
grover = QuantumCircuit(inputs, calcs, oracle, c)

grover.h(inputs)
grover.barrier()

grover = adder(grover, inputs[0], inputs[1], calcs[1], calcs[0]) # calcs[1]がcarry calcs[0]がsum
grover.x(calcs[0])
grover.mct(calcs, oracle)
grover.barrier()

grover.measure(inputs, c[:2])
grover.measure(oracle, c[2])

f:id:sakumadaisuke32:20210411165525p:plain
図3. a + b = 2 の解を検出する場合

この測定結果は下記です。

f:id:sakumadaisuke32:20210411165633p:plain
図4. 図3の回路の測定結果

この時点ではすべて等確率ですが、MCToracleビットが反転しているのは入力が11の場合(図4の一番右)のみです。

逆演算(Uncomputation)

満たす解を確認するだけだったらここまでで良いのですが、Groverのアルゴリズムを実行するには、逆演算(Uncomputation)という処理が必要です。
測定したい入力のビット(inputs)とオラクルのビット(oracle)以外のancilaとして用いたビット(ここではcaluculation)へ掛けたゲートの効果を「無かったことにする」処理です。
掛けた時と逆順にゲートを掛ければよいです。

inputs = QuantumRegister(2, 'input')
calcs = QuantumRegister(2, 'calculation')
oracle = QuantumRegister(1, 'oracle')
c = ClassicalRegister(2, 'output')

grover = QuantumCircuit(inputs, calcs, oracle, c)
grover.h(inputs)
grover.barrier()

grover = adder(grover, inputs[0], inputs[1], calcs[1], calcs[0]) # calcs[1]がcarry calcs[0]がsum
grover.x(calcs[0])
grover.mct(calcs, oracle)
grover.x(calcs[0])
grover = adder_r(grover, inputs[0], inputs[1], calcs[1], calcs[0])

f:id:sakumadaisuke32:20210411171433p:plain
図5. 逆演算をかけた場合

逆演算そのものの数理的な解説は「量子コンピューティング 量子アルゴリズムから機械学習まで(嶋田義皓著)」の2.5節を参照してください。 ここでは実際に逆演算を仮にしないまま位相反転増幅した場合にうまくいかないことを示します。

計算する回路全体の構成は下記です。

inputs = QuantumRegister(2, 'input')
calcs = QuantumRegister(2, 'calculation')
oracle = QuantumRegister(1, 'oracle')
outputs = ClassicalRegister(2, 'output')
grover = QuantumCircuit(inputs, calcs, oracle, outputs)

grover.h(inputs)
grover.x(oracle)
grover.h(oracle)
grover.barrier()

grover = adder(grover, inputs[0], inputs[1], calcs[1], calcs[0]) # calcs[1]がcarry calcs[0]がsum

grover.x(calcs[0])
grover.mct(calcs, oracle)
grover.barrier()

grover.h(oracle)
grover.x(oracle)

grover.h(inputs)
grover.x(inputs)
grover.barrier(inputs)
grover.h(inputs[1])
grover.cx(inputs[0], inputs[1])
grover.h(inputs[1])
grover.barrier(inputs)
grover.x(inputs)
grover.h(inputs)

grover.measure(inputs, outputs)

f:id:sakumadaisuke32:20210417121428p:plain
図6. 逆演算しなかった場合

図6の回路でマーキングした直後の量子状態は下記です。(計算はご自身でやってみてください。)


\begin{aligned}
|absco \gt &= \frac{1}{2} ( |00 \gt_{ab} |10 \gt_{sc} + |01 \gt_{ab} |00 \gt_{sc} + |10 \gt_{ab} |00 \gt_{sc} - |11 \gt_{ab} |11 \gt_{sc}) \\
&\otimes \frac{1}{\sqrt{2}} ( |0 \gt_{o} - |1\gt_{o})
\end{aligned}

ここで、 aはinput0,  bはinput1,  sはsum(calculation0),  cはcarry(calculation1), そして ooracleの各ビットを示します。 求めたい入力である

\displaystyle{
|11 \gt_{ab} |11 \gt_{sc}
}

の状態がマーキングされています。

そして、この後位相反転増幅した後の最終状態は下記です。


\begin{aligned}
|absco \gt &= \frac{1}{4} \{ (- |00 \gt_{ab} + |01 \gt_{ab} + |10 \gt_{ab} + |11 \gt_{ab}) |10 \gt_{sc} \\
&+ 2(|00 \gt_{ab} + |11 \gt_{ab}) |00 \gt_{sc} \\
&-2 |00 \gt_{ab} |11 \gt_{sc} \} \otimes |0 \gt_{o}
\end{aligned}

sumとcarryとのエンタングルメントを解除していないので、各項が分かれたままきれいにまとまらなくなってしまっています。 これを |ab \gtを中心に整理しなおすと


\begin{aligned}
|absco \gt &= \frac{1}{4} \{ (|00 \gt_{ab} (- |10 \gt_{sc} + 2|00 \gt_{sc} - 2 |11 \gt_{sc}) \\
&+ |01 \gt_{ab} |10 \gt_{sc} \\
&+ |10 \gt_{ab} |10 \gt_{sc}\\
&+ |11 \gt_{ab} (|10 \gt_{sc} + 2 |00 \gt_{sc}) \} \otimes |0 \gt_{o} \\
\end{aligned}

となり、[tex: |00 \gt{ab}, |11 \gt{ab}]の確率振幅の大きさが[tex: |01 \gt{ab}, |10 \gt{ab}]のものより少し大きくなるということになります。

実際に測定すると下図となり、推測した通りの傾向になり、そして解を得ることができません。

f:id:sakumadaisuke32:20210417124847p:plain
図7. 逆演算しなかった場合の確率分布

逆演算の必要性を理解していただけると幸いです。

計算結果

逆演算と位相反転増幅も含めた全体の回路は下記です。

a + b = 0を満たす解を求める

inputs = QuantumRegister(2, 'input')
calcs = QuantumRegister(2, 'calculation')
oracle = QuantumRegister(1, 'oracle')
outputs = ClassicalRegister(2, 'output')
grover = QuantumCircuit(inputs, calcs, oracle, outputs)

grover.h(inputs)
grover.x(oracle)
grover.h(oracle)
grover.barrier()

grover = adder(grover, inputs[0], inputs[1], calcs[1], calcs[0]) # calcs[1]がcarry calcs[0]がsum

grover.x(calcs)
grover.mct(calcs, oracle)
grover.x(calcs)

grover = adder_r(grover, inputs[0], inputs[1], calcs[1], calcs[0])
grover.barrier()

grover.h(oracle)
grover.x(oracle)

grover.h(inputs)
grover.x(inputs)
grover.barrier(inputs)
grover.h(inputs[1])
grover.cx(inputs[0], inputs[1])
grover.h(inputs[1])
grover.barrier(inputs)
grover.x(inputs)
grover.h(inputs)

grover.measure(inputs, outputs)

f:id:sakumadaisuke32:20210411172142p:plain
図8. a + b = 0を満たす解を求める回路

この測定結果は下図で、確かに | 00 \gt が増幅されています。

f:id:sakumadaisuke32:20210411172252p:plain
図9. 図8の回路の測定結果

a + b = 2を満たす解を求める

inputs = QuantumRegister(2, 'input')
calcs = QuantumRegister(2, 'calculation')
oracle = QuantumRegister(1, 'oracle')
outputs = ClassicalRegister(2, 'output')
grover = QuantumCircuit(inputs, calcs, oracle, outputs)

grover.h(inputs)
grover.x(oracle)
grover.h(oracle)
grover.barrier()

grover = adder(grover, inputs[0], inputs[1], calcs[1], calcs[0]) # calcs[1]がcarry calcs[0]がsum

grover.x(calcs[0])
grover.mct(calcs, oracle)
grover.x(calcs[0])

grover = adder_r(grover, inputs[0], inputs[1], calcs[1], calcs[0])
grover.barrier()

grover.h(oracle)
grover.x(oracle)

grover.h(inputs)
grover.x(inputs)
grover.barrier(inputs)
grover.h(inputs[1])
grover.cx(inputs[0], inputs[1])
grover.h(inputs[1])
grover.barrier(inputs)
grover.x(inputs)
grover.h(inputs)

grover.measure(inputs, outputs)

f:id:sakumadaisuke32:20210411172558p:plain
図10. a + b = 2を満たす解を求める回路

この測定結果は下図で、確かに | 11 \gt が増幅されています。

f:id:sakumadaisuke32:20210411172647p:plain
図11. 図10の回路の測定結果

まとめ

  • ラクル回路には求めたい解の条件を埋め込む
  • ラクル回路でマーキングした後は逆演算でエンタングルメントを解除する
  • 位相反転増幅で入力ビットの振幅を増幅

これでひとまずグローバーアルゴリズムの一連の解説を終わります。 オラクル回路を問題に合わせて埋め込めことができれば、解ける問題の数・種類が広がります。 良い量子プログラミングライフを!


量子コンピューティング 基本アルゴリズムから量子機械学習まで [ 情報処理学会出版委員会 ]

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

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

シリーズ予定:

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

対象とする読者

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

前提とする知識

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

バージョン情報

目次

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の状態が一番確率振幅の大きさが大きくなっています。

qasm_simulatorで動作確認

図2の回路で、測定をした場合の確率分布が下記です。

f:id:sakumadaisuke32:20210314110925p:plain
図6. 3量子ビットGrover回路の測定結果

まとめ

グローバーのアルゴリズム基礎(2):位相反転増幅回路

本日は前回の続きで、位相反転増幅回路について解説します。

シリーズ予定:

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

前回はオラクル回路によって正解の状態の位相(正負)を反転させてマーキング(目印)しました。
ただ、この段階では測定しても確率振幅の大きさ自体は変わっていないので、他の正解でない状態の確率とは等確率のため見分けがつきません。
今回の位相反転増幅回路はマーキングした状態の確率振幅の大きさを大きくし、他の状態の確率振幅の大きさを相対的に低くする手法です。

(余談ですが、確率振幅「の大きさ」と言っているのは、確率振幅は複素数だからです。
通常の水波・地震波・振り子などの振幅は実数ですから「の大きさ」は不要です。)

対象とする読者

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

前提とする知識

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

バージョン情報

目次

位相反転増幅回路の役割と仕組み

役割と仕組み

Groverのアルゴリズムの概略は下記でした。

f:id:sakumadaisuke32:20210214150611p:plain
図1. Grover回路の概要

前回のオラクル回路はこの図の U_fで、今回は右側の H^{\otimes n} (2 | 0 ^n \gt \lt 0 ^n | - I_n) H^{\otimes n}の「位相反転増幅回路」について解説します。

この回路の役割はオラクル回路でマーキングした状態に量子状態を近づけることです。

f:id:sakumadaisuke32:20210227161153p:plain
図2. 位相反転増幅の概念図. マーキングした状態 |\Psi_{t'} \gtから元の状態 | s \gtを軸に反転する.(QuantumChallenge2020の問題より)

ラクル回路でターゲットの状態にマーキングした状態 |\Psi_{t'} \gtに対して、元の状態 | s \gtを軸に反転させる操作を行い、ターゲットの状態 | w \gtに近づけます。

この操作 U_sは数式で表すと下記です。

\displaystyle{
U_s |\Psi_{t'} \gt = (2 |s \gt \lt s | - I ) |\Psi_{t'} \gt
}

数式からの理解は嶋田さんの本に記載がありますのでここでは図形的に理解を図ります。
このユニタリー変換  U_s を日本語にすると、「ベクトル  | \Psi_{t'} \gt のうちベクトル  | s \gt に沿った成分を2倍して、反対向きのベクトルを加算する」なのですが図解すると、下図の通りひし形になります。
(4辺が等しい平行四辺形はひし形)

f:id:sakumadaisuke32:20210227170018p:plain
図3. ユニタリー変換 U_sの図解

これを最初に思いついた人はすごいとしか思えませんが、結果的にこの演算をすると、ベクトル | \Psi_{t'} \gt |s>を軸に反転することになります。
そして | s \gt | s' \gtの角度を \thetaとしたときに |s>から 2 \thetaだけ目的の状態に近づきます。

ラクル回路と位相反転増幅回路の組み合わせを適切な回数繰り返すことによって少しずつターゲットの状態 | w \gtに近づきます。
ただ、繰り返しすぎると通り過ぎてしまうのでちょうどいい回数で止めることに注意です。

実装例

この式の演算を量子ゲートに変換すると下記です。

qr = QuantumRegister(2, 'data')
cr = ClassicalRegister(2, 'output')
grover = QuantumCircuit(qr, cr)

# 重ね合わせ状態の生成
grover.h(qr)

# オラクル
grover.cz(qr[0], qr[1])

# 位相反転増幅
grover.h(qr)
grover.x(qr)
grover.h(qr[-1])
grover.cx(qr[0], qr[1])
grover.h(qr[-1])
grover.x(qr)
grover.h(qr)

# 測定
grover.measure(qr, cr)

f:id:sakumadaisuke32:20210214152132p:plain
図4. Grover回路の全体図. 位相反転増幅回路は緑枠で囲った箇所.

位相反転増幅回路の中身

位相反転増幅 U_sが何をしているのか分解していこうと思いますが、この式は次のように書き換えられます。
(簡単のために2量子ビット系にします。)


\begin{aligned}
U_s &= 2 |s >< s| - I \\
&= 2 H^{\otimes 2} | 00 >< 00 | H^{\dagger \otimes 2} - H^{\otimes 2}  I H^{\dagger \otimes 2} \\
&= H^{\otimes 2} (2 | 00 >< 00 | - I ) H^{\dagger \otimes 2} \\
&= H^{\otimes 2} (2 | 00 >< 00 | - I ) H^{\otimes 2}
\end{aligned}

両サイドのHゲートはそのままHゲートで良いのですが、間の  2 | 00 \gt \lt 00 | - I の実装方法がミソです。
図4の内 2 | 00 \gt \lt 00 | - I の操作に対応する箇所だけ抜き出したのが下記です。

f:id:sakumadaisuke32:20210227190437p:plain
図5.  2 | 00 \gt \lt 00 | - I に対応する回路

 2 | 00 \gt \lt 00 | - I の操作は  | 00 >を軸にそれ以外の状態の確率振幅の位相を反転する操作です。
(図3の | s \gt | 00 \gtに置き換えてください。)
この実装が図5と言われていますが、この回路そのものは「すべてのビットがゼロの状態(今回なら | 00 >)の位相だけ反転して、他の状態はそのまま」という操作をします。
まるで中身は逆ですがこれで合っています。
グローバル位相を考慮すれば、それ以外の確率振幅については等価な状態となっていることを以下で示します。

例えば、今回の2量子ビットの場合、Hゲート通過後でこの回路に入る直前の状態は下記です。

\displaystyle{
|s' > = \frac{1}{2} (|00> + |01> + |10> - |11>)
}

まず元の理論の確認をするため、ここに 2 | 00 \gt \lt 00 | - I を掛けます。


\begin{aligned}
(2 | 00 \gt \lt 00 | - I) |s' \gt &= (2 | 00 \gt \lt 00 | - I) \frac{1}{2} (|00 \gt + |01 \gt + |10 \gt - |11 \gt) \\
&= \frac{1}{2} (2 | 00 \gt - | 00 \gt - | 01 \gt - | 10 \gt + | 11 \gt ) \\
&= \frac{1}{2} (| 00 \gt - | 01 \gt - | 10 \gt + | 11 \gt ) \equiv |s'' \gt
\end{aligned}

一方、回路の挙動を確認するため図5の量子回路を通過した後の状態( | 00 \gt のみ位相を反転)は下記です。


\begin{aligned}
|s''' > &= \frac{1}{2}  (-|00> + |01> + |10> - |11>) \\
&= - \frac{1}{2}  (|00> - |01> - |10> + |11>) \equiv - |s'' \gt
\end{aligned}

ということで、グローバル位相は \piだけ異なりますが、他は同一であることが分かります。

この実装になっている理由は僕は分かりませんが、元々実現したかった「 | 00 >以外の確率振幅の位相を反転する」という操作を量子回路で実装する方法が思い当たりませんでした。
もしかしたらあるのかもしれませんが、難しいのかもしれません。

ちなみに、仮にどうにかして「 | 00 >以外の確率振幅の位相を反転する」という回路が実現できたとして、その後の量子状態は | s'' \gtです。

元々の数学上で行いたかった操作は下図の U_sですが、実装では U'_sのような変換をしているのだと思います。
結果的に目的の状態 | w \gtを取り出せているので良いです。

f:id:sakumadaisuke32:20210227192316p:plain
図6. 図5の量子回路が行っている操作 \hat{U'}_s

最後に再びHゲートを全体に掛けます。


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

ということで、マーキングした状態 |11\gtを取り出すことができました。
量子回路で実装した場合は位相が \pi違うだけで、測定して観測される状態は同一です。

まとめ

  • 位相反転増幅回路を通すことでオラクル回路でマーキングした増幅を取り出すことができる。
  • 位相反転増幅回路の内部では元々の数学上の操作とは異なるが、グローバル位相を考慮すれば入出力は目的を達成している。

グローバーのアルゴリズム基礎(1):オラクル回路

今回からは4回のシリーズでGrover探索アルゴリズムの解説をします。 Grover探索アルゴリズムもQuantumChallenge2020に出題され、仕組みと使い方を覚えると汎用性が高いアルゴリズムです。 QuantumChallenge2020公式

GroverのアルゴリズムはDeutsch–Jozsaアルゴリズム・Shorのアルゴリズムと並んで量子コンピュータの構想初期から注目され続けている量子アルゴリズムです。 公式をはじめとして様々なところで既に解説はされているのですが、数式と実装する量子回路との対応、どうしてそのゲートを使うのかなどの理由も深堀し直して初学者の人がイチからコーディングできるように、また応用を利かせられるような解説を目指します。

次の予定でいます。

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

今回は「オラクル回路」という、探索する状態の選択肢の中から「これが答えだ」「これを見つけたい」という状態を見つけてマーキングする回路です。 ちなみに「オラクル(Oracle)」とは日本語で「神託(神のお告げ)」(Wikipedia)ということで、与えられた問題の答えを予め知っているから、というのが言葉の由来です。

対象とする読者

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

前提とする知識

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

バージョン情報

目次

Groverのアルゴリズムの概要

概要の説明はQuantumChallenge2020問題の前半の解説(日本語)が絵も有ってまとまっていて良いと思います。

ここでは詳細な説明は割愛しますが、要素は次の2点です。

  • 答えに対応する状態をマーキング
  • そのマーキングされた状態の確率だけを増幅

最終的に状態を測定すると探していた状態が100%に近い確率で測定される、というものです。

簡略した量子回路は下記で書かれます。

f:id:sakumadaisuke32:20210214150611p:plain
図1. Grover回路の概要(QuantumChallenge2020の問題より)

この図の意味

2量子ビットの場合は下記が実装の1例になります。

f:id:sakumadaisuke32:20210214152132p:plain
図2. 2量子ビットGroverアルゴリズムの実装例

以降ではオラクル回路をどうして上図のように実装するのか、数式と量子回路の対応を確かめます。

ラクル回路の役割と仕組み

役割と仕組み

f:id:sakumadaisuke32:20210214152640p:plain
図3. オラクル回路の役割(QuantumChallenge2020の問題より)

左右の図について説明します。

 |w>が探している状態です。

 |s>が探索する全状態の重ね合わせ状態で、下式です。

\displaystyle{
|s> = \frac{1}{\sqrt{N}} \sum_{i=0}^{N-1}{|i>}
}

 Nが状態の数です。
上式は一般式ですが、
例えば2量子ビットの場合は取りうる全状態は  |00>, |01>, |10>, |11> の 22 = 4 個、
3量子ビットの場合は  |000>, |001>, |010>, |011>, |100>, |101>, |110>, |111> の 23 = 8 個です。
上式とは  |i=0> = |00>, |i=1> = |01>, |i=2> = |10>, |i=3> = |11> という風に対応します。

左の図では、 |w>に直交する( |w>との内積がゼロになる)状態を仮に |s'>としています。
( |s>の重ね合わせの内、 |w>の係数がゼロな状態が相当)
この状態(ベクトル)を軸に反転した状態を |\psi >とします。
この結果、 |\Psi _{t'}> |s>の重ね合わせの内、 |w>の係数だけマイナスになったものです。

この |\Psi _{t'} >の状態の確率振幅を図示したのが右の図です。
状態 |w>の確率振幅だけが他と絶対値は同じで符号が反転しています。

例えば2量子ビットの場合

\displaystyle{
|s> = \frac{1}{2} (|00> + |01> + |10> + |11>)
}

 |w> = |11>の時、

\displaystyle{
|\Psi_{t'}> = \frac{1}{2} (|00> + |01> + |10> - |11>)
}

です。

というように、全状態の重ね合わせ状態 |s>から探している状態 |w>の確率振幅の符号を反転(=位相を反転)させた |\Psi_{t'}>に変換することがオラクル回路の役割です。
またこの目的の状態の位相を反転することを「マーキングする」と言うことがあります。

ラクル回路の実装例

では数学的にやりたいことはわかったので、量子回路で実装しましょう。

やりたいこと:

  • 取りうるすべての状態を生成
  • 答えの状態の位相を反転するようにゲート操作

図2に対応する実装は下記です。

qr = QuantumRegister(2, 'data')
cr = ClassicalRegister(2, 'output')
grover = QuantumCircuit(qr, cr)

grover.h(qr)
grover.cz(qr[0], qr[1])
grover.measure(qr, cr)

grover.draw(output='mpl')

f:id:sakumadaisuke32:20210214182626p:plain
図4. 図2に対応するオラクル回路

この時の状態ベクトルの確率振幅を確認しましょう。
4つ目の状態だけ符号がマイナスになっていますね!
ということで、目指した状態をこの回路で作れることは分かりました。

qr = QuantumRegister(2, 'data')
cr = ClassicalRegister(2, 'output')
grover = QuantumCircuit(qr, cr)

grover.h(qr)
grover.cz(qr[0], qr[1])

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

state = result.get_statevector()
print(state)

>> array([ 0.5+0.j,  0.5+0.j,  0.5+0.j, -0.5+0.j])

ではこの時この回路では何が起きているのか、また |11>以外の状態の位相を反転させたいときにはどうしたらいいのか確認しましょう。

ラクル回路の中身

まずはHゲート(アダマールゲート)は取りうる全状態の重ね合わせを作ります。

qr = QuantumRegister(2, 'data')
cr = ClassicalRegister(2, 'output')
grover = QuantumCircuit(qr, cr)

grover.h(qr)

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

state = result.get_statevector()
print(state)
grover.draw(output='mpl')

>> [0.5+0.j 0.5+0.j 0.5+0.j 0.5+0.j]

f:id:sakumadaisuke32:20210214184817p:plain
図5. 全状態の重ね合わせの生成

f:id:sakumadaisuke32:20210214190136p:plain
図6. 図5の時のBloch球

Bloch球では2つのビットとも (1 / \sqrt{2}) (|0> + |1>)になっていて、また確率振幅も等しいですね。

つまり、下記の数式の状態になっています。


\begin{aligned}
H|0>_0 \otimes H|0>_1 &= \frac{1}{\sqrt{2}} (|0>_0 + |1>_0)\otimes \frac{1}{\sqrt{2}} (|0>_1 + |1>_1) \\
&= \frac{1}{2} (|00> + |01> + |10> + |11>) \\
&\equiv |s>
\end{aligned}

これはよくやるので多分当たり前になっていますが、「全状態の重ね合わせを作りたい」という目的に対して、Hゲートで1量子ビットの重ね合わせ状態を作り、それぞれの量子ビットを相互作用させず独立にしておくことで直積(総当たりの組み合わせ)で表現できるようにすることで目的を達成させている例です。

ではこの後のCZゲートが何しているのか確認しましょう。
そもそもCZゲートは制御ビットが |1>の時に標的ビットにZゲートを掛ける(位相を反転する)というものです。

制御ビットが |0>の場合
分かりやすくするため、標的ビットにはHゲートを掛けておきます。

qr = QuantumRegister(2, 'data')
qc = QuantumCircuit(qr)

qc.h(qr[1])
qc.cz(qr[0], qr[1])

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

state = result.get_statevector()
print(state)
qc.draw(output='mpl')
>> [ 0.70710678+0.j  0.        +0.j  0.70710678+0.j -0.        +0.j]

f:id:sakumadaisuke32:20210214191608p:plain
図7. 制御ビットが |0>の場合のCZ

f:id:sakumadaisuke32:20210214191740p:plain
図8. 図7の状態のBloch球

標的ビットの状態ベクトルに変化はないですね。
ちなみに確率振幅の配列は左から順に |00>, |10>, |01>, |11>です。
だから数式としては下記の状態になって、Zゲートがかからずそのままということです。


\begin{aligned}
|0>_0 \otimes H|0>_1 &\stackrel{\mathrm{CZ}}{\longrightarrow} |0>_0 \otimes \frac{1}{\sqrt{2}} \left\{ |0>_1 +  \exp^{i (\phi + 0 \dot \pi)}|1>_1 \right\} \\
&= \frac{1}{\sqrt{2}} (|00> + |01>)
\end{aligned}

制御ビットが |1>の場合

qr = QuantumRegister(2, 'data')
qc = QuantumCircuit(qr)

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

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

state = result.get_statevector()
print(state)
qc.draw(output='mpl')
>> [ 0.        +0.j  0.70710678+0.j  0.        +0.j -0.70710678+0.j]

f:id:sakumadaisuke32:20210214192456p:plain
図9. 制御ビットが |1>の時のCZ

f:id:sakumadaisuke32:20210214192559p:plain
図10. 図9の時のBloch球

標的ビットの状態ベクトルがZ軸を軸に180度回っています。
つまり、標的ビットの位相が反転したということです。
確率振幅も |11>だけ符号がマイナスになっています。


\begin{aligned}
|1>_0 \otimes H|0>_1 &\stackrel{\mathrm{CZ}}{\longrightarrow} |1>_0 \otimes \frac{1}{\sqrt{2}} \left\{ |0>_1 +  \exp^{i (\phi + 1 \dot \pi)}|1>_1 \right\} \\
&= \frac{1}{\sqrt{2}} (|10> - |11>)
\end{aligned}

そして制御ビットを重ね合わせにすると、4状態の内 |11>の符号だけ反転するということです。
ここではCZゲートの働きがやりたいことをちょうど満たしてくれているから使っているということです。

ラクル回路で符号を反転させる状態を変える

では例えば4状態の内、 |00>の符号を反転させるには?
CZゲートの後に制御ビット・標的ビットともにXゲートを掛ければ良いです。
イメージはCXゲートで制御ビットが|0>の時に標的ビットを反転させるのに近いです。

実装としてはこちらです。

qr = QuantumRegister(2, 'data')
qc = QuantumCircuit(qr)

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

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

state = result.get_statevector()
print(state)
qc.draw(output='mpl')
>> [-0.5+0.j  0.5+0.j  0.5+0.j  0.5+0.j]

f:id:sakumadaisuke32:20210214222108p:plain
図11.  |00>の符号を反転させるオラクル回路

数式上では


\begin{aligned}
&\frac{1}{2} \mathrm{X_0} (|0>_0 + |1>_0) \otimes \mathrm{X_1} (|0>_1 + |1>_1) = \frac{1}{2} (|0>_0 + |1>_0) \otimes (|0>_1 + |1>_1) \\
&\stackrel{\mathrm{CZ}}{\longrightarrow} \frac{1}{2} \left[|0>_0 \otimes \left\{ |0>_1 +  \exp^{i (0 \dot \pi)}|1>_1 \right\} + |1>_0 \otimes \left\{ |0>_1 +  \exp^{i (1 \dot \pi)}|1>_1 \right\}\right]\\
&= \frac{1}{2} (|00> + |01> + |10> - |11>) \\
&\stackrel{\mathrm{X^{\otimes 2}}}{\longrightarrow} (|11> + |10> + |01> - |00>)
\end{aligned}

もう一例、 |01>の符号を反転させたい場合は制御ビットだけCZの後にXゲートを掛ければ良いです。

qr = QuantumRegister(2, 'data')
qc = QuantumCircuit(qr)

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

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

state = result.get_statevector()
print(state)
qc.draw(output='mpl')
>> [ 0.5+0.j  0.5+0.j -0.5+0.j  0.5+0.j]

f:id:sakumadaisuke32:20210214222211p:plain
図12.  |01>の符号を反転させるオラクル回路

このCZゲートを使ったオラクル回路についてはデフォルトでは |11> (一般にはすべてのビットが1の状態)の符号が反転し、状態を |0>にするビットについてCZゲートの後にXゲートを掛ければよいです。

他にCXゲートを使ったオラクル回路がありますが、その解説はまたの機会に譲ります。

まとめ

今回はGroverのアルゴリズムの要素の一つであるオラクル回路の実装の一つでCZゲートを使ったものを解説しました。
CZゲートの制御ビットが1の時に標的ビットの位相を反転する性質がオラクル回路の実装に都合がいいということを示しました。
他の実装にはCXゲートを使ったものがあります、更にそれ以外があるかはわかりませんが、都合がいいゲートを組み合わせてプログラミングをしましょう。

参考

Groverのアルゴリズムの解説ページ