やること
誰かが言いました。
連立方程式は組合せ最適化問題である [要出典]
言ってないんだよなぁ。
今回はQUBOで連立方程式(文字が3つあるやつ)を解いてみます。
考え方
基本の条件設定「n個の量子ビットからm個を1にする」から考えていきます。
例えば、「3個の量子ビットから2個を1にする」は次の式です。
H = (q0 + q1 + q2 - 2)**2
これはつまり、3個の量子ビットを足した値が2になるということです。それぞれの量子ビットは [0, 1] のどちらかを取ります。合計値が2であればエネルギーは最小で、合計値が2からずれるとエネルギーが高くなります。
ここで、量子ビット名を [x, y, z] に変えて、それぞれに係数 [5, -1, 2] をかけてみます。
H = (5*x - y + 2*z - 7)**2
実はこれで方程式「5x – y + 2z = 7」が完成しました。「5x – y + 2z」の部分が7のときにエネルギーは最小で、7からずれるとエネルギーが高くなります。
このような方程式を複数設定します。エネルギーが最小の量子ビットの状態を求めるということは、すなわち連立方程式を解くということです。数学の問題ではピッタリの解が存在するはずなので、アニーリングでは理論上の最小エネルギーの解が得られるはずです。
その他の設定可能な条件式も気になる方は過去のまとめ記事をご参照ください↓(宣伝)
解が 0, 1 だけの場合
次の3元連立1次方程式を解きます。
(1)
解が0または1であることが分かっていることを前提とします。解が0と1だけなので [x, y, z] をそのまま量子ビットで表せます。
from pyqubo import Binary
from dwave_qbsolv import QBSolv
#量子ビットを用意する
x = Binary('x')
y = Binary('y')
z = Binary('z')
#連立方程式の設定
H = 0
H += ( 5*x - y +2*z - 7)**2
H += (-3*x +4*y + z + 2)**2
H += ( x -2*y -4*z + 3)**2
#コンパイル
model = H.compile()
qubo, offset = model.to_qubo()
print(f'offset\n{offset}')
#分割サンプリング
response = QBSolv().sample_qubo(qubo, num_repeats=500)
#レスポンスの確認
print('response')
for (sample, energy, occurrence) in response.data():
print(f'Sample:{list(sample.values())} Energy:{energy} Occurrence:{occurrence}')
offset
62.0
response
Sample:[1, 0, 1] Energy:-62.0 Occurrence:501
解が [1, 0, 1] と出ました。オフセットと解のエネルギーに着目すると、マイナスオフセットが理論上の最小エネルギーで、それに一致するエネルギーの解が得られています。この解は3つの方程式を完全に満たしていることを意味しています。
解が 0, 1, 2, 3 だけの場合
次の3元連立1次方程式を解きます。
(2)
解が0~3の整数であることが分かっていることを前提とします。この場合、[x, y, z] をそのまま量子ビットで表すことができません。そこで [x0, x1] という量子ビットを用意して、「x = 2*x0 + x1」のように2進数表現(2bit)します。x0が2^1の桁、x1が2^0の桁です。
from pyqubo import Binary
from dwave_qbsolv import QBSolv
#量子ビットを用意する
x0 = Binary('x0')
x1 = Binary('x1')
y0 = Binary('y0')
y1 = Binary('y1')
z0 = Binary('z0')
z1 = Binary('z1')
#x,y,zを2進数(2bit)で表す
x = 2*x0 + x1
y = 2*y0 + y1
z = 2*z0 + z1
#連立方程式の設定
H = 0
H += ( x + y + z - 6)**2
H += (2*x +3*y -2*z - 11)**2
H += (3*x - y + z - 4)**2
#コンパイル
model = H.compile()
qubo, offset = model.to_qubo()
print(f'offset\n{offset}')
#分割サンプリング
response = QBSolv().sample_qubo(qubo, num_repeats=500)
#レスポンスの確認
print('response')
for (sample, energy, occurrence) in response.data():
print(f'Sample:{list(sample.values())} Energy:{energy} Occurrence:{occurrence}')
#2進数から戻して確認
x0, x1 = list(sample.values())[:2]
y0, y1 = list(sample.values())[2:4]
z0, z1 = list(sample.values())[4:]
print(f'x = {2*x0 + x1}')
print(f'y = {2*y0 + y1}')
print(f'z = {2*z0 + z1}')
offset
173.0
response
Sample:[1, 0, 1, 1, 0, 1] Energy:-173.0 Occurrence:501
x = 2
y = 3
z = 1
正しい解が得られました!
解が 0~255 の整数の場合
最後に、次の3元連立1次方程式にチャレンジします。
(3)
解が0~255の整数であることが分かっていることを前提とします。お察しの通り、[x0, x1, x2, x3, x4, x5, x6, x7] という量子ビットを用意して、「x = 128*x0 + 64*x1 + 32*x2 + 16*x3 + 8*x4 + 4*x5 + 2*x6 + x7」という2進数表現(8bit)をします。3元の連立方程式に24量子ビットを使います。
from pyqubo import Binary
from dwave_qbsolv import QBSolv
#量子ビットを用意する
x0 = Binary('x0')
x1 = Binary('x1')
x2 = Binary('x2')
x3 = Binary('x3')
x4 = Binary('x4')
x5 = Binary('x5')
x6 = Binary('x6')
x7 = Binary('x7')
y0 = Binary('y0')
y1 = Binary('y1')
y2 = Binary('y2')
y3 = Binary('y3')
y4 = Binary('y4')
y5 = Binary('y5')
y6 = Binary('y6')
y7 = Binary('y7')
z0 = Binary('z0')
z1 = Binary('z1')
z2 = Binary('z2')
z3 = Binary('z3')
z4 = Binary('z4')
z5 = Binary('z5')
z6 = Binary('z6')
z7 = Binary('z7')
#x,y,zを2進数(8bit)で表す
x = 128*x0 + 64*x1 + 32*x2 + 16*x3 + 8*x4 + 4*x5 + 2*x6 + x7
y = 128*y0 + 64*y1 + 32*y2 + 16*y3 + 8*y4 + 4*y5 + 2*y6 + y7
z = 128*z0 + 64*z1 + 32*z2 + 16*z3 + 8*z4 + 4*z5 + 2*z6 + z7
#連立方程式の設定
H = 0
H += (10*x +14*y +4*z - 5120)**2
H += ( 9*x +12*y +2*z - 4230)**2
H += ( 7*x + 5*y +2*z - 2360)**2
#コンパイル
model = H.compile()
qubo, offset = model.to_qubo()
print(f'offset\n{offset}')
#分割サンプリング
response = QBSolv().sample_qubo(qubo, num_repeats=1000)
#レスポンスの確認
print('response')
for (sample, energy, occurrence) in response.data():
print(f'Sample:{list(sample.values())} Energy:{energy} Occurrence:{occurrence}')
#2進数から戻して確認
x0,x1,x2,x3,x4,x5,x6,x7 = list(sample.values())[:8]
y0,y1,y2,y3,y4,y5,y6,y7 = list(sample.values())[8:16]
z0,z1,z2,z3,z4,z5,z6,z7 = list(sample.values())[16:]
print(f'x = {128*x0+64*x1+32*x2+16*x3+8*x4+4*x5+2*x6+x7}')
print(f'y = {128*y0+64*y1+32*y2+16*y3+8*y4+4*y5+2*y6+y7}')
print(f'z = {128*z0+64*z1+32*z2+16*z3+8*z4+4*z5+2*z6+z7}')
offset
49676900.0
response
Sample:[1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0] Energy:-49676900.0 Occurrence:9
x = 130
y = 230
z = 150
Sample:[1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0] Energy:-49676871.0 Occurrence:9
x = 130
y = 229
z = 154
Sample:[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0] Energy:-49676814.0 Occurrence:9
x = 127
y = 232
z = 152
サンプリング数を1000に増やしていますが、何回も実行してたまたま正解が出たときの結果を載せています。一応ちゃんと解けましたね!
さいごに
最後のは正解が得られる確率が低かったです。これはアニーリングアルゴリズム(アニーリングマシン)の性能によるところでしょうか。神様が理想的なアニーリングマシンを用意してくれたら、ただの1回のアニーリングで最適解が得られるはずなのですが。