量子コンピュータや量子アニーリングのイジング定式化を機械学習で自動生成

はじめに

研究者はずっと計算しているので計算に慣れてますが、一般人は定式化、二次関数と言われると蕁麻疹が出ます。量子アニーリング使いたいけど、定式化がめんどくさい、人材が獲得できないという人のために、定式化をデータから自動化する仕組みを考えてみたいと思います。

モチベーション

量子コンピュータでの量子断熱計算、QAOA、最近ではGrover Adaptive Searchなど、イジングモデルと呼ばれる物理モデルに社会問題を実装し、それを数理最適化していきます。また、量子アニーリングマシンでも定式化はとても課題のようです。

天才的に適当に定式化ができる人がいる一方で、そんなに都合よく量子アニーリングに人材を割けないという人もいると思うので、定式化を自動化してみます。

定式化

イジングモデルは、符号はプラス、マイナスありますが、

E = -\sum h_iq_i - \sum J_{ij}q_iq_j

の形に定式化します。この定式化は、

コスト関数
制約条件

にわけて設計して組み合わせます。コスト関数だけの問題、制約条件だけの問題、両方組み合わせる問題などもあります。

また、それらを組み合わせるときに調整変数というハイパーパラメータも出ますが、これも曲者ですが、その解決方法は今度紹介します。

定式化をデータで自動設計

パンがなければケーキを食べればいいじゃない。人がいなければAIで自動化すればいいじゃないですかということです。トライしてみます。正直アドリブで記事を書いてるので、制約条件が出るような難しい問題ではなく、簡単な問題で再現してみたいです。

適当な式を

まず３量子ビットの適当な式を別のイジング式に移してみます。

E = 2q_0 + q_2 -3 q_0 q_1

みたいな式を移してみます。

1.初期化

まずは移す先の式を用意します。こちらはランダムで係数を準備

E_{trial} = aq_0+bq_1+cq_2+dq_0q_1+eq_0q_2+fq_1q_2

となりました。課題は係数が多いことですかね。今後はこの問題を解決できるよう頑張りますが、今回は量子ビット数3に対して、3+3c2 = 3+3 = 6個の係数が出ます。

そして、ランダムでやります。

E_{trial} = 0.2q_0+0.4q_1+0.5q_2+0.2q_0q_1+0.3q_0q_2+0.4q_1q_2

はい。適当な式ができました。

2.最適化

SGDで最適化します。入力データは量子ビットの0か1。ラベルは移す元のエネルギー値です。最適化パラメータはaからfまでのイジング係数になります。

数値微分でSGDやってみます。

まずは、適当な入力データを作り、01を量子ビットに入力し、元のイジング式のエネルギー値を出します。

E = 2q_0 + q_2 -3 q_0 q_1

なので、001とかを入れると、

E=1

になります。

次に、移す先のイジング式にも01を入れます。

E_{trial} = 0.2q_0+0.4q_1+0.5q_2+0.2q_0q_1+0.3q_0q_2+0.4q_1q_2

に入れたものと、係数を一つ選びそれを微小変化させて数値微分を取ります。それを勾配として係数を更新してみましょう。

読み込み

ツールを読み込みます

import numpy as np
import matplotlib.pyplot as plt

コピー元のイジング式とランダムで用意したものを準備します。

#元のイジング式
E = np.array([2,-3,0,0,0,0,0,0,1])

#更新するもの
Et = np.array([0.2,0.2,0.3,0,0.4,0.4,0,0,0.5])

次に学習データを作ります。

量子ビットの数値配列を準備し、それに対応するエネルギーを準備しておきます。

#量子ビットの数字配列を準備し、それに対応するコストも求めておく
q = np.array([[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]])
Eq = np.zeros(8)

for i in range(8):
    Eq[i] = q[i]@E.reshape(3,3)@q[i]

Eq

array([ 0.,  1.,  0.,  1.,  2.,  3., -1.,  0.])

これは、こんな感じに

array([ 0.,  1.,  0.,  1.,  2.,  3., -1.,  0.])

SGD

数値微分を使ってSGD。f'(x) = (f(x+h) - f(x)) / h 学習率を決めて、更新をかける。
f(x)はイジングのエネルギーを求めて、ラベルとの差を損失関数で誤差を計算 x = x - e*f'(x)

必要なパラメータを準備します。

h = 0.01
e = 0.0001
rep = 10000
arr = []

Eu = np.zeros(9)

そして、SGDをかけて、損失関数を計算しながら更新をかけます。

for k in range(rep):
    x = np.random.randint(0,8)
    loss = np.abs((q[x]@Et.reshape(3,3)@q[x]) - Eq[x])
    arr.append(loss)
    
    for i in range(9):
        Ett = np.copy(Et)
        Ett[i] += h
        Eu[i] = (np.abs((q[x]@Ett.reshape(3,3)@q[x]) - Eq[x]) - loss)/h
    
    for i in range(9):
        Et[i] -= Eu[i]*e 
        
plt.plot(arr)
plt.show()

<Figure size 432x288 with 1 Axes>

まぁまぁですかね。

np.round(Et, 2)

array([ 0.2 , -0.31,  0.3 , -0.51, -0.05,  0.28, -0.  , -0.12,  0.87])

多少相互作用のところに不満は残りますが、エネルギー値は近いものになったようです。

array([ 1.99, -1.4 ,  0.15, -1.6 ,  0.  ,  0.2 , -0.15, -0.2 ,  0.99])

確認

最後にちょっと確認をします。コピー元のイジング式とのエネルギー値を比較します。

#最後に[ 0.,  1.,  0.,  1.,  2.,  3., -1.,  0.]と比較します。
Ef = np.zeros(9)

for i in range(8):
    Ef[i] = q[i]@Et.reshape(3,3)@q[i]
Ef

array([ 0.        ,  0.87381312, -0.05109759,  0.99254179,  0.1951    ,
        1.36051312, -0.67279759,  0.66244179,  0.        ])

これは、ある程度コピーできていると思います。応用として、イジング式のわからない事象をイジング式に定式化できます。

今後の発展

最適化アルゴリズムはSGDでしたが、これもハイパラなので、より良いのもありそうです。是非色々試してみてください。以上です。