温度計パズルを極低温で解く

温度計パズル

Androidのゲームで温度計パズルというのを見つけました。

見たらなんとなくルールが分かると思いますが、
- 温度計をひっくり返して、赤または白にする。ただし、以下の制約を満たすこと
- 横方向の赤の数、縦方向の赤の数が左端、上端に書かれているのでそれに合わせること
- 温度計を赤にする場合は、先端の丸い所から順に赤にしなければならない(途中から赤になるとか、どこか途切れるとかはダメ)

やってたら、めんどくさくなってきて、いい方法ないかな、と思って、量子アニーリングを試してみることにしました。

量子アニーリング

量子アニーリングについては、以前私が書いた記事やBrainPadの方が書かれた記事などをご覧ください。

ざっくり言うと、以下の数式の、行列Qを作ると、Hの値が小さくなるようにxの値を勝手に調整してくれるのが量子アニーリングです。

H=\sum_{i,j}Q_{ij}x_ix_j

パット見は分かりにくいですが、Qを調節すると、xについての2次方程式が自在に作れます。(もちろん、1次方程式も作れます)

これでどうやって問題を解くのか?

まず、xの値で何を表現するかを決めます。今回の問題はかなり量子アニーリングに当てはめやすい問題でした。温度計が赤だったら1, 白だったら0と置くことができます。
(そう簡単にはいかない場合がほとんどです。今回はラッキーです)

5x5の場合、こんな風にxを割り当ててみましょう。

ビット配置

x0  x1  x2  x3  x4
x5  x6  x7  x8  x9
x10 x11 x12 x13 x14
x15 x16 x17 x18 x19
x20 x21 x22 x23 x24

今回のようなパズルを量子アニーリングで解くときの基本は、「制約をクリアしていたら、Hの値は増えないけれど、クリアしていなかったら、Hの値が増える」というような状況を作ることです。
そうなるようにQを作っていきます。
続いて見ていきましょう。

横方向、縦方向の数

数合わせは、量子アニーリングの定式化で非常によく使われる頻出テクニックです。
なので、少し丁寧めに解説します。

例えば、1列目の赤の数を2に合わせる制約式を考えましょう。
白は0、赤は1だったので、1列目の赤の数を数えるには、

$(x0 + x5 + x10 + x15 + x20)$

となります。それがちょうど2になるので、試しに2を引いてみましょう。

$(x0 + x5 + x10 + x15 + x20 - 2)$

この数式がちょうど0のとき、1列目の赤の数の合計が2になります。
量子アニーリングは数式を最小化するので、これだけだと、x0, x5, x10, x15, x20が全部ゼロのとき-2で最小になってしまいます。

ところが、この数式を2乗すると、x0, x5, x10, x15, x20の合計がちょうど2のときに最小になります。

$(x0 + x5 + x10 + x15 + x20 - 2)^2$

4行目の赤の数を4つにするのも同様です。4行目のビットを全部足し合わせて、4つにしたいので4を引いて、全体を2乗します。

$(x15 + x16 + x17 + x18 + x19 - 4)^2$

こんな感じで、数の制約を作ることができます。
全部の行、全部の列についてやると、5x5の場合は10個の2次方程式ができますが、これらを全部足し算して、ひとつの2次方程式にしてしまいます。
(全部の制約が満たされているときのみ、値がゼロになり、最小です)

温度計を丸い方から赤くしていく制約

問題を、以下の表のように表すことにします。

	2	2	3	2	2
1	o	←	←	o	↑
2	↓	↑	o	→	o
3	↓	o	←	←	o
4	↓	←	←	←	o
1	↓	o	→	→	→

oが温度計の丸で、矢印が棒の部分です。
温度計を丸い方から赤くするという制約は、言い換えると、
「矢印の位置に対応するビットが1なのに、矢印をひとつ逆方向に辿った位置に対応するビットが0なのはNG」となります。

例えば←oとあるとき、以下のような制約になります。

←	o	OK?
0	0	OK
0	1	OK
1	0	NG
1	1	OK

NGの場合のみエネルギーが高くなる状況を作るにはどうすればいいでしょう?
上の表とにらめっこして、NGのときのみ1になる論理式を考えます。
ビットを左からa, bと置くと、NOT(a) AND bがそれに当たります。

数式で書くと、NOT(a)は(1-a)に、ANDは掛け算になります。
つまり、NOT(a) AND bは(1-a)*bあるいはb - a*bになります。
(実際に数字を当てはめると、NGのときのみ1に、それ以外は0になっていることが確認できます)

矢印の数だけ、この制約ができますが、全部足し合わせちゃって大丈夫です。

PyQUBOを使ってみる

今回、PyQUBOを使ってやってみることにしました。

入力ファイル

問題を、以下に示す形式で標準入力から与えています。

1行目: 行数 列数
2行目: (最初に半角スペースを入れて)各列の赤の数
3行目以降: 1文字目に行内の赤の数、2文字目以降、矢印とoで温度計の配置

問題1

5 5
 22322
1o←←o↑
2↓↑o→o
3↓o←←o
4↓←←←o
1↓o→→→

問題2

5 5
 12232
2↑o→o↑
3↑↑o↓↑
3↑↑↓↓o
2ooo→→
0o→→→→

コード

以下のように解きました。
D-WAVEの実機を使うには、D-WAVE Leapにログインしてトークンを得て、コードに埋め込んでください。

solver.py

from pprint import pprint
from pyqubo import Binary, solve_qubo
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite

# このへんで問題の入力を受け取ってる
w, h = map(int, input().split())
vert = [int(c) for c in input().strip()]
horiz = []
board = []
for _ in range(h):
    line = input()
    horiz.append(int(line[0]))
    board.append(line[1:])

# PyQUBOでビットを作ってる
q = [Binary(f'q{i:>02}') for i in range(w * h)]

# 盤面の座標から対応するビットの番号を得る
def mapping(x, y):
    return x + y * w

# 縦方向の赤の数、横方向の赤の数に関する制約を作ってる。
# 全部足し算して、ほしい数を引いて、2乗
h1 = 0
for i, v in enumerate(horiz):
    h1 += (sum(q[mapping(j, i)] for j in range(w)) - v) ** 2
for i, v in enumerate(vert):
    h1 += (sum(q[mapping(i, j)] for j in range(h)) - v) ** 2

# 温度計の、丸から順に埋めていく制約を作ってる。
h2 = 0
for y, line in enumerate(board):
    for x, v in enumerate(line):
        dx = 0
        dy = 0
        if v == 'o':
            continue
        elif v == '↑':
            dy = -1
        elif v == '→':
            dx = 1
        elif v == '↓':
            dy = 1
        elif v == '←':
            dx = -1
        else:
            raise ValueError('Unknown value ' + v)
        h2 += q[mapping(x, y)] - q[mapping(x - dx, y - dy)] * q[mapping(x, y)]

# 作った2種類の制約を適当なバランスで混ぜ合わせる。
# なんとなく1:2にしたが、このへん調整したらもっと解けるようになるかも(?)
H = h1 + 2 * h2
# 作った制約を、行列Qに変換するのはPyQUBOがやってくれる
model = H.compile()
qubo, offset = model.to_qubo()

# PyQUBO内蔵シミュレータで解く
result_sim = solve_qubo(qubo)
# 以下、表示用に出力を整えている
result_sim = [result_sim[k] for k in sorted(result_sim.keys())]
print('Simulated annealing result')
print(' ' + ''.join(map(str, vert)))
for y, (hor, line) in enumerate(zip(horiz, board)):
    line = list(line)
    for x in range(w):
        if not result_sim[y * w + x]:
            line[x] = ' '
    print(str(hor) + ''.join(line))

# D-WAVEで解く

# D-WAVE使うときの設定
endpoint = 'https://cloud.dwavesys.com/sapi'
solver_name = 'DW_2000Q_2_1'
token = "ご自身のD-WAVEのトークンを入れてください"

print('Quantum annealing result')
sampler = EmbeddingComposite(DWaveSampler(endpoint=endpoint, token=token, solver=solver_name))
response = sampler.sample_qubo(qubo, num_reads=20)
result_dw = response.first.sample
# 以下、表示用に出力を整えている
result_dw = [result_dw[k] for k in sorted(result_dw.keys())]
print(' ' + ''.join(map(str, vert)))
for y, (hor, line) in enumerate(zip(horiz, board)):
    line = list(line)
    for x in range(w):
        if not result_dw[y * w + x]:
            line[x] = ' '
    print(str(hor) + ''.join(line))

結果

問題1

Simulated annealing result
 22322
1o    
2↓ o  
3  ←←o
4 ←←←o
1 o   
Quantum annealing result
 22322
1o    
2↓ o  
3  ←←o
4 ←←←o
1 o   

シミュレータ、実機ともに答えが合っていそうです。

問題2

Simulated annealing result
 12232
2 o o 
3  o↓↑
3  ↓↓o
2oo   
0     
Quantum annealing result
 12232
2 o o 
3  o↓↑
3  ↓↓o
2oo   
0     

これもクリア。

大きい問題に挑戦

できたので、サイズをデカくしてみましょう。

問題7x7

7 7
 3313123
0←←←←←←o
3o→→→→←o
2←←←←←←o
3↑oo↑o→↑
4↑↓↓↑oo↑
2o↓↓↑↓↓↑
2←o↓o↓↓o

7x7結果

Simulated annealing result
 3313123
0       
3o→    o
2     ←o
3 oo o  
4↑↓ ↑ o 
2o  ↑   
2   o  o
Quantum annealing result
 3313123
0       
3o     o
2     ←o
3 o  o  
4↑↓ ↑ o 
2o  ↑   
2   o  o

いけてるねー。

問題7x7その2

7 7
 2252534
6oo→→→→→
5↓o→→→→→
3←←←←oo↑
2↑↑oo↑↓o
3↑↑oo↑↓o
3↑↑↓↓o→↓
1←←←←←o↓

7x7その2結果

Simulated annealing result
 2252534
6oo→→→→ 
5↓o→→→  
3    oo↑
2  o   o
3  o ↑ o
3  ↓ o ↓
1     o 
Quantum annealing result
 2252534
6oo→→→→→
5↓o→→→→ 
3  ←←o  
2    ↑  
3  o ↑ o
3  ↓ o ↓
1      ↓

おやおや、D-WAVEの結果、怪しいですね。
アニーリングは運ゲーなんですが、5x5はだいたい合ってたのが、この辺から運ゲー要素が強まってきています。

もうちょっと増やしてみましょう。

問題8x8

8 8
 14154684
5↓↓o←←←←o
4↑↑↓oo→→→
4↑↑↓↓o→→→
5↑↑↓↓o→→→
4↑↑↓↓o↑o→
2↑↑↓↓↓↑↑↑
3↑↑↓↓↓↑↑↑
6ooo→→ooo

8x8結果

Simulated annealing result
 14154684
5   ←←←←o
4   oo→→ 
4   ↓o→→ 
5   ↓o→→→
4 ↑ ↓  o→
2 ↑    ↑ 
3 ↑   ↑↑ 
6ooo  ooo
Quantum annealing result
 14154684
5 ↓ ←←←←o
4   oo→→ 
4   ↓ →→ 
5 ↑ ↓o→→ 
4 ↑ ↓  o→
2    ↓ ↑ 
3     ↑↑↑
6oo   ooo

シミュレータは合っていそうです。
D-WAVEくん！！！！しっかりして！！！

9x9もやってみましょう。

問題9x9

9 9
 686663754
5↑oo→→→→→→
6o↓↑o→→→→→
3o↓↑o→→→→→
8↓↓↑oo↑o→→
6↓↓↑↓↓↑o↑o
7↓↓o↓↓↑↓↑↓
6↓↓o→↓↑↓↑↓
7↓o→→→o↓o↓
3←←←←←←←o↓

結果9x9

Simulated annealing result
 686663754
5↑oo→→    
6o↓ o→→→  
3o↓ o→    
8↓↓↑oo o→→
6↓↓↑ ↓ o o
7↓↓o ↓ ↓↑↓
6 ↓o→  ↓↑↓
7 o→→ o↓o 
3     ←←o 
Quantum annealing result
 686663754
5 oo→→→   
6o↓  → →→ 
3o↓    →→ 
8↓↓↑oo↑o→→
6 ↓↑↓↓↑o o
7↓↓o↓↓ ↓↑↓
6↓↓o→↓ ↓ ↓
7↓o→→→  o↓
3  ←   ←  

どっちもできてないやないか……

まとめ

QUBOを作って、解いてみたところ、小さい問題は解けたものの、大きい問題では運ゲー要素が高まり、なんとも微妙な結果になりました。

シミュレータよりも先に、実機の方がギブアップになった印象があります。
量子アニーリングのよさを信じている方にとっては「なんで?」と思われるかもしれませんが、大変残念ながら「実機で解けない問題がシミュレータでは解ける」という話は、この界隈では非常にあるあるです。

制約項を混ぜ合わせるバランスなどの調整が悪かったのかもしれませんし、実機自体がまだ発展途上なのかもしれませんし、D-WAVEに適した問題はこういうのではなく別のところにあるのかもしれません。また、シミュレータでも9x9のパズル問題を解けておらず、そもそもQUBOを使って解く方法自体が適切だったのか、という疑問も残ります。(量子アニーリングでは、QUBOかそれに等価なイジングモデルで解く必要がありますが、通常のコンピュータではQUBOを使って解く必要はありません)

量子アニーリングの今後の益々のご活躍とご発展に期待しましょう。