魔方陣

3x3の魔方陣とは,1から9までの各整数をちょうど1回ずつ含む3x3の行列であり,すべての行,すべての列,および2つの対角線の和が15になるものです. 以下に例を示します:

8 1 6
3 5 7
4 9 2

魔方陣を求めるための定式化

3x3の魔方陣 \(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\) のバイナリ配列です. 以下の4つの制約を課します.

1. ワンホット制約（各セルに1つの値）: 各セル \((i,j)\) について,\(x_{i,j,0}, x_{i,j,1}, \ldots,x _{i,j,8}\) のうちちょうど1つが1でなければなりません:

\[ \begin{aligned} c_1(i,j): & \sum_{k=0}^8 x _{i,j,k}=1 & (0\leq i,j\leq 2) \end{aligned} \]

2. 各値 \(k+1\) はちょうど1つのセルに現れなければなりません:

\[ \begin{aligned} c_2(k): & \sum_{i=0}^2\sum_{j=0}^2x _{i,j,k}=1 & (0\leq k\leq 8) \end{aligned} \]

3. 各行および各列の和は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} \)$

4. 対角線と反対角線の和 2つの対角線の和も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})\) は有効な3x3の魔方陣を表します.

魔方陣のためのHi-QUBOプログラム

以下のHi-QUBOプログラムは,これらの制約を実装し,魔方陣を求めます:

magic-program1.cpp

#include <qbpp/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::toExpr(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::toExpr(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::EasySolver(f);
  auto sol = solver.search({{"target_energy", 0}});
  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を達成します.

f に対するEasy Solverオブジェクト solver を作成し,目標エネルギーを0に設定することで,実行可能（最適）解が見つかり次第,探索が終了します. 返された解は sol に格納されます. 最後に,qbpp::onehot_to_int() を使用してワンホット表現を整数に変換します.この関数は \(\{0,1, \ldots, 8\}\) の整数からなる \(3\times 3\) の配列を返します.各要素に \(1\) を加えて結果の方陣を出力します.

このプログラムは以下の出力を生成します:

8 1 6
3 5 7
4 9 2

変数の部分固定

左上のセルに値2が割り当てられた解を求めたいとします. ワンホット符号化では,値2は \(k=1\) に対応するため,次のように固定します:

\[\begin{split} \begin{aligned} x_{0,0,k} &=1 & {\rm if\,\,} k=1\\ x_{0,0,k} &=0 & {\rm if\,\,} k\neq 1 \end{aligned} \end{split}\]

さらに,制約 \(c_2\) により各数 \(k+1\) はちょうど1回だけ出現するため,固定することで他のセルが値2を取れなくなります. したがって,次のようにも固定できます:

\[\begin{split} \begin{aligned} x_{i,j,1} &=0 & {\rm if\,\,} (i,j)\neq (0,0)\\ \end{aligned} \end{split}\]

これらの固定割り当てにより,残りのバイナリ変数の数が減少し,局所探索ベースのソルバーにとって有利になることが多いです.

変数の部分固定を用いた魔方陣のHi-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});
      }
    }

  auto g = qbpp::replace(f, ml);
  g.simplify_as_binary();

  auto solver = qbpp::EasySolver(g);
  auto sol = solver.search({{"target_energy", 0}});

  auto full_sol = qbpp::Sol(f).set(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() を使用して固定割り当てを追加しています. 次に,qbpp::replace(f, ml) を呼び出して固定値を代入し,元の f を変更せずに新しい式 g を生成します. ml に含まれる変数は g から消えます. Easy Solverを g に適用し,解 sol にはそれらの固定変数は含まれません. 最後に,ゼロ初期化された qbpp::Sol(f) に set(sol) と set(ml) をチェーンして完全な解を構築します. 得られた full_sol は完全な魔方陣を表します.

このプログラムは以下の出力を生成します:

2 7 6
9 5 1
4 3 8

左上のセルが意図通り2であることが確認できます.

einsum を使った簡潔な制約構築

上のプログラムでは制約 c2,c3,c4 を 3 重 for ループで構築していますが, これらはいずれも \(3 \times 3 \times 9\) の二値配列 x に対するテンソル縮約 そのものなので,qbpp::einsum を使えば各制約を 1 行で書けます:

  auto vals = qbpp::array({1, 2, 3, 4, 5, 6, 7, 8, 9});

  // c2: 各値 k がちょうど 1 回現れる. temp[k] = Σ_{i,j} x[i,j,k]
  auto c2 = qbpp::sum(qbpp::einsum<1>("ijk->k", x) == 1);

  // c3: row[i] = Σ_{j,k} (k+1) x[i,j,k]、column[j] = Σ_{i,k} (k+1) x[i,j,k]
  auto row    = qbpp::einsum<1>("k,ijk->i", vals, x);
  auto column = qbpp::einsum<1>("k,ijk->j", vals, x);
  auto c3     = qbpp::sum(row == 15) + qbpp::sum(column == 15);

  // c4: 対角線 Σ_k (k+1) Σ_i x[i,i,k]   — "ii" で軸 0 と軸 1 を結合
  auto diag      = qbpp::einsum<0>("k,iik->", vals, x);

  // 反対角線: 同じパターンだが軸 1 を反転する必要がある.
  // slice と concat で反転 — 単一要素スライスは軸を size 1 で残すので、
  // 軸 1 で順番を入れ替えて concat すれば 3×3×9 の反転配列ができる.
  auto x_flip = qbpp::concat(
      qbpp::concat(x(qbpp::all, qbpp::slice(2), qbpp::all),
                   x(qbpp::all, qbpp::slice(1), qbpp::all), 1),
      x(qbpp::all, qbpp::slice(0), qbpp::all), 1);
  auto anti_diag = qbpp::einsum<0>("k,iik->", vals, x_flip);
  auto c4 = (diag == 15) + (anti_diag == 15);

各 subscript の読み方:

• "ijk->k"（c2）— i と j を縮約,k を残す.

• "k,ijk->i"（row）— vals（軸 k）と x（軸 i, j, k）の間で j, k を縮約,i を残す.

• "k,ijk->j"（column）— row と同じだが j を残す.

• "k,iik->"（対角線）— x 内で ii のラベル繰り返しが軸 0 と軸 1 を結合し（x[i,i,k]）,結果はスカラー（k と対角の i 両方を縮約）.

反対角は x[i, n-1-i, k] が必要で,これは einsum の subscript で直接表現 できないため,まず slice と concat で x の軸 1 を 反転します.その後は同じ "k,iik->" パターンで反対角和が得られます.

得られる QUBO 式は for ループ版と完全に同じですが,各制約が数式定義そのまま の 1 行で書けます.サイズが大きい場合は einsum 内部でマルチスレッド並列に 式が構築されるため高速です.

3x3の魔方陣とは,1から9までの各整数をちょうど1回ずつ含み,すべての行,列,および2つの対角線の和が15となる3x3の行列です. 以下に例を示します：

8 1 6
3 5 7
4 9 2

魔方陣を求めるための定式化

one-hotエンコーディングを用いて,3x3の魔方陣 \(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\) のバイナリ配列です. 以下の4つの制約を課します.

1. one-hot制約（各セルに1つの値）： 各セル \((i,j)\) について,\(x_{i,j,0}, x_{i,j,1}, \ldots,x _{i,j,8}\) のうちちょうど1つが1でなければなりません：

\[ \begin{aligned} c_1(i,j): & \sum_{k=0}^8 x _{i,j,k}=1 & (0\leq i,j\leq 2) \end{aligned} \]

2. 各値 \(k+1\) はちょうど1つのセルに現れなければなりません：

\[ \begin{aligned} c_2(k): & \sum_{i=0}^2\sum_{j=0}^2x _{i,j,k}=1 & (0\leq k\leq 8) \end{aligned} \]

3. 各行と各列の和は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} \)$

4. 対角線と反対角線の和 2つの対角線の和も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})\) は有効な3x3の魔方陣を表します.

魔方陣のPyQBPPプログラム

以下のPyQBPPプログラムはこれらの制約を実装し,魔方陣を求めます：

magic-program1.py

import pyqbpp as qbpp

x = qbpp.var("x", shape=(3, 3, 9))

c1 = qbpp.sum(qbpp.vector_sum(x) == 1)

temp = qbpp.expr(shape=9)
for i in range(3):
    for j in range(3):
        for k in range(9):
            temp[k] += x[i][j][k]
c2 = qbpp.sum(temp == 1)

row = qbpp.expr(shape=3)
column = qbpp.expr(shape=3)
for i in range(3):
    for j in range(3):
        for k in range(9):
            row[i] += (k + 1) * x[i][j][k]
            column[j] += (k + 1) * x[i][j][k]
c3 = qbpp.sum(row == 15) + qbpp.sum(column == 15)

diag = 0
for k in range(9):
    diag += (k + 1) * (x[0][0][k] + x[1][1][k] + x[2][2][k])
anti_diag = 0
for k in range(9):
    anti_diag += (k + 1) * (x[0][2][k] + x[1][1][k] + x[2][0][k])
c4 = (diag == 15) + (anti_diag == 15)

f = c1 + c2 + c3 + c4
f.simplify_as_binary()

solver = qbpp.EasySolver(f)
sol = solver.search(target_energy=0)
for i in range(3):
    for j in range(3):
        val = next(k for k in range(9) if sol(x[i][j][k]) == 1)
        print(val + 1, end=" ")
    print()

このプログラムでは,\(3\times 3\times9\) のバイナリ変数配列 x を定義します. 次に,4つの制約式 c1,c2,c3,c4 を構築し,それらを f にまとめます. 式 f はすべての制約が満たされたとき最小エネルギー0を達成します.

f に対するEasy Solverオブジェクト solver を作成し,search() に target_energy=0 を渡します.これにより,実行可能（最適）解が見つかり次第,探索が終了します. 得られたone-hotエンコーディングは,sol(x[i][j][k]) == 1 となるインデックス k を見つけることでデコードされます.

このプログラムの出力は以下の通りです：

8 1 6
3 5 7
4 9 2

変数の部分的固定

左上のセルに値2を割り当てた解を求めたいとします. one-hotエンコーディングでは,値2は \(k=1\) に対応するため,以下を固定します：

\[\begin{split} \begin{aligned} x_{0,0,k} &=1 & {\rm if\,\,} k=1\\ x_{0,0,k} &=0 & {\rm if\,\,} k\neq 1 \end{aligned} \end{split}\]

さらに,制約 \(c_2\) が各数 \(k+1\) がちょうど1回現れることを強制するため,この固定は他のどのセルも値2を取れないことを直ちに意味します. したがって,以下も固定できます：

\[\begin{split} \begin{aligned} x_{i,j,1} &=0 & {\rm if\,\,} (i,j)\neq (0,0)\\ \end{aligned} \end{split}\]

これらの固定された割り当ては残りのバイナリ変数の数を減らし,局所探索ベースのソルバーにとって有益な場合が多いです.

変数を部分的に固定した魔方陣のPyQBPPプログラム

上記のプログラムを以下のように修正します：

magic-program2.py

import pyqbpp as qbpp

x = qbpp.var("x", shape=(3, 3, 9))

c1 = qbpp.sum(qbpp.vector_sum(x) == 1)

temp = qbpp.expr(shape=9)
for i in range(3):
    for j in range(3):
        for k in range(9):
            temp[k] += x[i][j][k]
c2 = qbpp.sum(temp == 1)

row = qbpp.expr(shape=3)
column = qbpp.expr(shape=3)
for i in range(3):
    for j in range(3):
        for k in range(9):
            row[i] += (k + 1) * x[i][j][k]
            column[j] += (k + 1) * x[i][j][k]
c3 = qbpp.sum(row == 15) + qbpp.sum(column == 15)

diag = 0
for k in range(9):
    diag += (k + 1) * (x[0][0][k] + x[1][1][k] + x[2][2][k])
anti_diag = 0
for k in range(9):
    anti_diag += (k + 1) * (x[0][2][k] + x[1][1][k] + x[2][0][k])
c4 = (diag == 15) + (anti_diag == 15)

f = c1 + c2 + c3 + c4
f.simplify_as_binary()

ml = {x[0][0][k]: 1 if k == 1 else 0 for k in range(9)}
ml.update({x[i][j][1]: 0 for i in range(3) for j in range(3) if not (i == 0 and j == 0)})

g = qbpp.replace(f, ml)
g.simplify_as_binary()

solver = qbpp.EasySolver(g)
sol = solver.search(target_energy=0)
full_sol = qbpp.Sol(f).set(sol, ml)

for i in range(3):
    for j in range(3):
        val = next(k for k in range(9) if full_sol(x[i][j][k]) == 1)
        print(val + 1, end=" ")
    print()

このコードでは,固定された割り当てを含む辞書 ml を作成します. 次に,元の式 f に対する解オブジェクト full_sol を作成します. replace(f, ml) を呼び出すと,固定された値が f に代入され,ml に含まれる変数は g から消えます. その結果,ソルバーが返す解 sol にはそれらの固定された変数が含まれません. 最後に,set() を使って sol と ml を full_sol にマージすることで,完全な割り当てを再構築します. 再構築された解 full_sol は完全な魔方陣を表します.

このプログラムの出力は以下の通りです：

2 7 6
9 5 1
4 3 8

意図通り,左上のセルが2であることを確認できます.

einsum を使った簡潔な制約構築

上のプログラムでは制約 c2,c3,c4 を 3 重 for ループで構築していますが, これらはいずれも \(3 \times 3 \times 9\) の二値配列 x に対するテンソル縮約 そのものなので,qbpp.einsum を使えば各制約を 1 行で書けます:

vals = qbpp.array([1, 2, 3, 4, 5, 6, 7, 8, 9])

# c2: 各値 k がちょうど 1 回現れる. temp[k] = Σ_{i,j} x[i,j,k]
c2 = qbpp.sum(qbpp.einsum("ijk->k", x) == 1)

# c3: row[i] = Σ_{j,k} (k+1) x[i,j,k]、column[j] = Σ_{i,k} (k+1) x[i,j,k]
row    = qbpp.einsum("k,ijk->i", vals, x)
column = qbpp.einsum("k,ijk->j", vals, x)
c3 = qbpp.sum(row    == 15) + \
     qbpp.sum(column == 15)

# c4: 対角線 Σ_k (k+1) Σ_i x[i,i,k]  — "ii" で軸 0 と軸 1 を結合
diag = qbpp.einsum("k,iik->", vals, x)

# 反対角線: 同じパターンだが軸 1 を反転する必要がある。
# スライスと concat で x を軸 1 に沿って反転する。
x_flip = qbpp.concat([x[:, 2:3, :], x[:, 1:2, :], x[:, 0:1, :]], axis=1)
anti_diag = qbpp.einsum("k,iik->", vals, x_flip)
c4 = (diag == 15) + (anti_diag == 15)

各 subscript の読み方:

• "ijk->k"（c2）— i と j を縮約,k を残す.

• "k,ijk->i"（row）— vals（軸 k）と x（軸 i, j, k）の間で j, k を縮約,i を残す.

• "k,ijk->j"（column）— row と同じだが j を残す.

• "k,iik->"（対角線）— x 内で ii のラベル繰り返しが軸 0 と軸 1 を結合し（x[i,i,k]）,結果はスカラー（k と対角の i 両方を縮約）.

反対角は x[i, n-1-i, k] が必要で,これは einsum の subscript で直接表現 できないため,Python のスライス構文（x[:, 2:3, :],x[:, 1:2, :], x[:, 0:1, :] — 単一要素スライスは軸を残します）と qbpp.concat(..., axis=1) で x の軸 1 を反転します.その後は同じ "k,iik->" パターンで反対角和が得られます.

得られる QUBO 式は for ループ版と完全に同じですが,各制約が数式定義そのまま の 1 行で書けます.サイズが大きい場合は,einsum の処理が C++ バックエンド内でマルチスレッド実行されるため,for ループ版で 1 反復ごとに 発生する Python の ctypes オーバーヘッドを回避でき大幅に高速になります.