量子アニーリングは巷のニュースでは最適化が得意と書いてありますが、実際にどのようなステップで最適化計算が進むのか、イメージとして確認をしたいと思います。その中で、最適化とサンプリングと呼ばれる2つの手法について確認をしたいと思います。
ステップ
基本のステップは
１、定式化
２、サンプリング
３、最適化
となります。早速みてみましょう。
１、定式化
まず量子アニーリングマシンで必ず行わないといけないことがあります。それが定式化です。イジングモデルと呼ばれるモデルを使っているので、それに準じた答えがでます。数式は量子ビット単体にかかるバイアスと、量子ビット間にかかるウェイトがあります。（文言は色々あります。）
E=−∑hiqi−∑JijqiqjE = -\sum h_i q_i -\sum J_{ij}q_iq_jE=−∑hiqi−∑Jijqiqj
のように、バイアスhiとウェイトJijを決めることで定式化できます。qは量子ビットの値を表し、自分では設定しません。答えとしてqの値が出ます。
２、サンプリング
次に量子アニーリングマシンが答えを出すメカニズムがあります。それは、理想的には、
p(x)=1Zexp(−E)p(x) = \frac{1}{Z} exp(-E)p(x)=Z1exp(−E)
のように、ある答えが出る確率は、イジングモデルのエネルギーに依存します。Zは規格化定数というもので、計算が大変ですが、ここでは、サンプリングを使ってこの確率分布を近似できます。
量子アニーリングマシンでは、この理想的に出る分布を仮定して、出る答えは確率的に決まります。必ずしも最小値が出るとは限らず、問題の設定によって変わります。そのため、何度も計算をして、確率分布を取る必要があります。
３、最適化
最後に、出た答えの中で一番小さいものを選ぶとそれが最適化です。最適化に至るには、自分で設定した問題の確率分布に応じて最適解が出る確率が変わることがわかります。
Advantageを使ってやってみる
適当に問題を決めてやってみます。
M=0.1
qubo = np.array([[0,1,1,-1],[0,0,1,-1],[0,0,0,1],[0,0,0,0]])*M
print(qubo)
#show as network
G = nx.from_numpy_matrix(qubo)
nx.draw_networkx(G)
plt.show()
こんな感じの接続になります。Mは適宜変えると色々変わります。ここでは、まず0.1にしてみます。
#to Advantage
bqm = dimod.BQM.from_numpy_matrix(qubo)
sampler = EmbeddingComposite(DWaveSampler(endpoint='https://cloud.dwavesys.com/sapi', token=token, solver='Advantage_system1.1'))
response = sampler.sample(bqm, num_reads=1000)
print(response)
1000回で計算してみます。
0 1 2 3 energy num_oc. chain_.
0 1 0 0 1 -0.1 278 0.0
1 0 1 0 1 -0.1 287 0.0
2 1 1 0 1 -0.1 434 0.0
3 1 0 0 0 0.0 1 0.0
['BINARY', 4 rows, 1000 samples, 4 variables]
結果はこうなりました。横方向が1回の計算結果です。1001から1101までの結果はトータルのエネルギー値が-0.1になっています。出現確率はばらつきがありますが、278から434になってます。エネルギー値が0で少し高めのものでも出ていますが、こちらは1回出現しています。
M=1にしてみると、
0 1 2 3 energy num_oc. chain_.
0 1 1 0 1 -1.0 470 0.0
1 0 1 0 1 -1.0 191 0.0
2 1 0 0 1 -1.0 339 0.0
['BINARY', 3 rows, 1000 samples, 4 variables]
M=5
0 1 2 3 energy num_oc. chain_.
0 1 0 0 1 -5.0 371 0.0
1 1 1 0 1 -5.0 389 0.0
2 0 1 0 1 -5.0 240 0.0
['BINARY', 3 rows, 1000 samples, 4 variables]
M=10
0 1 2 3 energy num_oc. chain_.
0 1 0 0 1 -10.0 244 0.0
1 1 1 0 1 -10.0 361 0.0
2 0 1 0 1 -10.0 394 0.0
3 0 0 0 0 0.0 1 0.0
['BINARY', 4 rows, 1000 samples, 4 variables]
少し行列を変えてみます。
M=1
qubo = np.array([[0,1,1,-2],[0,0,1,0],[0,0,0,-1],[0,0,0,0]])*M
0 1 2 3 energy num_oc. chain_.
0 1 0 0 1 -2.0 459 0.0
1 1 0 1 1 -2.0 539 0.0
2 0 0 1 0 0.0 1 0.0
3 1 1 1 1 0.0 1 0.0
['BINARY', 4 rows, 1000 samples, 4 variables]
こんな感じになりました。
M=0.02
qubo = np.array([[0,1,5,-2],[0,0,1,0],[0,0,0,-1],[0,0,0,0]])*M
0 1 2 3 energy num_oc. chain_.
0 1 0 0 1 -0.04 956 0.0
1 1 1 0 1 -0.02 10 0.0
2 0 0 1 1 -0.02 30 0.0
3 0 1 0 1 0.0 2 0.0
4 0 1 1 1 0.0 2 0.0
['BINARY', 5 rows, 1000 samples, 4 variables]
面白いですね。行列ないの最小値と最大値の差が大きい場合には影響が大きそうです。
もちょいやってみたところ、きちんとばらつきが出ました。
M=0.02
qubo = np.array([[0,1,10,-2],[0,0,1,0],[0,0,0,-1],[0,0,0,0]])*M
0 1 2 3 energy num\_oc. chain\_.
0 1 0 0 1 -0.04 706 0.0
1 1 1 0 1 -0.02 63 0.0
2 0 0 1 1 -0.02 129 0.0
3 0 0 1 0 0.0 23 0.0
4 0 1 0 0 0.0 11 0.0
5 1 0 0 0 0.0 20 0.0
6 0 1 0 1 0.0 10 0.0
7 0 0 0 0 0.0 12 0.0
8 0 1 1 1 0.0 6 0.0
9 0 0 0 1 0.0 18 0.0
10 0 1 1 0 0.02 2 0.0
['BINARY', 11 rows, 1000 samples, 4 variables]
M=1にしてみると、
0 1 2 3 energy num_oc. chain_.
0 1 0 0 1 -2.0 482 0.0
1 1 1 0 1 -1.0 44 0.0
2 0 0 1 1 -1.0 388 0.0
3 1 0 0 0 0.0 8 0.0
4 0 0 0 1 0.0 6 0.0
5 0 1 1 1 0.0 25 0.0
6 0 0 1 0 0.0 29 0.0
7 0 0 0 0 0.0 4 0.0
8 0 1 0 0 0.0 9 0.0
9 0 1 0 1 0.0 5 0.0
['BINARY', 10 rows, 1000 samples, 4 variables]
結構面白いです。みんなで色々やってみましょう。