AI要約
量子アニーリング（QUBO）を用いてパズル「スリザーリンク」を解きます。最大の難関である「一本のループにする」制約に対し、直線接続へのボーナス付与という独自の手法でアプローチ。離れ小島の発生を抑え、アニーリングで見事正解を導き出せるか検証します。
やること
久しぶりに量子アニーリングでパズルを解きます。
以前、「スリザーリンクのような線で囲む系はQUBOにできない」と言いました。（「スリザーリンク」で記事検索すると出てきます）
正確に言うと、線を繋ぐのは簡単だが「１本の線で囲む」というルールがQUBO化できないのです。
まあでも、できるところまでやってみましょうよ、ということで今回はスリザーリンクをQUBOで解いてみます。
スリザーリンクとは
いつもお世話になっている Puzzle Team さんから問題を拝借。いつも感謝です。
条件を満たすように格子点に辺を書き込んでいくパズルです。こちらが例題。
左が正しい答え、右が惜しい答えです。
数字は４近傍の辺の数を指定していて、辺は途切れることなく１本のループになっている必要があります。
なんだ基本の「ｎ個の量子ビットからｍ個を１にする」で手を繋いでいけば完成じゃないか、と思うのですが、それだけでは離れ小島の解もできてしまうのです。
これを防ぐ式を作るのが難しいんですね。「自分はあの量子ビットと連続的に繋がっているか？」の関係を記述する方法が思いつきません。（少なくともQUBOの範囲では）
おさらい
今回はこちらの３パターンの条件設定を使います。
基本「ｎ個の量子ビットからｍ個を１にする」
A3「2個の量子ビットが同時に1になったら報酬を与える」
D2 解が0と1だけの方程式（右辺が複数）
定式化
簡単のためにこちらの最小問題を試します。
このような解が想定されます。どちらも手の本数は正しい。
まず、辺に量子ビットを割り当てます。
格子点にも補助量子ビットを割り当てます。
制約条件１は、「各格子点から出る手の本数が０本または２本」です（D2の式）。これで線が途切れたり分岐したりすることはなくなります。
制約条件２は、「数字セルの４辺は数字の数だけ１にする」です。これは基本の式ですね。
これだけだと離れ小島の解も同率１位で出てきてしまいます。
そこでコストとして「各格子点とも直線接続ならプチボーナス」を与えます（A3の式）。２つの解を見比べてみると、一つの島にしたほうが直線接続が多いことがわかります。果たしてこれでいけるでしょうか？？
コードはバカみたいにベタ書きしているところがありますが、ご勘弁を。Pythonで格子の辺にIDを割り振る方法がないのが悪いんですよ（←お前が作れ）。
import itertools
import numpy as np
import matplotlib.pyplot as plt
from tytan import *
#問題
size = [1, 3]
rule = {(0, 1): 2}
#辺の数
num1 = 2 * (size[0] * size[1]) + size[0] + size[1]
print(f'num1 = {num1}')
#点の数
num2 = (size[0] + 1) * (size[1] + 1)
print(f'num2 = {num2}')
#量子ビット
q = symbols_list(num1, 'q{}')
#補助量子ビット
s = symbols_list(num2, 's{}')
#セルに量子ビット（辺）を対応付け
edge = {
(0, 0): [q[0], q[3], q[4], q[7]],
(0, 1): [q[1], q[4], q[5], q[8]],
(0, 2): [q[2], q[5], q[6], q[9]],
}
#各格子点から出る手の本数が０本または２本
Hconst1 = 0
Hconst1 += (q[0] + q[3] - 2*s[0])**2
Hconst1 += (q[0] + q[1] + q[4] - 2*s[1])**2
Hconst1 += (q[1] + q[2] + q[5] - 2*s[2])**2
Hconst1 += (q[2] + q[6] - 2*s[3])**2
Hconst1 += (q[3] + q[7] - 2*s[4])**2
Hconst1 += (q[4] + q[7] + q[8] - 2*s[5])**2
Hconst1 += (q[5] + q[8] + q[9] - 2*s[6])**2
Hconst1 += (q[6] + q[9] - 2*s[7])**2
#数字セルの４辺は数字の数だけ１にする
Hconst2 = 0
for i, (key, value) in enumerate(rule.items()):
Hconst2 += (np.sum(edge[key]) - value)**2
#直線接続ならプチボーナス
Hcost = 0
Hcost -= q[0]*q[1] + q[1]*q[2] + q[7]*q[8] + q[8]*q[9]
#合体
H = Hconst1 + Hconst2 + 0.01*Hcost
#コンパイル
qubo, offset = Compile(H).get_qubo()
print(f'offset = {offset}')
#サンプラー選択
solver = sampler.ArminSampler()
#サンプリング
result = solver.run(qubo, shots=500, T_num=1000)
#上位5件
for r in result[:5]:
print(f'Energy {r[1]}, Occurrence {r[2]}')
#辺の配列（上書きしていることに注意）
q, subs = Auto_array(r[0]).get_ndarray('q{}')
#再度、セルに量子ビット（辺）を対応付け
edge = {
(0, 0): [q[0], q[3], q[4], q[7]],
(0, 1): [q[1], q[4], q[5], q[8]],
(0, 2): [q[2], q[5], q[6], q[9]],
}
#表示
plt.figure(figsize=(8, 8))
ax = plt.axes()
#辺
for _, (key, value) in enumerate(edge.items()):
tmp = edge[key]
i = key[1]
j = key[0]
if tmp[0] == 1:
plt.plot([i, i+1], [j, j], '-k', lw=5)
if tmp[1] == 1:
plt.plot([i, i], [j, j+1], '-k', lw=5)
if tmp[2] == 1:
plt.plot([i+1, i+1], [j, j+1], '-k', lw=5)
if tmp[3] == 1:
plt.plot([i, i+1], [j+1, j+1], '-k', lw=5)
#格子点
for i in range(size[0] + 1):
for j in range(size[1] + 1):
plt.plot(j, i, 'ok', ms=10)
#数字
for _, (key, value) in enumerate(rule.items()):
i = key[0]
j = key[1]
plt.text(j+0.5, i+0.55, f'{value}', ha='center', va='center', weight='bold', size=28)
ax.set_aspect('equal')
plt.xlim(-1, size[1] + 1)
plt.ylim(size[0] + 1, -1)
plt.show()num1 = 10
num2 = 8
offset = 4.0
MODE: GPU
DEVICE: cuda:0
Energy -4.040000915527344, Occurrence 190
Energy -4.0, Occurrence 9
Energy -3.020000457763672, Occurrence 7
（以下略）
上位３つの解を表示しました。正解である第１位が-4.04、離れ小島の第２位は-4.0としっかり差がついています。ボーナスが４ヶ所なので。これはいまくいきそうですね！
例題１
では、冒頭の例題を解いてみましょう。
誰も見ないと思いますが、量子ビット番号を確認したい方はこれを参照してください。
問題サイズに合わせてコードを書き換えます。
#問題
size = [5, 5]
rule = {
(0, 0): 3,
(1, 4): 3,
(2, 1): 3,
(2, 2): 3,
(2, 3): 1,
(2, 4): 3,
(3, 1): 0,
(3, 4): 3,
(4, 1): 3,
}#セルに量子ビット（辺）を対応付け
edge = {
(0, 0): [q[0], q[5], q[6], q[11]],
(0, 1): [q[1], q[6], q[7], q[12]],
(0, 2): [q[2], q[7], q[8], q[13]],
(0, 3): [q[3], q[8], q[9], q[14]],
(0, 4): [q[4], q[9], q[10], q[15]],
(1, 0): [q[11], q[16], q[17], q[22]],
(1, 1): [q[12], q[17], q[18], q[23]],
(1, 2): [q[13], q[18], q[19], q[24]],
(1, 3): [q[14], q[19], q[20], q[25]],
(1, 4): [q[15], q[20], q[21], q[26]],
(2, 0): [q[22], q[27], q[28], q[33]],
(2, 1): [q[23], q[28], q[29], q[34]],
(2, 2): [q[24], q[29], q[30], q[35]],
(2, 3): [q[25], q[30], q[31], q[36]],
(2, 4): [q[26], q[31], q[32], q[37]],
(3, 0): [q[33], q[38], q[39], q[44]],
(3, 1): [q[34], q[39], q[40], q[45]],
(3, 2): [q[35], q[40], q[41], q[46]],
(3, 3): [q[36], q[41], q[42], q[47]],
(3, 4): [q[37], q[42], q[43], q[48]],
(4, 0): [q[44], q[49], q[50], q[55]],
(4, 1): [q[45], q[50], q[51], q[56]],
(4, 2): [q[46], q[51], q[52], q[57]],
(4, 3): [q[47], q[52], q[53], q[58]],
(4, 4): [q[48], q[53], q[54], q[59]],
}
print(f'edge = {edge}')
#各格子点から出る手の本数が０本または２本
Hconst1 = 0
Hconst1 += ( q[0] + q[5] - 2*s[0])**2
Hconst1 += (q[0] + q[1] + q[6] - 2*s[1])**2
Hconst1 += (q[1] + q[2] + q[7] - 2*s[2])**2
Hconst1 += (q[2] + q[3] + q[8] - 2*s[3])**2
Hconst1 += (q[3] + q[4] + q[9] - 2*s[4])**2
Hconst1 += (q[4] + q[10] - 2*s[5])**2
Hconst1 += (q[5] + q[11] + q[16] - 2*s[6])**2
Hconst1 += (q[6] + q[11] + q[12] + q[17] - 2*s[7])**2
Hconst1 += (q[7] + q[12] + q[13] + q[18] - 2*s[8])**2
Hconst1 += (q[8] + q[13] + q[14] + q[19] - 2*s[9])**2
Hconst1 += (q[9] + q[14] + q[15] + q[20] - 2*s[10])**2
Hconst1 += (q[10] + q[15] + q[21] - 2*s[11])**2
Hconst1 += (q[16] + q[22] + q[27] - 2*s[12])**2
Hconst1 += (q[17] + q[22] + q[23] + q[28] - 2*s[13])**2
Hconst1 += (q[18] + q[23] + q[24] + q[29] - 2*s[14])**2
Hconst1 += (q[19] + q[24] + q[25] + q[30] - 2*s[15])**2
Hconst1 += (q[20] + q[25] + q[26] + q[31] - 2*s[16])**2
Hconst1 += (q[21] + q[26] + q[32] - 2*s[17])**2
Hconst1 += (q[27] + q[33] + q[38] - 2*s[18])**2
Hconst1 += (q[28] + q[33] + q[34] + q[39] - 2*s[19])**2
Hconst1 += (q[29] + q[34] + q[35] + q[40] - 2*s[20])**2
Hconst1 += (q[30] + q[35] + q[36] + q[41] - 2*s[21])**2
Hconst1 += (q[31] + q[36] + q[37] + q[42] - 2*s[22])**2
Hconst1 += (q[32] + q[37] + q[43] - 2*s[23])**2
Hconst1 += (q[38] + q[44] + q[49] - 2*s[24])**2
Hconst1 += (q[39] + q[44] + q[45] + q[50] - 2*s[25])**2
Hconst1 += (q[40] + q[45] + q[46] + q[51] - 2*s[26])**2
Hconst1 += (q[41] + q[46] + q[47] + q[52] - 2*s[27])**2
Hconst1 += (q[42] + q[47] + q[48] + q[53] - 2*s[28])**2
Hconst1 += (q[43] + q[48] + q[54] - 2*s[29])**2
Hconst1 += (q[49] + q[55] - 2*s[30])**2
Hconst1 += (q[50] + q[55] + q[56] - 2*s[31])**2
Hconst1 += (q[51] + q[56] + q[57] - 2*s[32])**2
Hconst1 += (q[52] + q[57] + q[58] - 2*s[33])**2
Hconst1 += (q[53] + q[58] + q[59] - 2*s[34])**2
Hconst1 += (q[54] + q[59] - 2*s[35])**2#直線接続ならプチボーナス
Hcost = 0
Hcost -= q[0]*q[1] + q[1]*q[2] + q[2]*q[3] + q[3]*q[4]
Hcost -= q[11]*q[12] + q[12]*q[13] + q[13]*q[14] + q[14]*q[15]
Hcost -= q[22]*q[23] + q[23]*q[24] + q[24]*q[25] + q[25]*q[26]
Hcost -= q[33]*q[34] + q[34]*q[35] + q[35]*q[36] + q[36]*q[37]
Hcost -= q[44]*q[45] + q[45]*q[46] + q[46]*q[47] + q[47]*q[48]
Hcost -= q[55]*q[56] + q[56]*q[57] + q[57]*q[58] + q[58]*q[59]
Hcost -= q[5]*q[16] + q[16]*q[27] + q[27]*q[38] + q[38]*q[49]
Hcost -= q[6]*q[17] + q[17]*q[28] + q[28]*q[39] + q[39]*q[50]
Hcost -= q[7]*q[18] + q[18]*q[29] + q[29]*q[40] + q[40]*q[51]
Hcost -= q[8]*q[19] + q[19]*q[30] + q[30]*q[41] + q[41]*q[52]
Hcost -= q[9]*q[20] + q[20]*q[31] + q[31]*q[42] + q[42]*q[53]
Hcost -= q[10]*q[21] + q[21]*q[32] + q[32]*q[43] + q[43]*q[54]#サンプリング
result = solver.run(qubo, shots=20000, T_num=200000)num1 = 60
num2 = 36
offset = 64.0
MODE: GPU
DEVICE: cuda:0
Energy -64.08000183105469, Occurrence 10
Energy -64.05999755859375, Occurrence 4
Energy -63.100006103515625, Occurrence 7
（以下略）
おお、できてますね。正解が-64.08、離れ小島が-64.06となり0.02の差がついています。たしかに、正解から直線を３本捨てて１本増やすと惜しい解になるので理論通りのエネルギーです。
例題２
こちらの例題も試してみましょう。中央が正解で、離れ小島の解がたくさんあります。
#問題
size = [5, 5]
rule = {
(0, 1): 2,
(0, 2): 1,
(0, 3): 2,
(0, 4): 3,
(2, 3): 1,
(2, 4): 3,
(3, 2): 2,
(4, 3): 2,
(4, 4): 0,
}num1 = 60
num2 = 36
offset = 36.0
MODE: GPU
DEVICE: cuda:0
Energy -36.10000228881836, Occurrence 2
Energy -36.10000228881836, Occurrence 3
Energy -36.08000183105469, Occurrence 1
（以下略）
やったか…！？（）
ほとんどの離れ小島解はブロックできたのですが、第２位のパターンだけは同率で出てきてしまいました。やはり定式化が完璧ではなかったか・・・惜しい！
おわりに
まあ同率１位で出してからOpenCVで島の個数を数えればいいんですけどね～
QUBO一本でガッチリ正解を出すことはできませんでした。ほぼ解けたということで勘弁していただきましょう ( ˘ω˘ ; )