はじめに
ここでは、数独を、Julia上でSA（シミュレーティッド ・アニーリング：焼きなまし法）を使って解いてみたいと思います。
なお、数独は、既存のソルバーで解くことができ、様々な開放が紹介されています。例えば、JuMPで解く方法はこちら（Sudoku · JuMP）に掲載されています。そしてSAで解くよりも高速です。
ここでは、SAを使って数独を解くことができるか、という観点からの実験結果です。
実行環境について
以下では、「SageMaker Studio Lab」上で実行した結果を記しています。
「SageMaker Studio Lab」を使う方法については、別記事『JuliaをSageMaker Studio Labで使う』をご覧ください。
また、同じコードがWindows10上のJupyterLabでも動作することを確認しています。
数独について
数独（すうどく）とは、9×9の正方形の枠内を1〜9までの数字で埋めるパズルです。
次の条件に合うように数字を当てはめていき、すべてのマスを埋めることを目指します。
- いくつかのマスには、初期値としての数値が入っていて、これは変更できない。
- 使える数字は、1~9
- 縦・横の各列には、同じ数字が重複してはいけない。
- 9×9の枠は、3×3の枠（ブロック）が9個で構成されていて、ブロック内にも同じ数字が重複してはいけない。
詳細は、下記を参照してください。
SA法（シミュレーティド・アニーリング法、焼きなまし法）とは
解空間を探索する際には、最適解に近づいたかどうかを局所的に判定しながら進める方法が一般的なのですが、SAでは必ずしも最適解に近づく変化だけを行うわけではありません。ある確率で、最適解から遠ざかる変化も許容します。これは、最適解に到達する前に局所解に陥ってしまうのを防ぐ効果があります。
詳細は、下記を参照してください。
数独のエネルギー関数
数独の制約として、定められた範囲（横、縦、3×3領域）で1~9の数値が必ず１つずつ出現します。これを素直にエネルギー関数にしようとすると、かなり複雑になってしまいます。（不可能ではないが現実的ではない）
そこで今回は、『Sudoku · JuMP』のページで紹介されているように、1~9の数字を9個の要素を持つベクトルとして定義します。このベクトルの要素は0か1の値のみを持ち、1の値を持つインデックス（1~9）が、本来の数値を表すようにします。
このままでは、「数字を表現するベクトルに同じ数値は一つだけ」とするエネルギー関数を記述する必要があります。そこで、状態変化の方法を工夫することでワンホット表現（要素のうち一つだけが1で他はすべて0）であることを維持します。これによって「数字を表現するベクトルに同じ数値は一つだけ」という状態が常に保証されるため、関連するエネルギー関数は不要となります。
今回は次のような計算式を採用しました。この$H$の値が最小になる状態をSAで求めます。そして最小値0になった場合に数独の解が得られたことになります。
この計算式（$H$）では、同じ数字が、縦には1個、横には1個、そして3×3のブロックにも1個、という条件をすべて満たした場合に、最小値0になります。
$H = H_{row} + H_{col} + H_{block}$
$H_{row} = \displaystyle \sum_{num=1}^{9} (\displaystyle \sum_{row} (\displaystyle \sum_{i∈row} x_{i, num}-1)^{2} )$
$H_{col} = \displaystyle \sum_{num=1}^{9} (\displaystyle \sum_{col} (\displaystyle \sum_{i∈col} x_{i, num}-1)^{2} )$
$H_{block} = \displaystyle \sum_{num=1}^{9} (\displaystyle \sum_{block} (\displaystyle \sum_{i∈block} x_{i, num}-1)^{2} )$
プログラム
枠全体の挙動を定義します。
# Sudoku by Simulated Annealing
# シミュレーティッドアニーリング（焼きなまし法）による数独
struct Sudoku
# 枠全体：縦x横x数字
unit::Array{Int8, 3}
# 固定位置
fixed::Vector{Tuple{Int64,Int64}}
function Sudoku(fixed::Matrix{Int64})
# 初期化
u = zeros(Int8, 9, 9, 9)
f = []
for x = 1:9
for y = 1:9
if fixed[x, y] != 0
push!(f, (x, y))
# 固定値が指定されていればそれを設定
z = fixed[x, y]
u[x, y, z] = 1
else
# 指定がなければ乱数を設定
z = rand(1:9)
u[x, y, z] = 1
end
end
end
new(u, f)
end
end
function energy(board::Sudoku)
#
# エネルギー計算
# 解の場合に、エネルギー=0となる
#
e = 0
# row
e += sum([(sum([board.unit[x, y, z] for y = 1:9]) - 1)^2 for x = 1:9 for z = 1:9])
# col
e += sum([(sum([board.unit[x, y, z] for x = 1:9]) - 1)^2 for y = 1:9 for z = 1:9])
# block
e += sum([(sum([board.unit[x+i, y+j, z] for i = 0:2 for j = 0:2]) - 1)^2 for x = 1:3:7 for y = 1:3:7 for z = 1:9])
return e
end
function change(board::Sudoku, x, y, r)
# 指定されたrの分だけ要素をrotateする → 数字変化
# 常に1の値は一つのみで、ワンホット状態が維持される
for z1 = 1:9
if board.unit[x, y, z1] == 1
z2 = (z1 - 1 + r) % 9 + 1
board.unit[x, y, z1], board.unit[x, y, z2] = 0, 1
break
end
end
end
function display(board::Sudoku)
println()
for x = 1:9
s = []
for y = 1:9
for z = 1:9
if board.unit[x, y, z] == 1
push!(s, z)
break
end
end
end
println("[" * join(s, ", ") * "]")
end
end
function run(board::Sudoku, temp_start, temp_end, alpha, temp_iter)
now_energy = energy(board)
temperature = temp_start
while temperature > temp_end && now_energy > 0
for i = 1:temp_iter
# 一か所変化させる
cx = rand(1:9)
cy = rand(1:9)
# 固定値のマスが指定されたらやり直し
while (cx, cy) in board.fixed
cx = rand(1:9)
cy = rand(1:9)
end
r = rand(1:8)
change(board, cx, cy, r)
# エネルギー値の計算
new_energy = energy(board)
# 変化前とのエネルギー値の差分
de = new_energy - now_energy
# eval
if exp(- de / temperature) > rand() # 条件 (de < 0) 含む
# 変化を許す
now_energy = new_energy
#display(board)
#println("energy = ", now_energy, "temp = ", temperature)
# 解に到達したら終了
if now_energy == 0
return (true, now_energy, temperature)
end
else
# 変化を許さないので元に戻す
change(board, cx, cy, 9-r)
end
end
#display(board)
#println("energy = ", now_energy, " temp = ", temperature)
# 温度を減衰
temperature = temperature * alpha
end
return (false, now_energy, temperature)
end
実行例
固定値、および、各パラメータを指定して実行
fixed = [
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
]
# 開始温度
temp_start = 1000
# 減衰率
alpha = 0.99
# 終了温度
temp_end = 0.02
# 温度当たりの繰り返し数
temp_iter = 1000
board = Sudoku(fixed)
(result, ene, temp), runtime = @timed run(board, temp_start, temp_end, alpha, temp_iter)
display(board)
println("energy = ", ene, ", temp = ", temp)
if result
println("OK")
end
println(runtime, " sec") 実行すると、実行後の枠の状態、エネルギー値、温度、実行時間が表示されます。
上記を実行した場合、次のような結果が得られました。
[5, 3, 4, 6, 7, 8, 9, 1, 2]
[6, 7, 2, 1, 9, 5, 3, 4, 8]
[1, 9, 8, 3, 4, 2, 5, 6, 7]
[8, 5, 9, 7, 6, 1, 4, 2, 3]
[4, 2, 6, 8, 5, 3, 7, 9, 1]
[7, 1, 3, 9, 2, 4, 8, 5, 6]
[9, 6, 1, 5, 3, 7, 2, 8, 4]
[2, 8, 7, 4, 1, 9, 6, 3, 5]
[3, 4, 5, 2, 8, 6, 1, 7, 9]
energy = 0, temp = 0.7200126227424993
OK
25.856435285 sec この例では、最終的にエネルギーが0になって、最適解が得られていることがわかります。しかし、SA法では近似解を求めているので、必ずしも最適解にならない場合も出てきます。
その場合は、何度か繰り返すことで、「energy = 0」の状態を得ることができます。
備考
上記の実行例は、簡単な例だったので、何度かの繰り返しで、解を得ることができました。しかし、より難度の高い問題については、このパラメータ設定では、解を得るのは難しいようです。
推定される原因としては、人間にとって難度の高い問題は、局所的にエネルギーが低くなる状態（いわゆる局所解）が多く存在するのではないかと推定しています。
そうであるならば、各種パラメータを調整して、よりゆっくりとした焼きなましを行うことで、解を得ることができるのではないかと考えられます。
そのうえで、今回は、上記の問題に対して、解が得られるようにパラメータを調整しました。
まとめ
上記を実現したプログラムを、下記に格納してあります。
コメント