数独¶
このページでは、Hi-QUBO を使って数独パズルをイジングマシン（アニーリング）によって解く方法を紹介します。
数独とは¶
数独（Sudoku）は、9×9 グリッドに 1 〜 9 の数字を埋めるパズルです。
次のルールを満たすように空欄に数字を入れていきます。
各行には 1 〜 9 の数字が一度ずつ現れること
各列には 1 〜 9 の数字が一度ずつ現れること
3×3 のブロック内にも 1 〜 9 の数字が一度ずつ現れること
ヒントとして初期値がすでに何個か埋まっており、その状態から全てのマスを埋めると完成です。
イジングマシンによる解法の考え方¶
プログラムでは、数独のルールを 制約条件に変換し、イジングマシンに解かせます。
各セルに入る数字は「バイナリ変数（0/1）」として表現し、制約を満たす組み合わせを最適化によって見つけます。
制約条件の定式化¶
変数定義¶
\(y_{i,j,k} \in \{0,1\}\) を 3 次元配列で定義
\(i,j=0,...,8\) が行・列、\(k=0,...,8\) が数字（1〜9）を表す
\(y_{i,j,k}=1\) のとき、セル \((i,j)\) は \(k+1\) が入ることを意味する
これにより変数が \(9 \times 9 \times 9 = 729\) 個 作られます。
ルールからの制約¶
数独のルールは次のような one-hot 制約として定義できます：
制約 (a): 1 マスには 1 個の数字しか入らない
\[
\sum_{k=0}^8 y_{i,j,k} = 1
\]
制約 (b): 同じ行に同じ数字は 1 個
\[
\sum_{j=0}^8 y_{i,j,k} = 1
\]
制約 (c): 同じ列に同じ数字は 1 個
\[
\sum_{i=0}^8 y_{i,j,k} = 1
\]
制約 (d): 同じ \(3 \times 3\) ブロックに同じ数字は 1 個
\[
\sum_{(i,j)\in \text{block}} y_{i,j,k} = 1
\]
コード例¶
以下はサンプルコードです。
#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>
int main() {
const size_t num_size = 3;
const size_t use_num = num_size*num_size;
const size_t N = use_num;
// 初期配置を配列で表記
std::vector<std::vector<size_t>> initial =
{
{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}
};
// 変数の定義
auto y = qbpp::var("y", N, N, use_num);
// 初期配置を変数に反映するための配列
qbpp::MapList ml;
for (size_t i=0; i<N; i++){
for (size_t j=0; j<N; j++){
size_t k = initial[i][j];
if ( k != 0) {
ml.push_back({y[i][j][k-1], 1});
}
}
}
// 制約 (a): 1マスには1個の数字しか入らない
auto c_num_select = qbpp::sum(qbpp::vector_sum(y) == 1);
// 1~9が、各行と列とブロックの中に何回現れているか数える変数
auto row = qbpp::vector_sum(y, 1);
auto column = qbpp::vector_sum(y, 0);
auto block = qbpp::expr(use_num, N);
for(size_t k=0; k<use_num; k++){
for(size_t i=0; i<N; i++){
size_t i_q = i/num_size;
size_t i_r = i%num_size;
for(size_t j=0; j<N; j++){
size_t j_q = j/num_size;
size_t j_r = j%num_size;
block[k][i] += y[num_size*i_q+j_q][num_size*i_r+j_r][k];
}
}
}
// 制約 (b): 同じ行に同じ数字は 1 個
auto c_row = qbpp::sum(row==1);
// 制約 (c): 同じ列に同じ数字は 1 個
auto c_column = qbpp::sum(column==1);
// 制約 (d): 同じブロックに同じ数字は 1 個
auto c_block = qbpp::sum(block==1);
// 実際の数独の配列を作成する
auto num = qbpp::expr(use_num);
for (size_t k=0; k<use_num; k++){
num[k] = qbpp::toExpr(k+1);
}
auto x = qbpp::expr(N, N);
for(size_t i=0; i<N; i++){
for(size_t j=0; j<N; j++){
x[i][j] = qbpp::sum(num*y[i][j]);
}
}
// QUBO式の作成
auto f = c_row + c_column + c_block + c_num_select;
// 初期配列をQUBO式に反映(該当する要素を定数化)
auto g = qbpp::replace(f, ml);
// QUBO式の展開
f.simplify_as_binary();
g.simplify_as_binary();
// QUBO式を解く
auto solver = qbpp::easy_solver::EasySolver(g);
solver.time_limit(30);
solver.target_energy(0);
auto sol = solver.search();
// 初期条件を反映させた解を作成
qbpp::Sol full_sol(f);
full_sol.set(sol);
full_sol.set(ml);
// 結果の表示
const std::string hline = "+-------+-------+-------+";
for (size_t i = 0; i < N; i++) {
if (i % num_size == 0) {
std::cout << hline << std::endl;
}
for (size_t j = 0; j < N; j++) {
if (j % num_size == 0) {
std::cout << "| ";
}
std::cout << x[i][j](full_sol) << " ";
}
std::cout << "|" << std::endl;
}
std::cout << hline << std::endl;
}
import pyqbpp as qbpp
num_size = 3
use_num = num_size*num_size
N = use_num
# 初期配置を配列で表記
initial = [
[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]
]
# 変数の定義
y = qbpp.var("y", N, N, use_num)
# 初期配置を変数に反映するための配列
ml = qbpp.MapList()
for i in range(N):
for j in range(N):
k = initial[i][j]
if k != 0:
ml.add(y[i][j][k-1], 1)
# 制約 (a): 1マスには1個の数字しか入らない
c_num_select = qbpp.sum(qbpp.vector_sum(y) == 1)
# 1~9が、各行と列とブロックの中に何回現れているか数える変数
row = qbpp.vector_sum(y, 1)
column = qbpp.vector_sum(y, 0)
block = qbpp.expr(use_num, N)
for k in range(use_num):
for i in range(N):
i_q = i//num_size
i_r = i%num_size
for j in range(N):
j_q = j//num_size
j_r = j%num_size
block[k][i] += y[num_size*i_q+j_q][num_size*i_r+j_r][k]
# 制約 (b): 同じ行に同じ数字は 1 個
c_row = qbpp.sum(row==1)
# 制約 (c): 同じ列に同じ数字は 1 個
c_column = qbpp.sum(column==1)
# 制約 (d): 同じブロックに同じ数字は 1 個
c_block = qbpp.sum(block==1)
# 実際の数独の配列を作成する
num = qbpp.expr(use_num)
for k in range(use_num):
num[k] = qbpp.toExpr(k+1)
x = qbpp.expr(N, N)
for i in range(N):
for j in range(N):
x[i][j] = qbpp.sum(num*y[i][j])
# QUBO式の作成
f = c_row + c_column + c_block + c_num_select
# 初期配列をQUBO式に反映(該当する要素を定数化)
g = qbpp.replace(f, ml)
# QUBO式の展開
f.simplify_as_binary()
g.simplify_as_binary()
# QUBO式を解く
solver = qbpp.EasySolver(g)
solver.time_limit(60)
solver.target_energy(0)
sol = solver.search()
# 初期条件を反映させた解を作成
full_sol = qbpp.Sol(f).set(sol)
full_sol = full_sol.set(ml)
# 結果の表示
hline = "+-------+-------+-------+"
for i in range(N):
if i % num_size == 0:
print(hline)
for j in range(N):
if j % num_size == 0:
print("| ", end="")
print(f"{full_sol(x[i][j])}", end=" ")
print("|")
print(hline)
実行結果の例¶
+-------+-------+-------+
| 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 |
+-------+-------+-------+