魔方陣¶

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

1 6
5 7
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つの値）: 各セル $(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} \]

各値 $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} \]

各行および各列の和は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} $$
対角線と反対角線の和 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$ を加えて結果の方陣を出力します.

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

変数の部分固定

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

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

左上のセルが意図通り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の行列です. 以下に例を示します：

1 6
5 7
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つの制約を課します.

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} \]

各値 $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} \]

各行と各列の和は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} $$
対角線と反対角線の和 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 を見つけることでデコードされます.

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

1 6
5 7
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 は完全な魔方陣を表します.

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

7 6
5 1
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 オーバーヘッドを回避でき大幅に高速になります.