概要
MDR社が開発した、イジングモデルを扱った開発ツールWildqat(http://mdrft.com/wildqat/) のチュートリアルを行いました。
Wildqatはpythonをベースとしてイジングモデルという物理モデルを扱いながら、物理の知識がなくても計算を行うことのできる開発ツールです。
計算はローカルのマシンでシミュレーションを行い、計算を行うことができます。
(http://mdrft.com/wildqat/ から引用)
本記事ではそこで出た疑問やちょっとしたメモなどを書いておきます。
チュートリアルは以下から飛べます。
https://github.com/mdrft/Wildqat/tree/master/examples_ja
書いた人の知識
修士まで物性を用いた量子関係の研究(重ね合わせの保持時間)
プログラミングは素人レベル(半年だけSEしてました)
量子コンピュータの勉強も最近はじめたので、間違いなどあればコメントもしくはtwitterで指摘していただければ泣いて喜びます
チュートリアル1
QUBO
いきなりQUBOという単語が出てきます。
QUBOはQuadratic Unconstrained Binary Optimizationの略で、日本語では制約なし二次形式二値変数最適化というみたいです[1]。詳しい説明は[2]が参考になると思いますが、とりあえず使えればいいという観点で説明します。
Qiita内に素晴らしい説明があった[3]のでそれを引用します。
QUBOとは
H=∑Q_{ij}x_ix_j
の、Hを最小にするように $x_i∈0,1$ を発見する最適化問題でした。
QUBOを作るときに考えないといけないのは、 $x_i$ とは何を表現する量か、そのために$Q$をどのようにしたらいいか、です。
まず$H$ですが、$H$はコスト関数等と呼ばれます。
$H$を作る流れは以下のとおりです。
1. 最適化したい問題を$x_i \in 0, 1\ \ (i = 1, 2, \cdots)$を用いるように解釈し直す
2. 最適になれば最小になる式を探す←この式がコスト関数
QUBOは1.の"解釈し直す"部分に相当すると思っています。。
また、この手順で作ったコスト関数は$q_{ii}x_i$、または$q_{ij}x_ix_j$の形($q_{ii}, q_{ij}$は定数)で書け、この$q_{ii}, q_{ij}$を行列の形にしたものをQUBO行列と呼びます。
チュートリアル1ではQUBO行列は既に与えられているので、後はチュートリアルに従うだけで終了です。
# チュートリアル2
チュートリアル１で書いた、QUBOに解釈し直した後に、コスト関数を作ることについてメモを残しておきます。
といっても、解きたい等式を変形して
```math
f(x)-C = 0 (C:定数)
という形まで持っていって、これが極力成り立つような$x$が知りたいので、そのためには、
(f(x)-C)^2
の値を極力小さくしたいのです。ですので、この式を$E(x)$(:コスト関数)とおいてエネルギーっぽく解釈しているわけですね。
これがチュートリアル2における$x-2$に該当していることは言うまでもないかもしれません。
チュートリアル2中で
一方、xは量子ビットを使って
x = q0 + 2q1
という二進数表記ができますので、
という文章があり、理解できませんでしたが、以下のように捉えれば良いと思います。
取り敢えず解となる$x$が高々2bitで表せられる範囲だとする、つまり$0 \leq x \leq 3$であるとする。
すると、$x$は2量子ビット、{$q_0, q_1$}を用いて
x = q_0 + 2q_1
として表せられる。
これをコスト関数$E(x)$に代入すると
\begin{align*}
E(x) &= (x-2)^2 \\
&= (q_0+2q_1-2)^2 \\
&= q_0^2+4q_1^2+4+4q_0q_1-8q_1-4q_0
\end{align*}
で、この$q_0, q_1$はbitであり、バイナリ値$0, 1$をとるので$q_0, q_1$がそれぞれ0でも1でも２乗すればもとの数(つまり0か1)になることを考慮すれば
\begin{align*}
E(x) &= -3q_0+4q_0q_1-4q_1+4
\end{align*}
となりますね。
ここで今知りたいのは$E(x)$が最小になるような$x$の値で、$E(x)$自体の値は興味ありません。ですので$E(x)$に定数をいくらか足しても引いても掛けても問題ありません。ですので定数項は無視しちゃいます。
つまり
\begin{align*}
E(x) &= -3q_0+4q_0q_1-4q_1
\end{align*}
とします。
次にこれを行列表記したい、つまり
\begin{align*}
E(x) &= \left(
\begin{array}{c}
q_0 & q_1
\end{array}
\right)
\left(
\begin{array}{cc}
a & b \\
c & d
\end{array}
\right)
\left(
\begin{array}{c}
q_0 \\
q_1
\end{array}
\right)\\
&=\left(
\begin{array}{c}
aq_0+cq_1 & bq_0 +dq_1
\end{array}
\right)
\left(
\begin{array}{c}
q_0 \\
q_1
\end{array}
\right)\\
&=aq_0^2+cq_0q_1+bq_0q_1+dq_1^2\\
&=aq_0+cq_0q_1+bq_0q_1+dq_1 \ (\because q_0^2 = q_0, q_1^2 = q_1)
\end{align*}
としたいわけです。
そのためには、上記の$a, b, c, d$を要素に持つ2×2の行列がわかればよいのです。
ここで$q_0q_1$の項が2つ出てきますが、簡単のため項を1つだけにしたいので、$c=0$としましょう[4]。
すると
E(x) = aq_0+bq_0q_1+dq_1
となりますので、行列は
\left(
\begin{array}{cc}
a & b \\
c & d
\end{array}
\right) = \left(
\begin{array}{cc}
-3 & 4 \\
0 & -4
\end{array}
\right)
となります。
結果が[0, 1]となった、つまり$q_0 = 0, q_1 = 1$(チュートリアルは誤植です)となったので、求めたい$x$は
x = 0 + 2\times1 = 2
となりました。つまり$1+1=2$が計算できた、というわけですね。
チュートリアル3
コスト関数の展開
\begin{align*}
E &= \left\{ \sum_{i=1}^{N} n_i(2q_i-1)\right\}^2\\
&= \left\{ \sum_{i=1}^{N} 2n_iq_i-\sum_{i=1}^{N} n_i\right\}\left\{ \sum_{i=1}^{N} 2n_iq_i-\sum_{i=1}^{N} n_i\right\}\\
&= \left\{ \sum_{i=1}^{N} 2n_iq_i\right\}^2 -2\left\{ \sum_{i=1}^{N} 2n_iq_i\right\}\left\{ \sum_{i=1}^{N}n_i\right\} + \left\{ \sum_{i=1}^{N} 2^n_i\right\}^2\\
\end{align*}
紛らわしくなりますが$\sum_{i=1}^{N}n_i = n_{sum}$自体はN個の自然数を定めた時点で決まる定数項ですので、最後の項は無視できます。さらにチュートリアル2で書いたように定数で掛けても問題ないので、4で割って
\begin{align*}
E &= \left\{ \sum_{i=1}^{N} n_iq_i\right\}^2 -\left\{ \sum_{i=1}^{N} n_iq_i\right\}\left\{ \sum_{i=1}^{N}n_i\right\} \\
\end{align*}
一項目について、
\begin{align*}
\left\{ \sum_{i=1}^{N} n_iq_i\right\}^2&=\left(n_1q_1 + n_2q_2 + n_3q_3 + \cdots + n_Nq_N\right)\left(n_1q_1 + n_2q_2 + n_3q_3 + \cdots + n_Nq_N\right)\\
&= \left(n_1^2q_1^2+n_2^2q_2^2+n_3^2q_3^2 + \cdots n_N^2q_N^2\right) + 2\left(n_1n_2q_1q_2+n_1n_3q_1q_3+\cdots+n_1n_Nq_1q_N+n_2n_3q_2q_3+n_2n_4q_2q_4+\cdots +n_2n_Nq_2q_N + \cdots \right)\\
&= \left(n_1^2q_1+n_2^2q_2+n_3^2q_3 + \cdots n_N^2q_N\right) + 2\left(n_1n_2q_1q_2+n_1n_3q_1q_3+\cdots+n_1n_Nq_1q_N+n_2n_3q_2q_3+n_2n_4q_2q_4+\cdots +n_2n_Nq_2q_N + \cdots \right)\\
&= \sum_{i=1}^{N}n_i^2q_i + 2\sum_{i<j}n_in_jq_iq_j
\end{align*}
二項目については$\sum_{i=1}^{N}n_i = n_{sum}$とすることに注意すれば
\begin{align*}
E &= \sum_{i=1}^{N} (n_i^2-n_{sum}n_i)q_i +2\sum_{i<j}n_in_jq_iq_j \\
\end{align*}
ここからQUBO行列を求めるわけですが、チュートリアル2から急にN×N行列になって困惑するかもしれません。
しかし、対角項は$q_1, q_2, q_3, \cdots q_N$の係数、要素(QUBO)$_{12}$は$q_1q_2$の係数、要素(QUBO)$_{21}$は0、というように考えていけば大丈夫です。((QUBO)$_{nm}$はQUBO行列のn行m列の要素という意味で使っています)。するとチュートリアルにあるような行列が得られます。
結果は
>print(answer)
[1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0]
のように、最初に定義したnumbersの各indexにある数値がどちらのグループになるかで0か1かが割り振られています。
もう一度実行すると
[1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0]
となったりして、結果が変化しますね。面白いです。
チュートリアル4
- 例題を元にしながら考えたほうがわかりやすいと思うので、そうします。
コスト関数に関して、制約条件$Sx=b$は例題では
\left(
\begin{array}{ccc}
3 & 2 & 1\\
5 & 2 & 3
\end{array}
\right)
\left(
\begin{array}{c}
x_1\\
x_2\\
x_3
\end{array}
\right)
= \left(
\begin{array}{c}
3\\
5
\end{array}
\right)
\\
\Leftrightarrow
\left(
\begin{array}{c}
3-(3x_1+2x_2+x_3)\\
5-(5x_1+2x_2+3x_3)
\end{array}
\right)
=0
となります。つまり、左辺の各要素について0になれば良いので、
\begin{align}
(3-(3x_1+2x_2+x_3))^2 + (5-(5x_1+2x_2+3x_3))^2 \rightarrow 0
\end{align}
となるようにすればいい（=最小にする)、ということになり、コスト関数の第一項が出てきます。
コスト関数の第二項は最小化したい項、そのものですね。
ここからQUBO行列を求めるわけですが、まずコスト関数第一項目のQUBO行列を求めましょう。
b_1 = 3 \\
A_1 = 3x_1 + 2x_2 + x_3
としましょう。するとコスト関数一項目の$j=1$の項は
\begin{align}
(b_1 + A_1)^2 &= b_1^2 + 2b_1A_1+A_1^2\\
&\rightarrow 2b_1A_1 + A_1^2 \ \ (\because \ b_1 \in \textrm{ const.})
\end{align}
となります。
QUBO行列が
\left(
\begin{array}{ccc}
x_1 & x_1x_2 & x_1x_3\\
0 & x_2 & x_2x_3 \\
0 & 0 & x_3
\end{array}
\right)
の順に係数を入れれば良いので、$-2b_1A_1$については、numpyのdiagを使えば[5]
-2 * 3 * np.diag([3, 2, 1])
で表現できます。
$A_1^2$についてはwildqatの方で便利なsqrを用意してくれています。
これはwildqat.sqr([3, 2, 1])とすれば
\begin{align}
A_1^2 &= (3x_1 + 2x_2 + x_3)^2\\
&= 9x_1^2 + 4x_2^2 + x_3^2 + 12x_1x_2 + 4x_2x_3 + 6x_3x_1\\
&= 9x_1 + 4x_2 + x_3 +12x_1x_2 + 4x_2x_3 + 6x_3x_1
\end{align}
についてQUBO行列に変換して
\left(
\begin{array}{ccc}
9 & 12 & 6\\
0 & 4 & 4 \\
0 & 0 & 1
\end{array}
\right)
としてくれます。
以上を用いればQUBO行列の作成が
A = [[3,2,1],[5,2,3]]
b = [3,5]
a.qubo = np.zeros((3,3))
for i in range(len(b)):
a.qubo += -2*b[i]*np.diag(A[i]) + wq.sqr(A[i])
となるのも頷けるでしょう。
ここまでくればコスト関数の二項目はすぐに
matrix2 = np.diag([1,2,1])
に係数$B$をかければ良いことがわかります。
そしてこれをマニュアル通りに解けばよいのですがたまに
a.sa()
> [1, 0, 0]
となります。
a.plot()でコスト関数を見てみると、およそ-35で収束していることがわかりました。
また100回試してみると、うち16回で[1, 0, 0]になることがわかりました。
おそらくコスト関数の谷が[0, 1, 1]のものと[1, 0, 0]の２つあって初期条件によって[1, 0, 0]に収束することがあるんだと思います。
参考文献
[2] https://qiita.com/gyu-don/items/d8536b719f343cc365d6
[3] https://qiita.com/gyu-don/items/0c59423b12256821c38c
[4] https://qard.is.tohoku.ac.jp/T-Wave/?p=454
[5] http://d.hatena.ne.jp/nishiohirokazu/20111109/1320801173