Sudoku

このチュートリアルでは、Amplifyを用いたイジングマシンによる数独の解法について解説します。

数独の説明

数独(すうどく) は、以下のルールに従って \(9\times9\) のブロックに \(1\sim9\) の数字を入れるパズルです。

  • 空いているマスに \(1\sim9\) のいずれかの数字を入れる

  • 縦・横の各列、および \(9\times9\) のブロックの中に9個ある \(3\times3\) のブロックには重複した数字は入らない

まず、ヒントとして17個以上の数字が埋められた初期配置が与えられます (ヒントが16個以下の初期配置は解法を持ちえないことが証明されています)。上記のルールに従って、空いているマスに入る数字を確定していくことでゲームを進めることができます。ゲームの難易度が低い場合は、比較的簡単に次々と入り得る数字を確定できるマスを見つけることができますが、ゲームの難易度が上がると、そのような数字の確定が難しくなり、ある程度の経験を積まないとパズルを解き進めるのは難しくなります。

https://upload.wikimedia.org/wikipedia/commons/thumb/e/e0/Sudoku_Puzzle_by_L2G-20050714_standardized_layout.svg/361px-Sudoku_Puzzle_by_L2G-20050714_standardized_layout.svg.png

コンピュータを使って数独を解く方法として、深さ優先探索、確率的な方法、制約充足問題、exact cover problemなど用いた 様々なアルゴリズム が考案されています。これらのアルゴリズムに従って、機械的に数独を解くことができます。

このチュートリアルでは、組み合わせ最適化問題に特化したイジングマシンを使って数独を解く方法を紹介します。上記の数独のルールを制約条件として解釈し、それに伴ったコスト関数を定義し、コストが最も低い数字の組み合わせを見つけることで数独の解を見つけることができる仕組みになっています。よって、ここで行うべきことは数独のルールをどのような制約条件によって表すことができるかを考えることです。そのような制約条件を見つけてコスト関数を定義することができれば、あとは初期配置を与えるだけで、複雑なアルゴリズムを用いることなく、イジングマシンによって解が見つけることができます。

それでは次に、Amplifyを使って、数独を解くコードがどのように実装されるのか具体的に見てみましょう。

制約条件の定式化

アニーリングマシンにおける数独の制約条件の表現方法

次に、イジングマシンにおいて、数独のルールを表す制約条件を満たすコスト関数の作成方法を考えましょう。基本的には、二次二値多変数多項式を用いて制約条件を表現する方法を考えることとなります。ここでは、QUBO模型(各変数は0または1)を用いた一つの方法に着目して議論を進めます。

表すべき数独のルールは以下の3つとなります。

  1. 各行には \(1\sim9\) の数字が重複することなく入る

  2. 各列には \(1\sim9\) の数字が重複することなく入る

  3. \(3\times3\) のブロックには \(1\sim9\) の数字が重複することなく入る

まず、\(9\times9=81\) 個の各マスに、\(0\)\(1\) に値を取る変数を与えることを考えます。行と列を表すインデックスを \(i,j=0,\cdots,8\) とし、\(1,\cdots, 9\) 番目の行と列に対応させます。

1~3の制約を課さずに、数字のあらゆる重複を許す場合を考えると、全 \(9\times9 = 81\) マスには9個の数字が入り得ます。そこで、\(9\times9=81\) 個の変数を \(9\) セット考慮し、合計で \(9\times9\times9=729\) 個の変数を取り扱うことを考えます。これは、\(9\times9=81\) マスのレイヤを \(9\) 枚重ねるようなイメージです。ここで、レイヤのインデックスを \(k = 0\sim8\) とし、それぞれ数字の \(1\sim9\) に対応させます。行、列、レイヤのインデックスをそれぞれ \(i,j,k=0,\cdots,8\) とし、変数を \(q_{i,j,k}\) で表すと、例えば、\(3\) 列のマスに \(7\) が入る状態は \(q_{2,4,6}=1\)、そうでない場合は \(q_{2,4,6}=0\) と表現することができます。これらの変数を使うと、制約条件1~3はそれぞれ以下のone-hot制約として書き下すことができます。

\[\begin{split}\left\{ \begin{align} &\begin{split} (a) \quad &\sum_{j=0}^8 q_{i,j,k}=1 \end{split}\\ &\begin{split} (b) \quad &\sum_{i=0}^8 q_{i,j,k}=1 \end{split}\\ &\begin{split} (c) \quad &\sum_{i,j\in 3\times3\text{ブロック}}q_{i,j,k}=1 \end{split}\\ &\begin{split} (d) \quad &\sum_{k=0}^8 q_{i,j,k}=1 \end{split} \end{align} \right.\end{split}\]

制約条件 \((a)\)\((b)\)\((c)\) はそれぞれルール1、2、3に対応します。\((d)\) の制約条件は各マスには一つの数字しか入らないという基本的な条件に対応します。

初期配置

数独では、17個以上のいくつかのマスがすでに埋められた初期配置がヒントとして与えられます。ここでは、難しい問題とされる以下の初期配置を使います。以下の表記では、数字が埋められていないマスは \(0\) としました。

import numpy as np

# 初期配置をリストで表記
initial = np.array(
    [
        [2, 0, 5, 1, 3, 0, 0, 0, 4],
        [0, 0, 0, 0, 4, 8, 0, 0, 0],
        [0, 0, 0, 0, 0, 7, 0, 2, 0],
        [0, 3, 8, 5, 0, 0, 0, 9, 2],
        [0, 0, 0, 0, 9, 0, 7, 0, 0],
        [0, 0, 0, 0, 0, 0, 4, 5, 0],
        [8, 6, 0, 9, 7, 0, 0, 0, 0],
        [9, 5, 0, 0, 0, 0, 0, 3, 1],
        [0, 0, 4, 0, 0, 0, 0, 0, 0],
    ]
)

# 数独を成形して表示する
def print_sudoku(sudoku):
    print("\n\n")
    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)
2 0 5 | 1 3 0 | 0 0 4
0 0 0 | 0 4 8 | 0 0 0
0 0 0 | 0 0 7 | 0 2 0
---------------------
0 3 8 | 5 0 0 | 0 9 2
0 0 0 | 0 9 0 | 7 0 0
0 0 0 | 0 0 0 | 4 5 0
---------------------
8 6 0 | 9 7 0 | 0 0 0
9 5 0 | 0 0 0 | 0 3 1
0 0 4 | 0 0 0 | 0 0 0

制約条件の作成

変数の定義と初期配置の反映 まずは、Amplifyで提供されている gen_symbolsBinaryPoly を用いて、変数を用意します。

これによって、\(9^3=729\) 個の変数が三次元配列として用意されました。9, 9, 9 はそれぞれ、行・列・数値レイヤの要素数を表し、それらのインデックスを ijk として各要素には q[i][j][k] でアクセスできます。例えば 行番号 \(0\) 、列番号 \(0\) の9変数を表示する場合は次のようにします。

>>> q[0][0]
[q_0, q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8]

先ほど initial に格納された初期配置から数独のルールに従い、確定可能な未知変数の絞り込みを行います。例えば、i=1j=5 のマスにはすでに 8 (initial[1][5]=8) が入っているので、8 に対応した k=7 のレイヤの変数は q[1][5][7]=1 と指定されます。

ルール1、2により、k=7 のレイヤで i=1j=5 に対応したマスが属する行と列には同じ数字が入らないので、q[i][5][7]=0 (\(i\neq1\))、q[1][j][7]=0 (\(j\neq5\)) と変数の値をさらに確定させることができます。これは制約 \((a)\)\((b)\) を課すことにに対応します。

また、ルール3により、この数字が属する \(3\times3\) ブロック内で同じ数字は入らないので、\((i,j)\in\{(0,3), (0,4), (0,5), (1,3), (1,4), (2,3), (2,4), (2,5)\}\) において、q[i][j][7]=0 とすることができ、制約 \((c)\) を課したことになります。

さらに、数字が確定しているマスにその数字が一つだけ入るための制約 \((d)\) を課します。上記の例では、q[1][5][k]=0 ( \(k\neq7\) )となります。初期配置から与えらる全てのマスについて同様な操作を行うことで必要な変数を絞り込み、より少ない変数で計算を行うことができるようにします。

これで初期設定ができました。例として行番号 \(0\)、列番号 \(0\) の9変数を表示すると2番目の要素が1として確定、つまり \(2\) が確定していることが確認できます。

>>> q[0][0]
[0, 1, 0, 0, 0, 0, 0, 0, 0]

同様に行番号 \(0\)、列番号 \(1\) の9変数を表示すると、行・列・ブロック内に表れる数字 \(1,2,3,4,5,6\) が候補から外れていることが確認できます。

>>> q[0][1]
[0, 0, 0, 0, 0, 0, q_15, q_16, q_17]

制約条件の設定

次に、制約条件からコスト関数を定義します。先ほどの \((a)\sim(d)\) の one-hot 制約条件は Amplify の equal_to 関数を用いて表すことができます。まず、\((a)\) の各行には同じ数字が入らないという制約を表すコスト関数を定義してみます。行 i とレイヤ k で指定される全ての列に対する変数の和は \(\sum_{j=0}^{8}q[i][j][k]\) で与えられるので、これが \(1\) になる制約は次のように表されます。

同様にして、\((b)\) の同じ列に同じ数字が入らない制約条件と、\((d)\) の一つのマスには一つの数字しか入らない制約条件は以下のように表されます。

最後に \((c)\)\(3\times3\) の各ブロックには同じ数字が入らないという制約条件を表します。全てのレイヤに対して各 \(3\times3\) ブロック内で変数の和を取り、以下のようにして one-hot 制約を課します。

これで全ての制約条件が出そろったので、これらの制約条件を全て足し合わせます。

これで定式化に関する準備ができました。イジングマシンによって、全ての制約を満たす変数の組み合わせを見つけ出すことが出来れば、そのような変数の組み合わせは与えられた初期配置から導き出される数独の解となります。

イジングマシンの実行

先ほど作成した constraints を用いてイジングマシンを実行する準備を行います。まずイジングマシンのクライアントを作成し、timeout などのパラメーターを設定します。その後、ソルバーを作成しクライアントを設定します。

制約条件 constraintsBinaryQuadraticModel に与えることでの論理模型クラスとして定式化し、これを先ほど設定した solver に与えて実行します。

実行結果は values に格納されています。 decode_solution 関数に values と変数配列 q を与えることで変数配列に結果が代入されます。その後、全ての i, j に対して q[i][j][k] = 1 となる k を検索することで、それぞれの行と列における k + 1 を数独の解として取得できます。

最後に print_sudoku 関数を用いて解答を出力します。

>>> print_sudoku(answer)
2 8 5 | 1 3 9 | 6 7 4
6 7 3 | 2 4 8 | 5 1 9
4 1 9 | 6 5 7 | 3 2 8
---------------------
7 3 8 | 5 6 4 | 1 9 2
5 4 2 | 3 9 1 | 7 8 6
1 9 6 | 7 8 2 | 4 5 3
---------------------
8 6 1 | 9 7 3 | 2 4 5
9 5 7 | 4 2 6 | 8 3 1
3 2 4 | 8 1 5 | 9 6 7

数独の一般化

これまでは \(3 \times 3\) のブロックに区切られた \(9 \times 9\) マスの数独を取り扱いましたが、問題サイズを拡張した \(16\times16\)\(25\times25\) 等の数独にもイジングマシンは容易に対応できます。数独のマスの数を \(N\times N,\, (N\in\mathbb{Z})\)、区切られたブロックを \(n \times n\) (ただし \(N=n^2,\,(n\in\mathbb{Z})\)) とします。例えば、基本的な \(9\times9\) の数独の場合は、\(N=9\)\(n=3\) となります。

先ほどのコードを \(N\)\(n\) を用いて一般化し、\(16\times16\) マスの数独を例として解いてみます。

from amplify import (
    BinaryPoly,
    BinaryQuadraticModel,
    gen_symbols,
    Solver,
    decode_solution,
)
from amplify.constraint import equal_to
from amplify.client import FixstarsClient
import numpy as np

n = 4  # ブロックサイズ
N = n * n  # 全体のサイズ

# n = 4 の初期値
# 引用元: https://www.free-sudoku-puzzle.com/puzzle_fours/solve/3/238
initial = np.array(
    [
        [0, 0, 0, 7, 0, 0, 0, 1, 5, 0, 3, 16, 4, 0, 15, 0],
        [0, 11, 0, 0, 0, 0, 5, 0, 0, 2, 12, 6, 0, 0, 7, 14],
        [4, 0, 0, 0, 7, 8, 9, 0, 11, 0, 1, 15, 0, 0, 10, 0],
        [10, 0, 0, 0, 0, 0, 0, 15, 13, 0, 9, 7, 8, 0, 0, 1],
        [13, 0, 0, 16, 15, 0, 4, 9, 0, 0, 14, 0, 11, 0, 1, 0],
        [8, 0, 5, 0, 0, 0, 10, 0, 0, 0, 0, 0, 15, 0, 0, 0],
        [0, 0, 0, 11, 0, 0, 0, 8, 16, 7, 0, 9, 0, 0, 0, 0],
        [0, 0, 0, 14, 0, 0, 3, 0, 4, 0, 0, 5, 13, 0, 0, 0],
        [0, 0, 0, 0, 3, 0, 0, 14, 0, 0, 4, 0, 9, 12, 8, 15],
        [0, 0, 0, 0, 0, 1, 7, 10, 0, 15, 8, 11, 0, 0, 0, 0],
        [0, 0, 0, 0, 11, 12, 0, 0, 0, 0, 16, 0, 3, 5, 0, 0],
        [0, 0, 0, 0, 0, 16, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 6, 16, 0, 15, 1, 5, 0, 14, 2, 0, 0],
        [0, 0, 0, 3, 0, 0, 0, 0, 9, 0, 0, 14, 0, 1, 0, 4],
        [2, 0, 12, 0, 0, 0, 0, 0, 0, 16, 13, 0, 6, 0, 3, 5],
        [1, 0, 0, 0, 0, 15, 0, 0, 2, 11, 6, 12, 7, 9, 0, 10],
    ]
)

# 数独を成形して表示する
def print_sudoku(sudoku):
    print("\n\n")
    width = len(str(N))
    for i in range(len(sudoku)):
        line = ""
        if i % n == 0 and i != 0:
            print("-" * ((width + 1) * n * n + 2 * (n - 1)))
        for j in range(len(sudoku[i])):
            if j % n == 0 and j != 0:
                line += "| "
            line += str(sudoku[i][j]).rjust(width) + " "
        print(line)


q = gen_symbols(BinaryPoly, N, N, N)

for i, j in zip(*np.where(initial != 0)):
    k = initial[i, j] - 1

    for m in range(N):
        q[i][m][k] = BinaryPoly(0)  # 制約(a)
        q[m][j][k] = BinaryPoly(0)  # 制約(b)
        q[i][j][m] = BinaryPoly(0)  # 制約(d)
        q[n * (i // n) + m // n][n * (j // n) + m % n][k] = BinaryPoly(0)  # 制約(c)

    q[i][j][k] = BinaryPoly(1)  # 変数の値を1に指定する

# (a): 各行には同じ数字が入らない制約条件
row_constraints = [
    equal_to(sum(q[i][j][k] for j in range(N)), 1) for i in range(N) for k in range(N)
]

# (b): 各列には同じ数字が入らない制約条件
col_constraints = [
    equal_to(sum(q[i][j][k] for i in range(N)), 1) for j in range(N) for k in range(N)
]

# (d): 一つのマスには一つの数字しか入らない制約条件
num_constraints = [
    equal_to(sum(q[i][j][k] for k in range(N)), 1) for i in range(N) for j in range(N)
]

# (c): nxnブロック内には同じ数字が入らない制約条件
block_constraints = [
    equal_to(sum([q[i + m // n][j + m % n][k] for m in range(N)]), 1)
    for i in range(0, N, n)
    for j in range(0, N, n)
    for k in range(N)
]

constraints = (
    sum(row_constraints)
    + sum(col_constraints)
    + sum(num_constraints)
    + sum(block_constraints)
)

client = FixstarsClient()
client.url = "http://xxx.xxx.xxx.xxx"
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = 10000  # タイムアウト10秒

solver = Solver(client)
model = BinaryQuadraticModel(constraints)
result = solver.solve(model)
if len(result.solutions) == 0:
    raise RuntimeError("Any one of constaraints is not satisfied.")

values = result.solutions[0].values

q_values = decode_solution(q, values)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(N)])
>>> print_sudoku(initial)
 0  0  0  7 |  0  0  0  1 |  5  0  3 16 |  4  0 15  0
 0 11  0  0 |  0  0  5  0 |  0  2 12  6 |  0  0  7 14
 4  0  0  0 |  7  8  9  0 | 11  0  1 15 |  0  0 10  0
10  0  0  0 |  0  0  0 15 | 13  0  9  7 |  8  0  0  1
------------------------------------------------------
13  0  0 16 | 15  0  4  9 |  0  0 14  0 | 11  0  1  0
 8  0  5  0 |  0  0 10  0 |  0  0  0  0 | 15  0  0  0
 0  0  0 11 |  0  0  0  8 | 16  7  0  9 |  0  0  0  0
 0  0  0 14 |  0  0  3  0 |  4  0  0  5 | 13  0  0  0
------------------------------------------------------
 0  0  0  0 |  3  0  0 14 |  0  0  4  0 |  9 12  8 15
 0  0  0  0 |  0  1  7 10 |  0 15  8 11 |  0  0  0  0
 0  0  0  0 | 11 12  0  0 |  0  0 16  0 |  3  5  0  0
 0  0  0  0 |  0 16  0  0 |  0  5  0  0 |  0  0  0  0
------------------------------------------------------
 0  0  0  0 |  0  6 16  0 | 15  1  5  0 | 14  2  0  0
 0  0  0  3 |  0  0  0  0 |  9  0  0 14 |  0  1  0  4
 2  0 12  0 |  0  0  0  0 |  0 16 13  0 |  6  0  3  5
 1  0  0  0 |  0 15  0  0 |  2 11  6 12 |  7  9  0 10
>>> print_sudoku(answer)
12  2  8  7 | 14 11  6  1 |  5 10  3 16 |  4 13 15  9
 9 11  1 15 | 10  4  5 13 |  8  2 12  6 | 16  3  7 14
 4  3 16 13 |  7  8  9  2 | 11 14  1 15 |  5  6 10 12
10  5 14  6 | 16  3 12 15 | 13  4  9  7 |  8 11  2  1
------------------------------------------------------
13 12 10 16 | 15  5  4  9 |  6  3 14  8 | 11  7  1  2
 8  7  5  9 |  6 14 10 16 |  1 13 11  2 | 15  4 12  3
 3  4  2 11 | 12 13  1  8 | 16  7 15  9 | 10 14  5  6
 6  1 15 14 |  2  7  3 11 |  4 12 10  5 | 13  8  9 16
------------------------------------------------------
 5 16 11  1 |  3  2 13 14 |  7  6  4 10 |  9 12  8 15
14  6  9 12 |  5  1  7 10 |  3 15  8 11 |  2 16  4 13
15 10 13  2 | 11 12  8  4 | 14  9 16  1 |  3  5  6  7
 7  8  3  4 |  9 16 15  6 | 12  5  2 13 |  1 10 14 11
------------------------------------------------------
11  9  7 10 |  4  6 16 12 | 15  1  5  3 | 14  2 13  8
16 15  6  3 | 13 10  2  5 |  9  8  7 14 | 12  1 11  4
 2 14 12  8 |  1  9 11  7 | 10 16 13  4 |  6 15  3  5
 1 13  4  5 |  8 15 14  3 |  2 11  6 12 |  7  9 16 10