魔方陣
3×3 の魔方陣とは、1 から 9 までの整数をそれぞれ 1 回ずつ含む 3×3 の行列であり、各行、各列、および 2 本の対角線の和が 15 となるものです。 例を以下に示します：
8 1 6
3 5 7
4 9 2
魔方陣を求めるための定式化
3×3魔方陣 $S=(s_{i,j})$ ($0\leq i,j\leq 2$) を一値符号化を用いて定式化する。 二値変数 $x_{i,j,k}$ ($0\leq i,j\leq 2, 0\leq k\leq 8$) を導入し、以下のように定義する：
\[\begin{aligned} x_{i,j,k}=1 &\Longleftrightarrow & s_{i,j}=k+1 \end{aligned}\]したがって、$X=(x_{i,j,k})$ は $3\times 3\times 9$ の二値配列となる。 以下の四つの制約を課す。
- ワンホット制約（セルごとに1つの値）： 各セル $(i,j)$ において、$x_{i,j,0}, x_{i,j,1}, \ldots,x _{i,j,8}$ のうち、ちょうど1つだけが1でなければならない：
- 各値 $k+1$ は、正確に一つのセルに現れなければならない：
-
各行および各列の和は15でなければならない： \(\begin{aligned} c_3(i): & \sum_{j=0}^2\sum_{k=0}^8 (k+1)x _{i,j,k} = 15 &(0\leq i\leq 2)\\ c_3(j): & \sum_{i=0}^2\sum_{k=0}^8 (k+1)x _{i,j,k} = 15 &(0\leq j\leq 2) \end{aligned}\)
-
対角線と対角線上の和 二つの対角線の和も15でなければならない： \(\begin{aligned} c_4: & \sum_{k=0}^8 (k+1) (x_{0,0,k}+x_{1,1,k}+x_{2,2,k}) = 15 \\ c_4: & \sum_{k=0}^8 (k+1) (x_{0,2,k}+x_{1,1,k}+x_{2,0,k}) = 15 \end{aligned}\)
すべての制約が満たされる場合、割り当て $X=(x_{i,j,k})$ は有効な3×3の魔方陣を表す。
魔方陣のためのQUBO++プログラム
以下のQUBO++プログラムはこれらの制約を実装し、魔方陣を見つける：
#include "qbpp.hpp"
#include "qbpp_easy_solver.hpp"
int main() {
auto x = qbpp::var("x", 3, 3, 9);
auto c1 = qbpp::sum(qbpp::vector_sum(x) == 1);
auto temp = qbpp::expr(9);
for (size_t i = 0; i < 3; ++i)
for (size_t j = 0; j < 3; ++j)
for (size_t k = 0; k < 9; ++k) {
temp[k] += x[i][j][k];
}
auto c2 = qbpp::sum(temp == 1);
auto row = qbpp::expr(3);
auto column = qbpp::expr(3);
for (size_t i = 0; i < 3; ++i)
for (size_t j = 0; j < 3; ++j)
for (size_t k = 0; k < 9; ++k) {
row[i] += (k + 1) * x[i][j][k];
column[j] += (k + 1) * x[i][j][k];
}
auto c3 = qbpp::sum(row == 15) + qbpp::sum(column == 15);
auto diag = qbpp::Expr(0);
for (size_t k = 0; k < 9; ++k)
diag += (k + 1) * (x[0][0][k] + x[1][1][k] + x[2][2][k]);
auto anti_diag = qbpp::Expr(0);
for (size_t k = 0; k < 9; ++k)
anti_diag += (k + 1) * (x[0][2][k] + x[1][1][k] + x[2][0][k]);
auto c4 = (diag == 15) + (anti_diag == 15);
auto f = c1 + c2 + c3 + c4;
f.simplify_as_binary();
auto solver = qbpp::easy_solver::EasySolver(f);
solver.target_energy(0);
auto sol = solver.search();
auto result = qbpp::onehot_to_int(sol(x));
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
std::cout << result[i][j] + 1 << " ";
}
std::cout << std::endl;
}
}
このプログラムでは、$3\times 3\times9$ の二値変数配列を定義する xを定義する。
次に、4つの制約式 c1, c2, c3,および c4、そしてこれらを結合して fとします。
式 f は、すべての制約が満たされたときに最小エネルギー 0 を達成します。
Easy Solverオブジェクトsolverを作成し f を作成し、目標エネルギーを0に設定します。これにより、実行可能な（最適）解が見つかり次第、探索は終了します。
返された解は solに格納する。
最後に、one-hot表現を整数に変換するため qbpp::onehot_to_int()を使用し、整数の $3\times 3$ 配列を返します。
各要素に $1$ を加算して結果の正方形を出力します。
このプログラムは以下の出力を生成する:
8 1 6
3 5 7
4 9 2
変数の部分的な固定
左上セルに値2を割り当てる解を求める場合を考えます。 ワンホット符号化では値2は$k=1$に対応するため、
\[\begin{aligned} x_{0,0,k} &=1 & {\rm if\,\,} k=1\\ x_{0,0,k} &=0 & {\rm if\,\,} k\neq 1 \end{aligned}\]さらに、制約 $c_2$ は各数 $k+1$ が正確に1回だけ出現することを強制するため、これを固定すると 他のセルが値2を取ることは即座に不可能となる。 したがって、以下も固定できる：
\[\begin{aligned} x_{i,j,1} &=0 & {\rm if\,\,} (i,j)\neq (0,0)\\ \end{aligned}\]これらの固定値により残りの二進変数の数が減少し、局所探索ベースのソルバーにとって有益な場合が多い。
変数を部分的に固定した魔方陣のQUBO++プログラム
上記プログラムを以下のように修正する：
qbpp::MapList ml;
for (size_t k = 0; k < 9; ++k) {
if (k == 1)
ml.push_back({x[0][0][k], 1});
else
ml.push_back({x[0][0][k], 0});
}
for (size_t i = 0; i < 3; ++i)
for (size_t j = 0; j < 3; ++j) {
if (!(i == 0 && j == 0)) {
ml.push_back({x[i][j][1], 0});
}
}
qbpp::Sol full_sol(f);
f.replace(ml);
f.simplify_as_binary();
auto solver = qbpp::easy_solver::EasySolver(f);
solver.target_energy(0);
auto sol = solver.search();
full_sol.set(sol);
full_sol.set(ml);
auto result = qbpp::onehot_to_int(full_sol(x));
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
std::cout << result[i][j] + 1 << " ";
}
std::cout << std::endl;
}
このコードでは、 qbpp::MapList オブジェクト ml を作成し、 push_back()を使用します。
次に full_sol元の式に対する解オブジェクト fを作成します。 f.replace(ml) は固定値を f に直接置換されるため、 ml から消えます。結果として、 f。
結果として、ソルバーが返す解 sol が返す解にはそれらの固定変数は含まれない。
最後に、 sol と ml を結合して full_sol 経由で set()によってとをに統合することで完全な割り当てを再構築する。
再構築された解 full_sol は完全な魔方陣を表します。
このプログラムは以下の出力を生成する:
2 7 6
9 5 1
4 3 8
左上のセルが意図通り2であることを確認できます。