数独¶
このチュートリアルでは、Amplifyを用いたイジングマシンによる数独の解法について解説します。
数独の説明¶
数独（すうどく）は、以下のルールに従って$9\times9$のブロックに$1\sim9$の数字を入れるパズルです。
- 空いているマスに$1\sim9$のいずれかの数字を入れる
- 縦・横の各列、および$9\times9$のブロックの中に9個ある$3\times3$のブロックには重複した数字は入らない
まず、ヒントとしていくつかの数字が埋められた初期配置が与えられます。上記のルールに従って、空いているマスに入る数字を確定していくことでゲームを進めることができます。ゲームの難易度が低い場合は、比較的簡単に次々と入り得る数字を確定できるマスを見つけることができますが、ゲームの難易度が上がると、そのような数字の確定が難しくなり、ある程度の経験を積まないとパズルを解き進めるのは難しくなります。
コンピュータを使って数独を解く方法として、深さ優先探索、確率的な方法、制約充足問題、exact cover problemなど用いた様々なアルゴリズムが考案されています。これらのアルゴリズムに従って、機械的に数独を解くことができます。
このチュートリアルでは、組み合わせ最適化問題に特化したイジングマシンを使って数独を解く方法を紹介します。上記の数独のルールを制約条件として解釈し、それに伴ったコスト関数を定義し、コストが最も低い数字の組み合わせを見つけることで数独の解を見つけることができる仕組みになっています。よって、ここで行うべきことは数独のルールをどのような制約条件によって表すことができるかを考えることです。そのような制約条件を見つけてコスト関数を定義することができれば、あとは初期配置を与えるだけで、複雑なアルゴリズムを用いることなく、イジングマシンによって解が見つけることができます。
それでは次に、Amplifyを使って、数独を解くコードがどのように実装されるのか具体的に見てみましょう。
制約条件の定式化¶
アニーリングマシンにおける数独の制約条件の表現方法¶
次に、イジングマシンにおいて、数独のルールを表す制約条件を満たすコスト関数の作成方法を考えましょう。基本的には、二次二値多変数多項式を用いて制約条件を表現する方法を考えることとなります。ここでは、QUBO模型（各変数は0または1）を用いた一つの方法に着目して議論を進めます。
表すべき数独のルールは以下の4つとなります。
- 各マスには $1\sim9$ の数字のいずれかが入る
- 各行には $1\sim9$ の数字が 1 つずつ重複することなく入る
- 各列には $1\sim9$ の数字が 1 つずつ重複することなく入る
- 各 $3\times3$ のブロックには$1\sim9$の数字が 1 つずつ重複することなく入る
まず、ルール 0 「各マスには $1\sim9$ の数字のいずれかが入る」について、各マスに入っている数字をいくつかのバイナリ変数 ($0$ または $1$ の値を取る変数) を用いて表すことを考えます。
ある 1 マスについて、バイナリ変数が $9$ 個あれば、そのマスに何の数字が入っているかを表すことができます。それぞれのバイナリ変数は、「そのマスに $1$ が入っているかどうか」「そのマスに $2$ が入っているかどうか」$\cdots$「そのマスに $9$ が入っているかどうか」をそれぞれ表します。たとえば以下の表は、数字 $3$ が入っている場合の例です。
| 変数名 | 変数 0 | 変数 1 | 変数 2 | 変数 3 | 変数 4 | 変数 5 | 変数 6 | 変数 7 | 変数 8 |
|---|---|---|---|---|---|---|---|---|---|
| 数字 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
| 変数の値 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
$9\times9$ のマス目に対しても同様に考えて、$9\times9=81$ 個の変数を $9$ セット、合計で $9\times9\times9=729$ 個の変数を用意します。これは、$9\times9=81$ マスのレイヤを$9$枚重ねるようなイメージです。各レイヤは、「それぞれのマスに $1$ が入っているかどうか」「それぞれのマスに $2$ が入っているかどうか」$\cdots$「それぞれのマスに $9$ が入っているかどうか」を表します。$i$ 行 $j$ 列かつ $k$ 番目のレイヤにある変数を $q_{i, j, k}$ と書くことにします。ただし、インデックスはそれぞれ 0 から始まる ($i,j,k=0, 1,\ldots,8$) ものとします。
例えば、$3$ 行 $5$ 列のマスに $7$ が入る状態は $q_{2,4,6}=1$、そうでない場合は $q_{2,4,6}=0$ と表現することができます。
それぞれのマスについて、そのマスに対応する各レイヤの $9$ 個の変数のうちいずれかちょうど 1 つが $1$ である必要があることに注意してください。さもないとそのマスに入る数字が一意に定まらないためです。つまり、
$$ (a) \quad \sum_{k=0}^8 q_{i,j,k}=1 $$
が必要です。
次に、ルール 1～3 を上記のバイナリ変数を用いて表します。それぞれ以下のone-hot制約 (a)、(b)、(c) として書き下すことができます。
$$ \left\{ \begin{align*} (b) \quad &\sum_{j=0}^8 q_{i,j,k}=1 \\ (c) \quad &\sum_{i=0}^8 q_{i,j,k}=1 \\ (d) \quad &\sum_{i,j\in 3\times3\text{ブロック}}q_{i,j,k}=1 \\ \end{align*} \right. $$
たとえばルール 1 である「各行には $1\sim9$ の数字が 1 つずつ重複することなく入る」は、どの $1\sim 9$ の数字についても「各行にはその数字がちょうど 1 つ入る」と言い換えられます。これは、どのレイヤに対しても、そのレイヤの各行には $1$ がちょうど 1 つあるということなので、ルール 1 は one-hot 制約 (a) として書けることが分かります。ルール 2, 3 についても同様です。
初期配置¶
数独ではいくつかのマスがすでに埋められた初期配置がヒントとして与えられます。ここでは、以下の初期配置を使います。以下では、数字が埋められていないマスは$0$としました。
import numpy as np
# 初期配置をリストで表記
initial = np.array(
[
[5, 3, 0, 0, 7, 0, 0, 0, 0],
[6, 0, 0, 1, 9, 5, 0, 0, 0],
[0, 9, 8, 0, 0, 0, 0, 6, 0],
[8, 0, 0, 0, 6, 0, 0, 0, 3],
[4, 0, 0, 8, 0, 3, 0, 0, 1],
[7, 0, 0, 0, 2, 0, 0, 0, 6],
[0, 6, 0, 0, 0, 0, 2, 8, 0],
[0, 0, 0, 4, 1, 9, 0, 0, 5],
[0, 0, 0, 0, 8, 0, 0, 7, 9],
]
)
# 数独を成形して表示する
def print_sudoku(sudoku):
for i in range(len(sudoku)):
line = ""
if i == 3 or i == 6:
print("---------------------")
for j in range(len(sudoku[i])):
if j == 3 or j == 6:
line += "| "
line += str(sudoku[i][j]) + " "
print(line)
print_sudoku(initial)
from amplify import VariableGenerator
gen = VariableGenerator()
q = gen.array("Binary", shape=(9, 9, 9))
これによって、$9^3=729$ 個の変数が三次元配列として用意されました。9, 9, 9
はそれぞれ、行・列・レイヤの要素数を表し、それらのインデックスをi、j、kとして各要素にはq[i, j, k]でアクセスできます。
# 行番号 $1$、列番号 $2$、レイヤ $3$ の変数
print(q[1, 2, 3])
# 行番号 $0$、列番号 $0$ の9変数
print(q[0, 0])
# 行番号 $2$、数値レイヤ $5$ の9変数
print(q[2, :, 5])
先ほど initial に格納された初期配置から、いくつかの変数をただちに確定することができます。例えば、i=1、j=4
のマスにはすでに数字 9 が入っているので (initial[1, 4] == 9)、対応する変数 q[1, 4, 8] は
$1$
に確定できます。インデックスが 0 から始まる都合上、数字 9 に対応するレイヤはレイヤ $8$ であることに注意してください。
また、i=1、j=4 のマスには数字 9 以外の数字は入らないので、i=1
の行、j=4 の列でレイヤ番号が $8$ 以外であるような変数 (たとえば q[1, 4, 7]) の値は $0$ に確定できます。
さらに、ルール 1 により、i=1 の行の j=4 以外のマスには数字 9 以外は入らないので、
i=1
の行の k=8 のレイヤで行番号 $j$ が $4$ 以外であるような変数 (たとえば q[1, 5, 8]) の値は $0$ に確定できます。
同様に考えて、ルール 2、ルール 3 により、k=8 のレイヤで i=1、j=4 のマスと同じ列あるいは同じ
$3\times 3$
のブロックに所属するマスに対応する変数の値は $0$ に確定できます。たとえば q[0, 4, 8] や q[2, 5, 8] の値は $0$ です。
for i, j in zip(*np.where(initial != 0)):
# 初期配置の i 行 j 列には数字が既に入っている
k = initial[i, j] - 1
q[i, j, :] = 0 # i 行 j 列のマスには数字 k+1 以外は入らない
q[i, :, k] = 0 # 同じ行には数字 k+1 は入らない
q[:, j, k] = 0 # 同じ列には数字 k+1 は入らない
for m in range(9):
# 同じブロックには数字 k+1 は入らない
q[(3 * (i // 3) + m // 3), (3 * (j // 3) + m % 3), k] = 0
q[i, j, k] = 1 # i 行 j 列のマスには数字 k+1 が入る
これで初期設定ができました。例として行番号 $0$、列番号 $0$ の 9 変数を表示すると5番目の要素が 1 に、それ以外が 0 に確定されていることが確認できます。これは左上のマスが 5 であることに対応します。
q[0, 0]
同様に行番号 $0$、列番号 $2$ の9変数を表示すると、3 番目の変数と 5 番目以降の変数が 0 であることが分かります。つまり、一番上の行の左から 3 番目のマスには、1, 2, 4 のいずれかの数字が入ります。
q[0, 2]
制約条件の設定¶
次に、制約条件 $(a)\sim(d)$ を定義します。これらはすべていくつかのバイナリ変数のうち 1 つだけが 1 である one-hot 制約条件となっています。
まず、$(a)$ の一つのマスには一つの数字しか入らないという制約を表す制約条件を作成します。
$$ (a) \quad \sum_{k=0}^8 q_{i,j,k}=1 $$
変数配列 q の行 (axis=0) と列 (axis=1) が同じでレイヤ (axis=2)
が異なる
$9$ つの変数のうち $1$ つのみが $1$ である制約なので、one_hot 関数の axis パラメータに 2 を与えることで、$9\times
9$
個の制約条件を一気に作成することができます。
from amplify import one_hot
# (a): 一つのマスには一つの数字しか入らない制約条件
layer_constraints = one_hot(q, axis=2)
print(len(layer_constraints))
同様にして、$(b)$ の同じ行に同じ数字が入らない制約条件と、$(c)$ の同じ行に同じ数字が入らない制約条件は以下のように表されます。
$(b)$ は列 (axis=1) とレイヤ (axis=2) が同じで行 (axis=0) が異なる $9$
つの変数のうち
$1$ つのみが $1$ である制約であり、 $(c)$ は行 (axis=0) とレイヤ (axis=2) が同じで列
(axis=1) が異なる $9$ つの変数のうち $1$ つのみが $1$ である制約なので、one-hot 関数の
axis
パラメータにそれぞれ 1 と 0 を与えます。
# (b): 各行には同じ数字が入らない制約条件
row_constraints = one_hot(q, axis=1)
# (c): 各列には同じ数字が入らない制約条件
col_constraints = one_hot(q, axis=0)
最後に$(c)$ の $3\times3$ の各ブロックには同じ数字が入らないという制約条件を表します。各レイヤごとに、各 $3\times3$ ブロック内で変数の和を取り、以下のようにして one-hot 制約を課します。
from amplify import sum as amplify_sum
# (c): 3x3ブロック内には同じ数字が入らない制約条件
block_constraints = amplify_sum(
one_hot(amplify_sum([q[i + m // 3, j + m % 3, k] for m in range(9)]))
for i in range(0, 9, 3)
for j in range(0, 9, 3)
for k in range(9)
)
これで全ての制約条件が出そろったので、これらの制約条件を全て足し合わせます。
constraints = layer_constraints + row_constraints + col_constraints + block_constraints
これで定式化に関する準備ができました。イジングマシンによって、全ての制約を満たす変数の組み合わせを見つけ出すことが出来れば、そのような変数の組み合わせは与えられた初期配置から導き出される数独の解となります。
イジングマシンの実行¶
先ほど作成した constraints を用いてイジングマシンを実行する準備を行います。まずイジングマシンのクライアントを作成し、timeout
などのパラメーターを設定します。その後、ソルバーを作成しクライアントを設定します。
from amplify import solve, FixstarsClient
from datetime import timedelta
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=1000) # タイムアウト 1000 ms
# ローカル環境等で使用する場合は、コメントを外して Fixstars Amplify AE のアクセストークンを入力してください
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
# 求解の実行
result = solve(constraints, client)
if len(result) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
実行結果は result に格納されています。以下のようにすることで求解結果から変数配列 q の値が得られます。
q_values = q.evaluate(result.best.values).astype(int)
print(q_values[0, 0, 4]) # 左上のマスが 5 であるか
print(q_values[4, 4, 0]) # 5 行 5 列のマス (盤面のちょうど真ん中のマス) が 1 であるか
解から盤面を復元します。あるマスに対応する 9 変数 について、各変数は「そのマスに $1$ が入っているかどうか」「そのマスに $2$ が入っているかどうか」$\cdots$「そのマスに $9$
が入っているかどうか」を表すので、そのマスの値を知るためには 9 つの変数の値と長さ 9 の 1 次元ベクトル [1, 2, 3, 4, 5, 6, 7, 8, 9]
との内積をとればよいです。
つまり、盤面は以下の式で2次元配列として表されます。
$$ \text{answer}[i, j] = \sum_{k=0}^{8} \text{q\_values}[i, j, k] \cdot k $$
answer = q_values @ np.arange(1, 10)
最後に print_sudoku 関数を用いて解答を出力します。正しい解が得られたことが分かります。
print_sudoku(answer)
