はじめに
これは、前回(optunaに数独を解かせてみた)の続編です。
前回は数独のルールすら知らない状態で、optunaを使って最適化することで正解にたどり着けるか試してみました。
今回は、縦・横・3x3のブロックで数字の重複がないというルールを知ったうえでoptunaに最適化させてみました。
もっと早くに投稿するつもりだったのですが、前回の投稿後に始めたリングフィット アドベンチャーがハードすぎて気力・体力がともになくなって、遅くなってしまいました。
方法
今回は、焼き鈍し法で解きます。
焼き鈍し法は、現在の値を少しランダムに変化させて次の値を作ります。(初回は元の値がないので、完全にランダムに作ります。)
作った新しい値から計算されるコストが減少するなら新しい値に移動します。
減少しない場合でも確率で移動したり、しなかったりします。
徐々にその確率を小さくしていくことで、局所解に収束します。
実装
optunaで焼き鈍しをする方法は公式ドキュメントにあります。
新しい値をサンプリングする部分を変更しただけで、後は公式ドキュメントのとおりです。
サンプラーを変更した以外は前回のコードのままです。
ここでは、新しい値をサンプリングする部分のみを示します。
コード全体は最後に示します。
params = {}
# サンプリングの順序をランダムにする
name_list = shuffle(list(search_space))
for param_name in name_list:
# 2次元の座標(i, j)
i, j = int(param_name[1]), int(param_name[2])
# 予め(i, j)がある縦・横・3x3のブロックの要素を取り出すマスクを作っておいた
# ただし、(i, j)の要素はマスクでは取り出さない
mask = self._mask_array[i, j]
tmp = self._state * mask
# 1~9の各要素がそれぞれいくつあるかを数える
cost = np.asarray([np.count_nonzero(tmp == v) for v in range(1, 10)])
probability = softmax(-cost * 5)
# 新しい値をサンプリング
new_value = np.random.choice(9, p=probability) + 1
# 新しい値を記録
params[param_name] = new_value
self._state[i, j] = new_value
return params
新しい値は、縦・横・3x3のブロックで数字の重複がないというルールに基づいてサンプリングされる確率が高くなるようにしています。
ここでルールについての知識を利用しています。
実験結果
左に*がある数字は元からヒントとして与えられているものです。
*5*3 2 4*7 6 9 1 8
*6 7 4*1*9*5 3 5 2
1*9*8 2 3 8 4*6 7
*8 1 9 5*6 4 7 2*3
*4 2 6*8 7*3 5 9*1
*7 5 3 9*2 1 8 4*6
9*6 1 3 5 7*2*8 4
2 8 7*4*1*9 6 3*5
3 4 5 6*8 2 1*7*9
loss: 4.0
惜しい。
トライアル数が115回目にこの答えに辿り着きました。
この後も、200回まで行いましたがダメでした。
一目ではどこが間違っているかわからないですが、中央で縦方向に7が重複しています。
中央の7が縦方向では間違っていますが、横方向では条件を満たしています。
そのため、複数の要素が同時に変わらないと正解に辿り着けないので、かなり深い局所解になっています。(コストがかなり高い状態を経由しないと正解に辿り着かない)
まとめ
局所解になってしまいました。
正解にたどり着くためには、複数の要素が同時に変わらないといけないので、現在の方法では難しいようです。
初期値を変えて何度もやり直せば、正解できるはずです。
今後の予定
今回は中央で縦方向に7が重複していましたが、7は元から与えられているものがあるので、サンプリングされるべきでないはずです。
(ルールに関する知識をどこまで利用するかという条件に関係してくる)
また、右から3列目、上から3行目の様に元から与えられている数字から値が確定するものもあります。
これらが今回の方法では全く考慮されていなかったので、次回はこれを解決したいと思います。
(もうoputunaは関係ない気がする)
コード
from itertools import product
import numpy as np
import optuna
import click
from scipy.special import softmax
from sklearn.utils import shuffle
# https://optuna.readthedocs.io/en/stable/tutorial/sampler.html
class SimulatedAnnealingSampler(optuna.samplers.BaseSampler):
def __init__(self, temperature=100):
self._rng = np.random.RandomState()
self._temperature = temperature
self._current_trial = None
self._state = None
self._mask_array = make_mask_array()
def sample_relative(self, study, trial, search_space):
if search_space == {}:
return {}
# 現在のtrialはstudy.trials[-1]
previous_trial = study.trials[-2]
if self._current_trial is None or previous_trial.value <= self._current_trial.value:
probability = 1.0
else:
probability = np.exp((self._current_trial.value - previous_trial.value) / self._temperature)
self._temperature *= 0.99
if self._rng.uniform(0, 1) < probability:
self._current_trial = previous_trial
if self._state is None:
self._state = np.empty([9, 9], dtype=np.int32)
for i, j in product(range(9), repeat=2):
name = 'p{}{}'.format(i, j)
if name in preset:
self._state[i, j] = preset[name]
for i, j in product(range(9), repeat=2):
name = 'p{}{}'.format(i, j)
if name in self._current_trial.params:
self._state[i, j] = self._current_trial.params[name]
params = {}
name_list = shuffle(list(search_space))
for param_name in name_list:
i, j = int(param_name[1]), int(param_name[2])
mask = self._mask_array[i, j]
tmp = self._state * mask
cost = np.asarray([np.count_nonzero(tmp == v) for v in range(1, 10)])
probability = softmax(-cost * 5)
new_value = np.random.choice(9, p=probability) + 1
params[param_name] = new_value
self._state[i, j] = new_value
return params
def infer_relative_search_space(self, study, trial):
return optuna.samplers.intersection_search_space(study)
def sample_independent(self, study, trial, param_name, param_distribution):
independent_sampler = optuna.samplers.RandomSampler()
return independent_sampler.sample_independent(study, trial, param_name, param_distribution)
def make_mask_array():
mask_array = np.zeros([9, 9, 9, 9])
for i, j in product(range(9), repeat=2):
mask = mask_array[i, j]
mask[i] = 1
mask[:, j] = 1
s, t = i // 3 * 3, j // 3 * 3
mask[s:s + 3, t:t + 3] = 1
# 自分自身を取り除く
mask[i, j] = 0
return mask_array
"""
-----------------
|5|3| | |7| | | | |
|-+-+-+-+-+-+-+-+-|
|6| | |1|9|5| | | |
|-+-+-+-+-+-+-+-+-|
| |9|8| | | | |6| |
|-+-+-+-+-+-+-+-+-|
|8| | | |6| | | |3|
|-+-+-+-+-+-+-+-+-|
|4| | |8| |3| | |1|
|-+-+-+-+-+-+-+-+-|
|7| | | |2| | | |6|
|-+-+-+-+-+-+-+-+-|
| |6| | | | |2|8| |
|-+-+-+-+-+-+-+-+-|
| | | |4|1|9| | |5|
|-+-+-+-+-+-+-+-+-|
| | | | |8| | |7|9|
-----------------
"""
preset = {'p00': 5, 'p01': 3, 'p04': 7,
'p10': 6, 'p13': 1, 'p14': 9, 'p15': 5,
'p21': 9, 'p22': 8, 'p27': 6,
'p30': 8, 'p34': 6, 'p38': 3,
'p40': 4, 'p43': 8, 'p45': 3, 'p48': 1,
'p50': 7, 'p54': 2, 'p58': 6,
'p61': 6, 'p66': 2, 'p67': 8,
'p73': 4, 'p74': 1, 'p75': 9, 'p78': 5,
'p84': 8, 'p87': 7, 'p88': 9}
def evaluate(answer):
tmp = np.reshape(answer, [3, 3, 3, 3])
loss = np.sum((
np.sum([np.count_nonzero(np.logical_not(np.any(answer == i, axis=0))) for i in range(1, 10)]),
np.sum([np.count_nonzero(np.logical_not(np.any(answer == i, axis=1))) for i in range(1, 10)]),
np.sum([np.count_nonzero(np.logical_not(np.any(tmp == i, axis=(1, 3)))) for i in range(1, 10)]),
))
return loss
def objective(trial):
candidate = (1, 2, 3, 4, 5, 6, 7, 8, 9)
answer = np.empty([9, 9], dtype=np.uint8)
for i, j in product(range(9), repeat=2):
key = 'p{}{}'.format(i, j)
if key in preset:
answer[i, j] = preset[key]
else:
answer[i, j] = trial.suggest_categorical(key, candidate)
return evaluate(answer)
def run(n_trials):
study_name = 'sudoku'
sampler = SimulatedAnnealingSampler()
study = optuna.create_study(study_name=study_name, storage='sqlite:///sudoku.db', load_if_exists=True,
sampler=sampler)
study.optimize(objective, n_trials=n_trials)
show_result(study.best_params, study.best_value)
df = study.trials_dataframe()
df.to_csv('tpe_result.csv')
def show_result(best_params, best_value):
for i in range(9):
for j in range(9):
key = 'p{}{}'.format(i, j)
if key in preset:
print('*{:1d}'.format(preset[key]), end='')
else:
print('{:2d}'.format(best_params[key]), end='')
print('')
print('loss: {}'.format(best_value))
@click.command()
@click.option('--n-trials', type=int, default=1000)
def cmd(n_trials):
run(n_trials)
def main():
cmd()
if __name__ == '__main__':
main()