量子コンピュータ
量子アニーリング
qubo

はじめに

ちょっとアドベントカレンダーの間が空いてるようでしたので、頭の体操で問題を作って見ました。せっかくなので量子アニーリングで1+1をして見たいと思います。四則演算はのやり方は下の方に書いてあります。

基本

四則演算を行う際には数字を量子ビット表記します。
求めたい値を$x = q_0+2q_1$とすると、

q_3+q_4 = x

が今回求めたい$q_3+q_4$の足し算です。早速式変形をします。途中式で$q^2 = q$となっているのは、量子ビットが{0,1}を取りうるため、$0^2=0$もしくは$1^2=1$が成り立つため、二次式を一次に直しています。

\begin{eqnarray*}
(x-q_3-q_4)^2 &=& (q_0+2q_1-q_3-q_4)^2\\
&=& q_0^2+4q_1^2+q_3^2+q_4^2+4q_0q_1-2q_0q_3-2q_0q_4-4q_1q_3-4q_1q_4+2q_3q_4\\
&=& q_0+4q_1+q_3+q_4+4q_0q_1-2q_0q_3-2q_0q_4-4q_1q_3-4q_1q_4+2q_3q_4\\
\end{eqnarray*}

この式の最小値0を求めると$q_3+q_4$が成り立つことと同じになるので、これが量子アニーリングで足し算をすることに相当します。

実際の量子アニーリングはビット表記を{0,1}ではなく、物理モデルに適用できるように{-1,1}を使用します。その際の変換は単純で$q=\frac{q'+1}{2}$を代入します。

\begin{eqnarray*}
&=& q_0'q_1'-\frac{1}{2}q_0'q_3'-\frac{1}{2}q_0'q_4'-q_1'q_3'-q_1'q_4'+\frac{1}{2}q_3'q_4'+\frac{1}{2}q_0'+q_1'-\frac{1}{2}q_3'-\frac{1}{2}q_4'+2\\
\end{eqnarray*}

これは、二次計画問題のQUBOと呼ばれる形なので、量子アニーリングでとけます。
ちなみにQUBOについては下記の過去の記事参照。

QUBOを使って日本地図の彩色問題やってみた
https://qiita.com/gyu-don/items/dc83392cb890db6c8d61

sqa3.py
from numpy import *

#run sqa algorithm with params

h = [0.5,1,-0.5,-0.5]
J = [[0,1,-0.5,-0.5],[1,0,-1,-1],[-0.5,-1,0,0.5],[-0.5,-1,0.5,0]]
c = 2
N = 4
m = 10
rep = 10
kT = 0.02
tau = 0.99
Ginit = 5
Gfin = 0.01
qarr = []
Earr = []

# simulated quantum annealing simulator using quantum monte carlo & metropolis

for j in range(rep):
        G = Ginit

        q = []
        for i in range(m):
                q.append(random.choice([-1,1],N))

        while G > Gfin:
                for i in range(N*m):
                        x = random.randint(N)
                        y = random.randint(m)
                        dE = -2*h[x]*q[y][x]*1.0/m
                        for k in range(N):
                                if k != x:
                                        dE += -2*J[x][k]*q[y][x]*q[y][k]*1.0/m
                        dE += -q[y][x]*(q[(m+y-1)%m][x]+q[(y+1)%m][x])*log(tanh(G/kT/m))*1.0/kT
                        if exp(-dE/kT)>random.rand():
                                q[y][x] = -q[y][x]
                G*=tau

        E = 0
        for a in range(N):
                E += h[a]*q[0][a]
                for b in range(a+1,N):
                        E += J[a][b]*q[0][a]*q[0][b]
        qarr.append(q)
        Earr.append(E+c)

print(array([qarr,Earr]))

結果を見ると、

result
[[ list([array([ 1, -1,  1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([ 1, -1,  1, -1]), array([ 1, -1,  1, -1]), array([ 1, -1,  1, -1]), array([ 1, -1,  1, -1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1])])
  list([array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1])])
  list([array([-1, -1, -1, -1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1])])
  list([array([-1, -1, -1, -1]), array([ 1, -1,  1, -1]), array([ 1, -1,  1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1])])
  list([array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([-1, -1, -1, -1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([-1, -1, -1, -1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1]), array([ 1, -1, -1,  1])])]
 [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]]

量子アニーリングのサンプリングで複数回試行したので複数の基底状態が出てきてしまいました。たとえば、$q_0=1,q_1=0,q_3=1,q_4=0$とか、$q_0=0,q_1=0,q_3=0,q_4=0$とか、$q_0=0,q_1=1,q_3=1,q_4=1$とか、$q_0=1,q_1=0,q_3=0,q_4=1$となりました。

これは、

q_3 + q_4 = 1 + 0 = q_0+2q_1 = 1+0 = 1\\
q_3 + q_4 = 0 + 0 = q_0+2q_1 = 0+0 = 0\\
q_3 + q_4 = 1 + 1 = q_0+2q_1 = 0+2 = 2\\
q_3 + q_4 = 0 + 1 = q_0+2q_1 = 1+0 = 1

となりました。式としてはあっていました。さらに基底状態がたくさんある計算式なので、あっている答えを同時に全部計算してしまったようです。計算したい数字は$q_3 = 1, q_4 = 1$でした。
上記のQUBOだと$q_3$と$q_4$に初期値をセットしてないので、$q_3+q_4$を満たす全ての解が出てしまいます。

ですので、$q_3$と$q_4$に初期値をセットします。ここでは、$q_3 = 1, q_4 = 1$をみたすため、
$q_3$と$q_4$に磁場項hを$\sigma_z$方向に1になるように作用させることで実現します。
適当に$q_3$と$q_4$の磁場項hに対応する部分に-5を作用させると、

sqa4.py
from numpy import *

#run sqa algorithm with params

h = [0.5,1,-5.5,-5.5]
J = [[0,1,-0.5,-0.5],[1,0,-1,-1],[-0.5,-1,0,0.5],[-0.5,-1,0.5,0]]
c = 2
N = 4
m = 10
rep = 10
kT = 0.02
tau = 0.99
Ginit = 5
Gfin = 0.01
qarr = []
Earr = []

# simulated quantum annealing simulator using quantum monte carlo & metropolis

for j in range(rep):
        G = Ginit

        q = []
        for i in range(m):
                q.append(random.choice([-1,1],N))

        while G > Gfin:
                for i in range(N*m):
                        x = random.randint(N)
                        y = random.randint(m)
                        dE = -2*h[x]*q[y][x]*1.0/m
                        for k in range(N):
                                if k != x:
                                        dE += -2*J[x][k]*q[y][x]*q[y][k]*1.0/m
                        dE += -q[y][x]*(q[(m+y-1)%m][x]+q[(y+1)%m][x])*log(tanh(G/kT/m))*1.0/kT
                        if exp(-dE/kT)>random.rand():
                                q[y][x] = -q[y][x]
                G*=tau

        E = 0
        for a in range(N):
                E += h[a]*q[0][a]
                for b in range(a+1,N):
                        E += J[a][b]*q[0][a]*q[0][b]
        qarr.append(q)
        Earr.append(E+c)

print(array([qarr,Earr]))

となり、結果は、

result
[[ list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])
  list([array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1]), array([-1,  1,  1,  1])])]
 [-10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0]]

のようになり、最終のコスト関数は作用させた-5*2の分だけ0-10となっていますが、

$q_0=0,q_1=1,q_3=1,q_4=1$となり、

実際に

1+1 = 2

が量子アニーリングだけで計算できました!

結局

量子アニーリングで計算しても1+1は2です!

四則演算について

気になる人もいると思うので、ついでに四則演算をサポートします。

2-1の引き算

2-1 = 1

を求めるためには、

(q_3+2q_4) - q_5 = x

となる$x$を求めます。$x = q_0$とすると、

\begin{eqnarray*}
q_3+2q_4-q_5 &=& x\\
x-q_3-2q_4+q_5 &=& 0\\
(x-q_3-2q_4+q_5)^2 &=& 0\\
(q_0-q_3-2q_4+q_5)^2 &=& 0
\end{eqnarray*}

こうなりましたのでは、下記の式の最小値をみたす量子ビットの組み合わせを求めれば引き算ができます。

(q_0-q_3-2q_4+q_5)^2

2*1の掛け算

2\times1 = 2

の計算をするためには、

\begin{eqnarray*}
(q_3+2q_4) \times q_5 &=& x\\
q_3q_5+2q_4q_5 &=& x\\
x-q_3q_5-2q_4q_5 &=& 0\\
q_0+2q_1-q_3q_5-2q_4q_5 &=& 0\\
(q_0+2q_1-q_3q_5-2q_4q_5)^2 &=& 0

\end{eqnarray*}

ですので、同様に、下記の式の最小値0をみたすそれぞれの量子ビットを求めればOK。

(q_0+2q_1-q_3q_5-2q_4q_5)^2

2/1の割り算

2\div 1 = 2

こちらは

(q_3+2q_4)\div q_5 = x

を解けばok

\begin{eqnarray*}
q_3+2q_4 &=& xq_5\\
xq_5-q_3-2q_4 &=& 0\\
(q_0+2q_1)q_5-q_3-2q_4 &=& 0\\
q_0q_5+2q_1q_5-q_3-2q_4 &=& 0\\
(q_0q_5+2q_1q_5-q_3-2q_4)^2 &=& 0
\end{eqnarray*}

ですので、

(q_0q_5+2q_1q_5-q_3-2q_4)^2

こちらの式を最小値0にする量子ビットの組み合わせを求めます。

多体問題

上記の式を見て気づいた方は掛け算と割り算には最小値問題に落とし込むときに4つの量子ビットの掛け算となる多体問題が出ることが確認できます。

この多体問題は現在の量子アニーリングマシンでは2体問題までしか直接は計算できないので、解くためにはさらに一工夫必要です。

多体問題の分解には通常グレブナー基底を用いて式展開します。