量子アニーリング（HOBO）でビルディングパズルを解いてみた

量子コンピュータ

2025/01/11

やること
TYTANについて
HOBOTANについて
ビルディングパズル（Skyscrapers）とは
定式化と予備実験
パズルへの拡張
もう一問
おわりに

やること

まだまだ行きます。HOBO対応の疑似アニーリングパッケージ「HOBOTAN」でペンシルパズルを解くシリーズ。

今回は「ビルディングパズル（Skyscrapers）」を解いてみましょう。非常に難しかったです。

※この記事はQUBOアニーリングの基礎講座を履修済みであることを前提としています。

TYTANについて

前身となる「TYTAN」はOSSの疑似量子アニーリングパッケージです。

GitHub - tytansdk/tytan: Python SDK for large QUBO problems

Python SDK for large QUBO problems. Contribute to tytansdk/tytan development by creating an account on GitHub.

日本量子コンピューティング協会が運営しているチュートリアル教材はこちら。

GitHub - tytansdk/tytan_tutorial

Contribute to tytansdk/tytan_tutorial development by creating an account on GitHub.

私のオンライン講義でよければこちらから見れます。

【第3世代】1/5　QUBOアニーリング　TYTANチュートリアル

discordコミュニティもあるので気軽に質問できます。

https://discord.gg/qT5etstPW8

HOBOTANについて

それを元にしてできたのが「HOBOTAN」。使い方はTYTANとほとんど同じです。分からないことがあれば何でも聞いてください。

GitHub - ShoyaYasuda/hobotan: tytanのHOBO版（テスト）

tytanのHOBO版（テスト）. Contribute to ShoyaYasuda/hobotan development by creating an account on GitHub.

ビルディングパズル（Skyscrapers）とは

またまた Puzzle Team さんから問題を拝借しました。いつもありがとうございます。

16個のマスを１～４の数字で埋めます。どの行も列も数字が重複してはいけません。マスに書き込む数字はビルの高さを表していて、ある方向から見たときにいくつのビルが見えるかが指定されています。例えば手前から順に３，２，４，１と書き込んだ場合、３，４のビルは見えますが２，１のビルは隠れて見えません。この場合は２棟が見えることになります。

答えはこちらです。

コンセプティスパズルさんでも遊べます。

定式化と予備実験

長らくQUBOで定式化できず困っていました。まぢムリ。。２次制約で収まる方法があるのでしょうか？もうね、諦めてHOBOで取り組むことにしました。

まずは簡略化して４つのビルを考えてみます。4×4の量子ビットを用意して次のように対応させます。

まず、一番手前の列（１列目）のビルが必ず見えるのは自明ですね。

問題は２列目以降です。２列目は次のように考えました。簡単に言うと「自分の前にあるビルが自分より低ければ見えるよね」です。もし「自分の高さが４」かつ「１列目のビルの高さが１～３」であれば見えます。または、「自分の高さが３」かつ「１列目のビルの高さが１～２」でも見えます。または、「自分の高さが２」かつ「１列目のビルの高さが１」でも見えます。この３パターンの”見える”を合計したものが、２列目のビルが見えるかどうかの判定（０or１）となります。

４列目についても考えてみましょう。発想は同じで、「自分の前にあるすべてのビルが自分より低ければ見えるよね」です。判定条件が長くなります。もし「自分の高さが４」かつ「１列目のビルの高さが１～３」かつ「２列目のビルの高さが１～３」かつ「３列目のビルの高さが１～３」であれば見えます。以下同文。

定式化において「かつ」はかけ算でした。４列目は「A かつ B かつ C かつ D」が登場するので４次の制約になります。（実際にはさらに２乗するので８次になる）

ということで実装してみます。解は理想的なエネルギーのだけ表示するようにしています。

from hobotan import *
import numpy as np
import matplotlib.pyplot as plt

#量子ビットの用意
q = symbols_list([4, 4], 'q{}_{}')

#式を用意
H = 0

#各行に１個、各列に１個
for i in range(4):
    H += (np.sum(q[i, :]) - 1)**2
for j in range(4):
    H += (np.sum(q[:, j]) - 1)**2

#手前から１番目のビルが見えるかどうか（見える：１、見えない：０）
num0 = 1

#手前から２番目のビルが見えるかどうか
num1 = 0
num1 += q[0, 1] * np.sum(q[1:, 0])
num1 += q[1, 1] * np.sum(q[2:, 0])
num1 += q[2, 1] * np.sum(q[3:, 0])

#手前から３番目のビルが見えるかどうか
num2 = 0
num2 += q[0, 2] * np.sum(q[1:, 0]) * np.sum(q[1:, 1])
num2 += q[1, 2] * np.sum(q[2:, 0]) * np.sum(q[2:, 1])
num2 += q[2, 2] * np.sum(q[3:, 0]) * np.sum(q[3:, 1])

#手前から４番目のビルが見えるかどうか
num3 = 0
num3 += q[0, 3] * np.sum(q[1:, 0]) * np.sum(q[1:, 1]) * np.sum(q[1:, 2])
num3 += q[1, 3] * np.sum(q[2:, 0]) * np.sum(q[2:, 1]) * np.sum(q[2:, 2])
num3 += q[2, 3] * np.sum(q[3:, 0]) * np.sum(q[3:, 1]) * np.sum(q[3:, 2])

#見えるビルの数
H += (num0 + num1 + num2 + num3 - 4)**2

#HOBOテンソルにコンパイル
hobo, offset = Compile(H).get_hobo()
print(f'offset\n{offset}')

#サンプラー選択
solver = sampler.MIKASAmpler()

#サンプリング
result = solver.run(hobo, shots=100)

#理想的な解のみ
for r in result:
    print(f'Energy {r[1]}, Occurrence {r[2]}')
    
    if (r[1] + offset)**2 < 0.1:

        #配列取得
        arr, subs = Auto_array(r[0]).get_ndarray('q{}_{}')
        
        #さくっと画像に
        img, subs = Auto_array(r[0]).get_image('q{}_{}')
        plt.figure(figsize=(3, 3))
        plt.imshow(img)
        plt.show()

tensor order = 8
----- compile -----
MODE: GPU
DEVICE: cuda:0
-------------------
offset
17.0
===== sampler =====
MODE: GPU
DEVICE: cuda:0
===================
.......... 100/2000 min=-15.0 mean=2.819999933242798
.......... 200/2000 min=-16.0 mean=2.069999933242798
.......... 300/2000 min=-16.0 mean=14.469999313354492
（中略）
.......... 1800/2000 min=-17.0 mean=14.579999923706055
.......... 1900/2000 min=-17.0 mean=13.639999389648438
.......... 2000/2000 min=-17.0 mean=13.269999504089355
Energy -17.0, Occurrence 15

やはり８次の式になったようです。「４棟見える」にしたので以下の１パターンが出ました。４棟見えるのはこのパターンしかないですね。

「３棟見える」だと以下の６通りが出ました。

「２棟見える」だと以下の11通り。

「１棟見える」は以下の６通り。

網羅できていると思います。全部で4!＝24通り。この予備実験の段階で必要十分な解が出ていることを入念に確認しておく必要があります。「エネルギーは理想的なのに間違った解が出る」とか「正しい解なのにエネルギーが理想にならない」ということがあってはいけません。そのまま先に進んでも永久にパズルは解けません。

パズルへの拡張

問題を再掲。

まず、数独と同じように4x4x4の量子ビットを用意しましょう（→量子アニーリング（QUBO）でナンプレ（数独）と不等号ナンプレを解く）。

数独と同じように「どの直線の４個もワンホット」にします。これの意味は、「ある行には高さNのビルが必ず一つ」とか「ある列には高さNのビルが必ず一つ」とか「あるマスに立つビルは必ず一つ」みたいなことです。

あとは周囲から見えるビルの数を設定していくのですが、予備実験のやつを関数化して、３次元配列から切り出した２次元配列を渡してやることにします。切り出す方向が難しいですね。うお～～（paizaで鍛えたPython力がうなりを上げる）

from hobotan import *
import numpy as np
import matplotlib.pyplot as plt

#量子ビットの用意
q = symbols_list([4, 4, 4], 'q{}_{}_{}')

#式を用意
H = 0

#各次元に超絶ワンホット
for i in range(4):
    for j in range(4):
        H += (np.sum(q[i, j, :]) - 1)**2
for i in range(4):
    for k in range(4):
        H += (np.sum(q[i, :, k]) - 1)**2
for j in range(4):
    for k in range(4):
        H += (np.sum(q[:, j, k]) - 1)**2

def visible_num(s):
    '''
    4x4の量子ビット配列を入力
    見える数 num0 + num1 + num2 + num3 を返す
    '''
    #手前から１番目のビルが見えるかどうか（見える：１、見えない：０）
    num0 = 1
    
    #手前から２番目のビルが見えるかどうか
    num1 = 0
    num1 += s[0, 1] * np.sum(s[1:, 0])
    num1 += s[1, 1] * np.sum(s[2:, 0])
    num1 += s[2, 1] * np.sum(s[3:, 0])
    
    #手前から３番目のビルが見えるかどうか
    num2 = 0
    num2 += s[0, 2] * np.sum(s[1:, 0]) * np.sum(s[1:, 1])
    num2 += s[1, 2] * np.sum(s[2:, 0]) * np.sum(s[2:, 1])
    num2 += s[2, 2] * np.sum(s[3:, 0]) * np.sum(s[3:, 1])
    
    #手前から４番目のビルが見えるかどうか
    num3 = 0
    num3 += s[0, 3] * np.sum(s[1:, 0]) * np.sum(s[1:, 1]) * np.sum(s[1:, 2])
    num3 += s[1, 3] * np.sum(s[2:, 0]) * np.sum(s[2:, 1]) * np.sum(s[2:, 2])
    num3 += s[2, 3] * np.sum(s[3:, 0]) * np.sum(s[3:, 1]) * np.sum(s[3:, 2])
    
    return num0 + num1 + num2 + num3

#問題の反映
#左から見る
aim = [1, 3, 2, 3]
for i in range(4):
    q_tmp = q[i, :, :].T
    H += (visible_num(q_tmp) - aim[i])**2

#右から見る
aim = [4, 1, 2, 2]
for i in range(4):
    q_tmp = q[i, :, :].T[:, ::-1]
    H += (visible_num(q_tmp) - aim[i])**2

#上から見る
aim = [1, 2, 3, 2]
for j in range(4):
    q_tmp = q[:, j, :].T
    H += (visible_num(q_tmp) - aim[j])**2

#下から見る
aim = [3, 2, 1, 2]
for j in range(4):
    q_tmp = q[:, j, :].T[:, ::-1]
    H += (visible_num(q_tmp) - aim[j])**2

#HOBOテンソルにコンパイル
hobo, offset = Compile(H).get_hobo()
print(f'offset\n{offset}')

#サンプラー選択
solver = sampler.MIKASAmpler()

#サンプリング
result = solver.run(hobo, shots=200, T_num=10000)

#上位3件
for r in result[:3]:
    print(f'Energy {r[1]}, Occurrence {r[2]}')

    #配列取得
    arr, subs = Auto_array(r[0]).get_ndarray('q{}_{}_{}')
    
    #ワンホットチェック
    check = np.sum(arr, axis=2)
    print(f'check\n{check}')
    
    #解答
    ans = 4 - np.argmax(arr, axis=2)
    print(f'ans\n{ans}')

tensor order = 8
----- compile -----
MODE: GPU
DEVICE: cuda:0
-------------------
offset
80.0
===== sampler =====
MODE: GPU
DEVICE: cuda:0
===================
.......... 200/10000 min=-24.0 mean=62.26499938964844
.......... 400/10000 min=-25.0 mean=72.12999725341797
.......... 600/10000 min=-28.0 mean=66.61000061035156
（中略）
.......... 9600/10000 min=-80.0 mean=204.5349884033203
.......... 9800/10000 min=-80.0 mean=182.93499755859375
.......... 10000/10000 min=-80.0 mean=162.47999572753906
Energy -80.0, Occurrence 2
check
[[1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]]
ans
[[4 3 2 1]
 [2 1 3 4]
 [3 4 1 2]
 [1 2 4 3]]
Energy -77.0, Occurrence 1
check
[[1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]]
ans
[[4 3 2 1]
 [3 2 1 4]
 [1 4 3 2]
 [2 1 4 3]]
Energy -76.0, Occurrence 1
check
[[1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]]
ans
[[4 3 2 1]
 [2 1 3 4]
 [3 4 1 2]
 [1 3 4 2]]

上位３件の解を表示しています。チェックとして各マスのビルが一つだけになっているかどうかを確認しつつ、答えの数字に戻しています。

#１番目の解
[[4 3 2 1]
 [2 1 3 4]
 [3 4 1 2]
 [1 2 4 3]]

一致していますね。お見事！

もう一問

コンセプティスパズルさんの問題も解いてみましょう。

tensor order = 8
----- compile -----
MODE: GPU
DEVICE: cuda:0
-------------------
offset
54.0
===== sampler =====
MODE: GPU
DEVICE: cuda:0
===================
.......... 200/20000 min=-28.0 mean=17.98499870300293
.......... 400/20000 min=-28.0 mean=25.154998779296875
.......... 600/20000 min=-28.0 mean=27.93000030517578
（中略）
.......... 19600/20000 min=-54.0 mean=36.70000076293945
.......... 19800/20000 min=-54.0 mean=30.65999984741211
.......... 20000/20000 min=-54.0 mean=30.049999237060547
Energy -54.0, Occurrence 2
check
[[1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]]
ans
[[2 1 4 3]
 [3 4 1 2]
 [4 2 3 1]
 [1 3 2 4]]