当社はCookieを使用して、お客様が当社のWebサイトでより良い体験を得られるようにしています。引き続き閲覧する場合は、当社グループのプライバシーポリシー 
に同意したことになります。
このチュートリアルでは、Amplifyを用いたイジングマシンによる数独の解法について解説します。
数独（すうどく）は、以下のルールに従って$9\times9$のブロックに$1\sim9$の数字を入れるパズルです。
まず、ヒントとして17個以上の数字が埋められた初期配置が与えられます (ヒントが16個以下の初期配置は解法を持ちえないことが証明されています)。上記のルールに従って、空いているマスに入る数字を確定していくことでゲームを進めることができます。ゲームの難易度が低い場合は、比較的簡単に次々と入り得る数字を確定できるマスを見つけることができますが、ゲームの難易度が上がると、そのような数字の確定が難しくなり、ある程度の経験を積まないとパズルを解き進めるのは難しくなります。
コンピュータを使って数独を解く方法として、深さ優先探索、確率的な方法、制約充足問題、exact cover problemなど用いた様々なアルゴリズムが考案されています。これらのアルゴリズムに従って、機械的に数独を解くことができます。
このチュートリアルでは、組み合わせ最適化問題に特化したイジングマシンを使って数独を解く方法を紹介します。上記の数独のルールを制約条件として解釈し、それに伴ったコスト関数を定義し、コストが最も低い数字の組み合わせを見つけることで数独の解を見つけることができる仕組みになっています。よって、ここで行うべきことは数独のルールをどのような制約条件によって表すことができるかを考えることです。そのような制約条件を見つけてコスト関数を定義することができれば、あとは初期配置を与えるだけで、複雑なアルゴリズムを用いることなく、イジングマシンによって解が見つけることができます。
それでは次に、Amplifyを使って、数独を解くコードがどのように実装されるのか具体的に見てみましょう。
次に、イジングマシンにおいて、数独のルールを表す制約条件を満たすコスト関数の作成方法を考えましょう。基本的には、二次二値多変数多項式を用いて制約条件を表現する方法を考えることとなります。ここでは、QUBO模型（各変数は0または1）を用いた一つの方法に着目して議論を進めます。
表すべき数独のルールは以下の3つとなります。
まず、$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$ 行 $5$ 列のマスに $7$ が入る状態は $q_{2,4,6}=1$、そうでない場合は $q_{2,4,6}=0$ と表現することができます。これらの変数を使うと、制約条件1～3はそれぞれ以下のone-hot制約として書き下すことができます。
$$ \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. $$制約条件$(a)$、$(b)$, $(c)$はそれぞれルール1、2、3に対応します。$(d)$の制約条件は各マスには一つの数字しか入らないという基本的な条件に対応します。
数独では、$17$個以上のいくつかのマスがすでに埋められた初期配置がヒントとして与えられます。ここでは、難しい問題とされる以下の初期配置を使います。以下の表記では、数字が埋められていないマスは$0$としました。
import numpy as np
# 初期配置をリストで表記
# 引用元: 数独問題集-2018年3月1日-最高級-( http://www.sudokugame.org/archive/printable.php?nd=4&y=2018&m=03&d=1 )
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):
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 BinarySymbolGenerator
gen = BinarySymbolGenerator()
q = gen.array(9, 9, 9)
これによって、$9^3=729$ 個の変数が三次元配列として用意されました。9, 9, 9
はそれぞれ、行・列・数値レイヤの要素数を表し、それらのインデックスをi、j、kとして各要素にはq[i, j, k]でアクセスできます。例えば
行番号 $0$、列番号 $0$ の9変数を表示する場合は次のようにします。
# 行番号 $0$、列番号 $0$ の9変数
print(q[0, 0])
# 行番号 $2$、数値レイヤ $5$ の9変数
print(q[2, :, 5])
先ほど initial
に格納された初期配置から数独のルールに従い、確定可能な未知変数の絞り込みを行います。例えば、i=1、j=5
のマスにはすでに 8 (initial[1][5]=8) が入っているので、8 に対応した k=7
のレイヤの変数は q[1][5][7]=1 と指定されます。
ルール1、2により、$k=7$ のレイヤで i=1 と j=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$)
となります。初期配置から与えらる全てのマスについて同様な操作を行うことで必要な変数を絞り込み、より少ない変数で計算を行うことができるようにします。
for i, j in zip(*np.where(initial != 0)): # 0 ではない行と列とレイヤのインデックスを取得
k = initial[i, j] - 1 # 値からインデックスに変換するため -1 することに注意
q[i, :, k] = 0 # 制約(a)
q[:, j, k] = 0 # 制約(b)
q[i, j, :] = 0 # 制約(d)
for m in range(9):
q[(3 * (i // 3) + m // 3), (3 * (j // 3) + m % 3), k] = 0 # 制約(c)
q[i, j, k] = 1 # 変数の値を1に確定させる
これで初期設定ができました。例として行番号 $0$、列番号 $0$ の9変数を表示すると2番目の要素が1として確定、つまり $2$ が確定していることが確認できます。
q[0, 0]
同様に行番号 $0$、列番号 $1$ の9変数を表示すると、行・列・ブロック内に表れる数字 $1,2,3,4,5,6$ が候補から外れていることが確認できます。
q[0, 1]
次に、制約条件からコスト関数を定義します。先ほどの $(a)\sim(d)$ の one-hot 制約条件は Amplify のone_hot
関数を用いて表すことができます。まず、$(a)$ の各行には同じ数字が入らないという制約を表すコスト関数を定義してみます。行 i とレイヤ k
で指定される全ての列に対する変数の和は $\sum_{j=0}^{8}q[i][j][k]$ で与えられるので、これが $1$ になる制約は次のように表されます。
from amplify import sum_poly
from amplify.constraint import one_hot
# (a): 各行には同じ数字が入らない制約条件
row_constraints = [one_hot(q[i, :, k]) for i in range(9) for k in range(9)]
同様にして、$(b)$ の同じ列に同じ数字が入らない制約条件と、$(d)$ の一つのマスには一つの数字しか入らない制約条件は以下のように表されます。
# (b): 各列には同じ数字が入らない制約条件
col_constraints = [one_hot(q[:, j, k]) for j in range(9) for k in range(9)]
# (d): 一つのマスには一つの数字しか入らない制約条件
num_constraints = [one_hot(q[i, j, :]) for i in range(9) for j in range(9)]
最後に$(c)$ の $3\times3$ の各ブロックには同じ数字が入らないという制約条件を表します。全てのレイヤに対して各 $3\times3$ ブロック内で変数の和を取り、以下のようにして one-hot 制約を課します。
# (c): 3x3ブロック内には同じ数字が入らない制約条件
block_constraints = [
one_hot(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 = (
sum(row_constraints)
+ sum(col_constraints)
+ sum(num_constraints)
+ sum(block_constraints)
)
これで定式化に関する準備ができました。イジングマシンによって、全ての制約を満たす変数の組み合わせを見つけ出すことが出来れば、そのような変数の組み合わせは与えられた初期配置から導き出される数独の解となります。
先ほど作成した constraints を用いてイジングマシンを実行する準備を行います。まずイジングマシンのクライアントを作成し、timeout
などのパラメーターを設定します。その後、ソルバーを作成しクライアントを設定します。
from amplify import Solver
from amplify.client import FixstarsClient
client = FixstarsClient()
client.parameters.timeout = 1000 # タイムアウト1秒
solver = Solver(client)
制約条件 constraints を BinaryQuadraticModel に与えることでの論理模型クラスとして定式化し、これを先ほど設定した
solver に与えて実行します。
from amplify import BinaryQuadraticModel
model = BinaryQuadraticModel(constraints)
result = solver.solve(model)
if len(result) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
values = result[0].values
実行結果は values に格納されています。 変数配列 qのdecode メンバ関数に values
を与えることで変数配列に結果が代入されます。その後、全ての i, j に対して q[i, j, k] = 1 となる
k を検索することで、それぞれの行と列における k + 1 を数独の解として取得できます。
from amplify import decode_solution
q_values = q.decode(values)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(9)])
最後に print_sudoku 関数を用いて解答を出力します。
print_sudoku(answer)
これまでは $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 BinarySymbolGenerator, BinaryQuadraticModel, Solver
from amplify.constraint import one_hot
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):
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 = BinarySymbolGenerator().array(N, N, N)
for i, j in zip(*np.where(initial != 0)):
k = initial[i, j] - 1
q[i, :, k] = 0 # 制約(a)
q[:, j, k] = 0 # 制約(b)
q[i, j, :] = 0 # 制約(d)
for m in range(N):
q[(n * (i // n) + m // n), (n * (j // n) + m % n), k] = 0 # 制約(c)
q[i, j, k] = 1 # 変数の値を1に指定する
# (a): 各行には同じ数字が入らない制約条件
row_constraints = [one_hot(q[i, :, k]) for i in range(N) for k in range(N)]
# (b): 各列には同じ数字が入らない制約条件
col_constraints = [one_hot(q[:, j, k]) for j in range(N) for k in range(N)]
# (d): 一つのマスには一つの数字しか入らない制約条件
num_constraints = [one_hot(q[i, j, :]) for i in range(N) for j in range(N)]
# (c): nxnブロック内には同じ数字が入らない制約条件
block_constraints = [
one_hot(sum([q[i + m // n, j + m % n, k] for m in range(N)]))
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.parameters.timeout = 10000 # タイムアウト10秒
solver = Solver(client)
model = BinaryQuadraticModel(constraints)
result = solver.solve(model)
if len(result) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
values = result[0].values
q_values = q.decode(values)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(N)])
print_sudoku(initial)
print_sudoku(answer)
