PennyLaneを用いた量子回路最適化
僕が学生時代に量子情報処理の研究に関わっていたのはこちらで自己紹介した通りで、
研究者ではありませんが今も量子情報処理・量子計算には興味があります。
そして今は機械学習エンジニア・データサイエンティストとしてソフトウェア開発に携わっています。
そんな自分が今量子コンピュータの界隈でできたらいいな・機会があれば取り組みたい・自分の経歴を生かすことができると思っていることが、量子機械学習です。
量子機械学習についてはQiitaでKuma氏がまとめてくださっています。
主に以下の目的で研究されています。
ただ量子計算も機械学習も両方わからないといけないのでなかなかハードです。
みんなが量子コンピュータ使える頃には使い方をちゃんと知っている人しか使いこなせないという参入障壁が高い状態になっていると困ります。
だから今のうちにやりたい、やりましょう!
ここからは量子機械学習を実装したフレームワークの一つであるPennyLaneのチュートリアルをベースに、要素技術を深掘りして解説していきます。
今回のシリーズは量子回路のパラメータ最適化(回転ゲートの回転角の調整)です。
こちらのチュートリアルです。
対象とする読者
前提とする知識
- 機械学習もしくは最適化という言葉を知っていること
- 量子ビット、波動関数、重ね合わせ状態、測定という言葉を知っていること
- 簡単な1量子ビットゲート(X, Y, Z, Rx, Ry, Rz)を理解していること
バージョン情報
- Python 3.9.5
- PennyLane 0.15.1
目次
PennyLaneとは
PannyLaneはカナダのXanado社が開発している量子機械学習を実装した(量子機械学習をするのに必要な・あったら便利な関数やクラスを用意してくれている)ライブラリの一つで、機械学習フレームワークであるPyTorchをベースにしています。
PyTorchで作られているので計算結果をPyTorchに連携できます。
ただし、現状では自分のPC内にシミュレータを構築して計算しているので、実機ではありません。
その他に量子機械学習を実装しているフレームワークは2021年6月現在把握してある限り以下です。
あまりないですね、、
中でもこのPennyLaneはチュートリアルが豊富で、ノートブックも用意されているのでとりあえず実行することができます!
コアとなる概念: 量子ノード(Quantum Node, QNode)
PennyLaneを使う上で重要な概念が量子ノード(Quantum Node, QNode)です。
下図のように量子ビットや量子ゲートのパラメータ(入力)、測定結果(出力)、そしてその間の何らかの量子回路からなる計算単位です。
このQNodeを一つの計算グラフとして計算結果を後続の処理に渡したり、組み合わせて量子ニューラルネットワークを構築するのでしょう。
今回のシリーズではQNodeは1つしか使わないので、入力から量子ゲートを経て出力までをセットにした関数と思っていただければ良いです。
量子回路の計算結果を目的関数として最小化問題を解く
古典の機械学習では解きたい課題を表す目的関数を定義し、様々なパラメータでその目的関数を計算し、目的関数の値が最小もしくは最大となるパラメータの組み合わせを見つけます。
量子機械学習でも同様に、量子回路での計算結果を使って目的関数を定義し、量子回路のパラメータを色々変えて目的関数を計算し、目的関数の値が最小もしくは最大となるパラメータの組み合わせを見つけます。
通常、TensorFlowやPyTorchなどの機械学習フレームワークでは最小とするように計算するようになっています。
今回の解く問題と量子回路
今回解こうとしている問題は「1つの量子ビットの状態をからに遷移させるための, ゲートの角度はいくらか」(図2)です。
これは明らかにXゲートもしくはYゲートをかければ良いのですが、あえてこれを探索的に見つけようという練習課題です。
量子回路は下記です。
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⟩
可視化のためにゲートとゲートの角度を指定していますが、これをこの後で探索します。
なお、角度の単位はradianです。
今回の誤差関数
PennyLaneでは測定にPauli演算子を使います。
今回はZ軸への射影測定を想定してPauliZ測定します。
演算子の期待値は状態については1, については-1であることを利用して、下記となります。
上で定義した量子回路そのままです。
def cost(x): return circuit(x)
誤差関数を数式で表すと、演算子を任意の状態で期待値を求めるということなので、細かい計算は割愛しますが下記になります。
の時に最大で, もしくはの時に最小で となります。
勾配降下法で最適化
勾配降下法で最適化します。
コードとしては以下です。
# 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))
やることは古典の機械学習と同じです。
誤差関数の測定値と勾配に基づいて、パラメータ更新をします。
誤差関数の勾配は下記です。
計算したらわかりますがの時に勾配もゼロとなりパラメータ更新されなくなってしまうので、初期値を少しずらしてあるのがミソです。
パラメータ設定と誤差測定を繰り返して数値的に最適なパラメータを探索する
上のコードを実行すると下記のように誤差関数の更新が行われて、最初に期待した通りにY軸周りに回転する場合が解として求まりました。
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で定義する
- 問題設定・誤差関数を量子回路で表現する
- パラメータ最適化は古典機械学習と同じく誤差関数の勾配から勾配降下法で行える
赤池情報量基準(AIC)を多項式の回帰で実装
こんにちは、さくまっちょです。
今回は本業のデータサイエンス関連で気になっていたこととして赤池情報量基準(AIC)を実際に実験用にデータを作って計算し その意味や威力を確認してみました。
機械学習・データサイエンスやっているとモデルの汎化性能を上げるために複雑なモデルを使わない方が良いという説明が
PRMLやはじパタなど各種機械学習系のテキストには必ずと言っていいほど載っています。
数式もちゃんと載っているのですが、いざ自分で計算してみようと思っても実装をしている記事などがなかなか出てこず、仕方がないので数式を調べながら実装してみました。
同じ疑問をお持ちの方のお役に立てば幸いです。
もしよければコメントいただきたいのですが普段みなさんってAICって計算しないのでしょうか?
もしくはいくつかのライブラリや計算クラスではAICやそれに準ずるモデル選択のメソッドが用意されているから自分では実装しないのでしょうか?
すべての実装についてはこちらの僕のノートブックをご覧ください。
対象とする読者
前提とする知識
- 回帰分析がわかること
- 基本的な数学(sin, logや線形代数)
バージョン情報
- Python 3.7.3
目次
赤池情報量基準(AIC)とは
Wikipediaによると
- 元統計数理研究所所長の赤池弘次が1971年に考案し1973年に発表した
- 統計モデルの良さを評価するための指標
- 「モデルの複雑さと、データとの適合度とのバランスを取る」ために使用される
- AIC最小のモデルを選択すれば、多くの場合、良いモデルが選択できる
です。
また、詳細な意味や議論については赤池氏の記事をぜひ読んでみてください。
よく比較されるBICとの違いについてもご本人が議論されています。
詳細な議論はともかくとして、AICを使うことによって様々に作ったモデルの中から「これが一番データに合ってそうだ」というモデルを決めるのに役立ちます。
多項式回帰での次数と汎化性能
PRMLやはじパタで説明されているように、僕もまずはデータを生成し多項式で回帰して汎化性能を確かめてみました。
作成したデータがこちらです。
これに0次から20次までの多項式で回帰した結果がこちらです。
6次くらいまでは正弦波っぽく描かれていますが、7次以降は赤点(訓練データ)との距離は縮んでそうですが歪みがあるのと青点(検証データ)との距離は離れてそうです。
実際に訓練データと検証データの回帰線とのRMSEを比較したのがこちらです。
3次までは訓練・検証ともRMSEは下がっていますね。
それ以降も横ばいな感じはありますが、訓練データは10次まで、検証データも9次まではRMSEが改善してそうです。
このグラフから得られる結果とグラフの形状から、なんとなく3次~6次の間でモデルを決定したいです。
ここで「なんとなく」ではなく定量的に次数を決めるための指標がAICです。
データの分布が正規分布を仮定した場合のAIC
AICの定義
Wikipediaによると
は尤度、はパラメータの数(今回は多項式の次数 + 1(定数項を含むため)) です。
尤度を計算するのに確率分布を導入しないといけなくて、このあたりの実装が厄介だったのですがこちらの41ページで正規分布の場合は簡単に整理できるという資料がありましたので、この通りに計算してみたのが次です。
データの分布が正規分布を仮定した場合のAIC
はデータの個数、は平均自乗誤差です。 実装はこちらになります。
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の式を用いて各次数について計算してみて最適次数を決定しました。
検証データのAICが最小になる3次にするのがよさそうと分かりました。
以上でAICの実装と使い方、そして実際に次数を決められることを示しました。
まとめ
qRAMの作り方
qRAM(Quantum Random Access Memory)を量子回路上で実装する方法を解説します。 こちらとこちらの論文をベースにエッセンスを解説し、これがあることで解ける問題の幅が広がることを紹介します。
対象とする読者
- 量子プログラミングに興味を持ち始めたばかりのエンジニア
- 量子計算に興味を持ち始めたばかりの高校生・大学生 など
前提とする知識
- 簡単な論理演算(AND, OR, NOT)が分かること
- 基本的な量子ゲート(X, CX, CCX)が分かること
バージョン情報
- Python 3.7.3
- Qiskit 0.23.1
目次
- qRAM(Quantum Random Access Memory)とは
- RAM(Random Access Memory)って何?
- 古典RAMと量子RAM(qRAM)の違い
- qRAMがあると何ができる?
- QRAMを実装するときの考え方
- QRAMの実装例
- まとめ
- 参考
qRAM(Quantum Random Access Memory)とは
定義
右辺に注目してもらえればいいのですが、入力の状態に対してデータがエンタングルして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')
ちゃんとエンタングルしているか測定して確認します。
状態を示す0/1の並びの左4桁がデータで右2桁が入力状態です。 上の表の通りにエンタングルしていますね。
QuantumChallenge2020の例
その他のqRAMの応用例としてQuantumChallenge2020のコードを紹介します。
まとめ
参考
グローバーのアルゴリズム基礎(4):まとめと簡単な応用例
すみません、忙しくて更新が遅くなりました。
このシリーズ最後はGroverのアルゴリズムの簡単な応用を紹介します。
(自分で考えたトイプロブレムです。)
シリーズ予定:
- オラクル回路
- 位相反転増幅回路
- 3量子ビット以上の場合
- まとめと簡単な応用例 <<今回
前回までで入力の量子ビットについて、特定の状態について振幅を増幅できるようになりました。
ただこれだけだと入力Nビットの何番目のビットが立った状態を増幅すればいいというだけなので、実行する前から答えがわかっている計算しかできません。
より実用に近づけるには、「入力Nビットについてある条件を満たすビットを知りたいのだけど、予めわからないのでGroverのアルゴリズムで探索・特定したい」という課題にアプローチできるようになる必要があります。
その一例はQuantumChallenge2020の出題の通りですが、これはQRAMとかやや複雑な数理があるので、今回は単純な問題設定にして解説します。
今回の問題設定は「a + b = 0もしくは2
の解となるa
, b
をGroverのアルゴリズムで見つける」です。
(= 1
となる解は増幅されないのでやりません。
理由はわからないのですが僕の理解が追いついていないからだと思うので、わかったらまたポストします。
もしくはご存じの方がいたらコメントください、、。)
対象とする読者
- 量子プログラミング・量子計算に興味を持って勉強し始めたエンジニア・高校生・大学生
- Groverのアルゴリズムについて改めて知りたい方 など
前提とする知識
- 簡単な論理演算(AND, OR, NOT)が分かること
- 基本的な量子ゲート(H, X, CX, CCX, Z, CZ)が分かること
- 線形代数(行列・ベクトルの積、直積、あとできれば線形変換)が分かること
- ベクトルのブラケット記法(など、もしくはDirac記法ともいう)が分かること
- 加算器の仕組みがわかること(わからなくとも入出力の結果を納得できること)
バージョン情報
- Python 3.7.3
- Qiskit 0.23.1
目次
問題設定
問題は「a + b = 0もしくは2
の解となるa
, b
をGroverのアルゴリズムで見つける」です。
この場合の答えはa = 0, b = 0
、a = 1, b = 1
で人間にとってはすぐわかりますが、これをGroverのアルゴリズムで解きます。
なお、a + b = 1
の解はa = 0, b = 1
もしくはa = 1, b = 0
でこちらも人間にとってはすぐわかりますが、Gorverのアルゴリズムで増幅されないので今回は対象外とします。
理由はわかってないのですが、複数の解が含まれる場合は増幅されないのかもしれません。
数理的に理解ができたら改めて解説しようと思います。
オラクル回路の構築
オラクル回路の演算子の確認
オラクル回路は見つけ出したい量子状態をマークして、その状態の位相を反転する回路でした。 一般に量子オラクルはマーキングしたい状態をとし、その個数をとすると下記のように表します。
例えばa + b = 0
を満たす入力はa = 0, b = 0
なので状態としてはを、
またa + b = 2
を満たす入力はa = 1, b = 1
なので状態としてはをマーキング(位相を反転)するような回路を構成できれば良いということです。
条件を満たす量子状態の振幅を反転
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])
この測定結果は下記です。
この時点ではすべて等確率ですが、MCTでoracle
ビットが反転しているのは入力が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])
この測定結果は下記です。
この時点ではすべて等確率ですが、MCTでoracle
ビットが反転しているのは入力が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])
逆演算そのものの数理的な解説は「量子コンピューティング 量子アルゴリズムから機械学習まで(嶋田義皓著)」の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)
図6の回路でマーキングした直後の量子状態は下記です。(計算はご自身でやってみてください。)
ここで、はinput0, はinput1, はsum(calculation0), はcarry(calculation1), そしてはoracleの各ビットを示します。 求めたい入力である
の状態がマーキングされています。
そして、この後位相反転増幅した後の最終状態は下記です。
sumとcarryとのエンタングルメントを解除していないので、各項が分かれたままきれいにまとまらなくなってしまっています。 これをを中心に整理しなおすと
となり、[tex: |00 \gt{ab}, |11 \gt{ab}]の確率振幅の大きさが[tex: |01 \gt{ab}, |10 \gt{ab}]のものより少し大きくなるということになります。
実際に測定すると下図となり、推測した通りの傾向になり、そして解を得ることができません。
逆演算の必要性を理解していただけると幸いです。
計算結果
逆演算と位相反転増幅も含めた全体の回路は下記です。
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)
この測定結果は下図で、確かにが増幅されています。
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)
この測定結果は下図で、確かにが増幅されています。
まとめ
これでひとまずグローバーのアルゴリズムの一連の解説を終わります。 オラクル回路を問題に合わせて埋め込めことができれば、解ける問題の数・種類が広がります。 良い量子プログラミングライフを!
グローバーのアルゴリズム基礎(3):3量子ビット以上の場合
前回まででグローバーのアルゴリズムの基本要素である「オラクル回路」「位相反転増幅回路」の理論と実装の対応を解説しました。
前回までは2量子ビットの場合でしたが、一般に3量子ビットの場合の実装方法を解説します。
シリーズ予定:
対象とする読者
- 量子プログラミング・量子計算に興味を持って勉強し始めたエンジニア・高校生・大学生
- Groverのアルゴリズムについて改めて知りたい方 など
前提とする知識
- 簡単な論理演算(AND, OR, NOT)が分かること
- 基本的な量子ゲート(H, X, CX, CCX, Z, CZ)が分かること
- 線形代数(行列・ベクトルの積、直積、あとできれば線形変換)が分かること
- ベクトルのブラケット記法(など、もしくはDirac記法ともいう)が分かること
バージョン情報
- Python 3.7.3
- Qiskit 0.23.1
目次
Groverのアルゴリズムのおさらい
オラクル回路(量子オラクル)
オラクル回路は見つけ出したい量子状態をマークして、その状態の位相を反転する回路でした。 一般に量子オラクルはマーキングしたい状態をとし、その個数をとすると下記のように表します。
2量子ビットの場合での場合、この量子オラクルを通過後の状態は下記になります。
位相反転増幅回路
マーキングした状態の確率振幅の大きさを増幅する位相反転増幅は下記です。
2量子ビットの場合でがマーキングした状態とすると、位相反転増幅後は
ということでマーキングした状態を取り出すことができました。
2量子ビットの場合の量子回路
上記のオラクル回路・位相反転増幅回路を組み合わせたGrover回路の一例は下図です。
3量子ビット以上の場合に生じる課題
結論から説明すると、オラクル回路が現状のCZゲート
による実装方法では3量子ビット以上で困ります。
3量子ビットでは初期状態(オラクル回路への入力)が次の式のようになります。
例えばがマーキングしたい状態だとして、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')
※4つ目の量子ビット('oracle'と名付けたはオラクル回路のための補助ビットで、そのあとの演算では用いられません。
以降、N量子ビットのGroverのアルゴリズムを実装する場合はMCTゲートを下に伸ばしていく形で、最低でも(N + 1)量子ビットあれば実装可能です。
オラクル回路が正しく動作していることの確認
図2の量子回路が正しくGroverのアルゴリズムの実装となっていることを確かめます。
まず、オラクル回路の入力状態は下記です。
第4量子ビット(一番右)は図2の補助ビット('oracle')です。
そして、マーキングしたい状態が (対応する入力は)の場合、オラクル回路の出力は
でマーキング状態だけ位相が反転するようになっています。
シミュレータで動作を確認します。
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')
この状態ベクトルの確率振幅のシミュレーション結果(状態: 確率振幅)は下記です。
(左から順に第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
だけ大きさが負になっており、目的の回路ができています。
位相反転増幅回路が正しく動作していることの確認
位相反転増幅回路の入力状態は下記です。
第4量子ビットには何も演算されないので位相反転増幅した結果は下記となります。
というように、マーキングしたの確率振幅が大きくなっています。
また、その大きさが規格化されていることも分かります。
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')
この状態ベクトルの確率振幅のシミュレーション結果(状態: 確率振幅)は下記です。 (左から順に第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の回路で、測定をした場合の確率分布が下記です。
まとめ
グローバーのアルゴリズム基礎(2):位相反転増幅回路
本日は前回の続きで、位相反転増幅回路について解説します。
シリーズ予定:
前回はオラクル回路によって正解の状態の位相(正負)を反転させてマーキング(目印)しました。
ただ、この段階では測定しても確率振幅の大きさ自体は変わっていないので、他の正解でない状態の確率とは等確率のため見分けがつきません。
今回の位相反転増幅回路はマーキングした状態の確率振幅の大きさを大きくし、他の状態の確率振幅の大きさを相対的に低くする手法です。
(余談ですが、確率振幅「の大きさ」と言っているのは、確率振幅は複素数だからです。
通常の水波・地震波・振り子などの振幅は実数ですから「の大きさ」は不要です。)
対象とする読者
- 量子プログラミング・量子計算に興味を持って勉強し始めたエンジニア・高校生・大学生
- Groverのアルゴリズムについて改めて知りたい方 など
前提とする知識
- 簡単な論理演算(AND, OR, NOT)が分かること
- 基本的な量子ゲート(H, X, CX, CCX, Z, CZ)が分かること
- 線形代数(行列・ベクトルの積、直積、あとできれば線形変換)が分かること
- ベクトルのブラケット記法(など、もしくはDirac記法ともいう)が分かること
バージョン情報
- Python 3.7.3
- Qiskit 0.23.1
目次
位相反転増幅回路の役割と仕組み
役割と仕組み
Groverのアルゴリズムの概略は下記でした。
前回のオラクル回路はこの図ので、今回は右側のの「位相反転増幅回路」について解説します。
この回路の役割はオラクル回路でマーキングした状態に量子状態を近づけることです。
オラクル回路でターゲットの状態にマーキングした状態に対して、元の状態を軸に反転させる操作を行い、ターゲットの状態に近づけます。
この操作は数式で表すと下記です。
数式からの理解は嶋田さんの本に記載がありますのでここでは図形的に理解を図ります。
このユニタリー変換 を日本語にすると、「ベクトル のうちベクトル に沿った成分を2倍して、反対向きのベクトルを加算する」なのですが図解すると、下図の通りひし形になります。
(4辺が等しい平行四辺形はひし形)
これを最初に思いついた人はすごいとしか思えませんが、結果的にこの演算をすると、ベクトルをを軸に反転することになります。
そしてとの角度をとしたときにからだけ目的の状態に近づきます。
オラクル回路と位相反転増幅回路の組み合わせを適切な回数繰り返すことによって少しずつターゲットの状態に近づきます。
ただ、繰り返しすぎると通り過ぎてしまうのでちょうどいい回数で止めることに注意です。
実装例
この式の演算を量子ゲートに変換すると下記です。
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)
位相反転増幅回路の中身
位相反転増幅が何をしているのか分解していこうと思いますが、この式は次のように書き換えられます。
(簡単のために2量子ビット系にします。)
両サイドのHゲートはそのままHゲートで良いのですが、間の の実装方法がミソです。
図4の内の操作に対応する箇所だけ抜き出したのが下記です。
の操作は を軸にそれ以外の状態の確率振幅の位相を反転する操作です。
(図3のをに置き換えてください。)
この実装が図5と言われていますが、この回路そのものは「すべてのビットがゼロの状態(今回なら)の位相だけ反転して、他の状態はそのまま」という操作をします。
まるで中身は逆ですがこれで合っています。
グローバル位相を考慮すれば、それ以外の確率振幅については等価な状態となっていることを以下で示します。
例えば、今回の2量子ビットの場合、Hゲート通過後でこの回路に入る直前の状態は下記です。
まず元の理論の確認をするため、ここに を掛けます。
一方、回路の挙動を確認するため図5の量子回路を通過した後の状態( のみ位相を反転)は下記です。
ということで、グローバル位相はだけ異なりますが、他は同一であることが分かります。
この実装になっている理由は僕は分かりませんが、元々実現したかった「以外の確率振幅の位相を反転する」という操作を量子回路で実装する方法が思い当たりませんでした。
もしかしたらあるのかもしれませんが、難しいのかもしれません。
ちなみに、仮にどうにかして「以外の確率振幅の位相を反転する」という回路が実現できたとして、その後の量子状態はです。
元々の数学上で行いたかった操作は下図のですが、実装ではのような変換をしているのだと思います。
結果的に目的の状態を取り出せているので良いです。
最後に再びHゲートを全体に掛けます。
ということで、マーキングした状態を取り出すことができました。
量子回路で実装した場合は位相が違うだけで、測定して観測される状態は同一です。
まとめ
- 位相反転増幅回路を通すことでオラクル回路でマーキングした増幅を取り出すことができる。
- 位相反転増幅回路の内部では元々の数学上の操作とは異なるが、グローバル位相を考慮すれば入出力は目的を達成している。
グローバーのアルゴリズム基礎(1):オラクル回路
今回からは4回のシリーズでGrover探索アルゴリズムの解説をします。 Grover探索アルゴリズムもQuantumChallenge2020に出題され、仕組みと使い方を覚えると汎用性が高いアルゴリズムです。 QuantumChallenge2020公式
GroverのアルゴリズムはDeutsch–Jozsaアルゴリズム・Shorのアルゴリズムと並んで量子コンピュータの構想初期から注目され続けている量子アルゴリズムです。 公式をはじめとして様々なところで既に解説はされているのですが、数式と実装する量子回路との対応、どうしてそのゲートを使うのかなどの理由も深堀し直して初学者の人がイチからコーディングできるように、また応用を利かせられるような解説を目指します。
次の予定でいます。
今回は「オラクル回路」という、探索する状態の選択肢の中から「これが答えだ」「これを見つけたい」という状態を見つけてマーキングする回路です。 ちなみに「オラクル(Oracle)」とは日本語で「神託(神のお告げ)」(Wikipedia)ということで、与えられた問題の答えを予め知っているから、というのが言葉の由来です。
対象とする読者
- 量子プログラミング・量子計算に興味を持って勉強し始めたエンジニア・高校生・大学生
- Groverのアルゴリズムについて改めて知りたい方 など
前提とする知識
- 簡単な論理演算(AND, OR, NOT)が分かること
- 基本的な量子ゲート(H, X, CX, CCX, Z, CZ)が分かること
- 線形代数(行列・ベクトルの積、直積、あとできれば線形変換)が分かること
- ベクトルのブラケット記法(など、もしくはDirac記法ともいう)が分かること
バージョン情報
- Python 3.7.3
- Qiskit 0.23.1
目次
Groverのアルゴリズムの概要
概要の説明はQuantumChallenge2020問題の前半の解説(日本語)が絵も有ってまとまっていて良いと思います。
ここでは詳細な説明は割愛しますが、要素は次の2点です。
- 答えに対応する状態をマーキング
- そのマーキングされた状態の確率だけを増幅
最終的に状態を測定すると探していた状態が100%に近い確率で測定される、というものです。
簡略した量子回路は下記で書かれます。
この図の意味
2量子ビットの場合は下記が実装の1例になります。
以降ではオラクル回路をどうして上図のように実装するのか、数式と量子回路の対応を確かめます。
オラクル回路の役割と仕組み
役割と仕組み
左右の図について説明します。
が探している状態です。
が探索する全状態の重ね合わせ状態で、下式です。
が状態の数です。
上式は一般式ですが、
例えば2量子ビットの場合は取りうる全状態は の 22 = 4 個、
3量子ビットの場合は の 23 = 8 個です。
上式とは という風に対応します。
左の図では、に直交する(との内積がゼロになる)状態を仮にとしています。
(の重ね合わせの内、の係数がゼロな状態が相当)
この状態(ベクトル)を軸に反転した状態をとします。
この結果、はの重ね合わせの内、の係数だけマイナスになったものです。
このの状態の確率振幅を図示したのが右の図です。
状態の確率振幅だけが他と絶対値は同じで符号が反転しています。
例えば2量子ビットの場合
での時、
です。
というように、全状態の重ね合わせ状態から探している状態の確率振幅の符号を反転(=位相を反転)させたに変換することがオラクル回路の役割です。
またこの目的の状態の位相を反転することを「マーキングする」と言うことがあります。
オラクル回路の実装例
では数学的にやりたいことはわかったので、量子回路で実装しましょう。
やりたいこと:
- 取りうるすべての状態を生成
- 答えの状態の位相を反転するようにゲート操作
図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')
この時の状態ベクトルの確率振幅を確認しましょう。
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])
ではこの時この回路では何が起きているのか、また以外の状態の位相を反転させたいときにはどうしたらいいのか確認しましょう。
オラクル回路の中身
まずは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]
Bloch球では2つのビットともになっていて、また確率振幅も等しいですね。
つまり、下記の数式の状態になっています。
これはよくやるので多分当たり前になっていますが、「全状態の重ね合わせを作りたい」という目的に対して、Hゲートで1量子ビットの重ね合わせ状態を作り、それぞれの量子ビットを相互作用させず独立にしておくことで直積(総当たりの組み合わせ)で表現できるようにすることで目的を達成させている例です。
ではこの後のCZゲートが何しているのか確認しましょう。
そもそもCZゲートは制御ビットがの時に標的ビットにZゲートを掛ける(位相を反転する)というものです。
制御ビットがの場合
分かりやすくするため、標的ビットには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]
標的ビットの状態ベクトルに変化はないですね。
ちなみに確率振幅の配列は左から順にです。
だから数式としては下記の状態になって、Zゲートがかからずそのままということです。
制御ビットがの場合
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]
標的ビットの状態ベクトルがZ軸を軸に180度回っています。
つまり、標的ビットの位相が反転したということです。
確率振幅もだけ符号がマイナスになっています。
そして制御ビットを重ね合わせにすると、4状態の内の符号だけ反転するということです。
ここではCZゲートの働きがやりたいことをちょうど満たしてくれているから使っているということです。
オラクル回路で符号を反転させる状態を変える
では例えば4状態の内、の符号を反転させるには?
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]
数式上では
もう一例、の符号を反転させたい場合は制御ビットだけ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]
このCZゲートを使ったオラクル回路についてはデフォルトでは (一般にはすべてのビットが1の状態)の符号が反転し、状態をにするビットについてCZゲートの後にXゲートを掛ければよいです。
他にCXゲートを使ったオラクル回路がありますが、その解説はまたの機会に譲ります。
まとめ
今回はGroverのアルゴリズムの要素の一つであるオラクル回路の実装の一つでCZゲートを使ったものを解説しました。
CZゲートの制御ビットが1の時に標的ビットの位相を反転する性質がオラクル回路の実装に都合がいいということを示しました。
他の実装にはCXゲートを使ったものがあります、更にそれ以外があるかはわかりませんが、都合がいいゲートを組み合わせてプログラミングをしましょう。