# 魔方陣 :::{container} prog-cpp 3x3の魔方陣とは,1から9までの各整数をちょうど1回ずつ含む3x3の行列であり,すべての行,すべての列,および2つの対角線の和が15になるものです. 以下に例を示します: ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` ## 魔方陣を求めるための定式化 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プログラムは,これらの制約を実装し,魔方陣を求めます: ```{literalinclude} /../programFiles/cppPrograms/example/puzzle/magic-program1.cpp :language: cpp :caption: magic-program1.cpp ``` このプログラムでは,$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$ を加えて結果の方陣を出力します. このプログラムは以下の出力を生成します: ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` ## 変数の部分固定 左上のセルに値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} $$ これらの固定割り当てにより,残りのバイナリ変数の数が減少し,局所探索ベースのソルバーにとって有利になることが多いです. ## 変数の部分固定を用いた魔方陣のHi-QUBOプログラム 上記のプログラムを以下のように修正します: ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` このコードでは,`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` は完全な魔方陣を表します. このプログラムは以下の出力を生成します: ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` 左上のセルが意図通り2であることが確認できます. ## `einsum` を使った簡潔な制約構築 上のプログラムでは制約 `c2`,`c3`,`c4` を 3 重 for ループで構築していますが, これらはいずれも $3 \times 3 \times 9$ の二値配列 `x` に対するテンソル縮約 そのものなので,[`qbpp::einsum`](../advanced/Einstein-sum.md) を使えば各制約を 1 行で書けます: ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` 各 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`](../advanced/slice-connection.md) で `x` の軸 1 を 反転します.その後は同じ `"k,iik->"` パターンで反対角和が得られます. 得られる QUBO 式は for ループ版と完全に同じですが,各制約が数式定義そのまま の 1 行で書けます.サイズが大きい場合は `einsum` 内部でマルチスレッド並列に 式が構築されるため高速です. ::: :::{container} prog-python 3x3の魔方陣とは,1から9までの各整数をちょうど1回ずつ含み,すべての行,列,および2つの対角線の和が15となる3x3の行列です. 以下に例を示します： ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` ## 魔方陣を求めるための定式化 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プログラムはこれらの制約を実装し,魔方陣を求めます： ```{literalinclude} /../programFiles/pythonPrograms/example/puzzle/magic-program1.py :language: python :caption: magic-program1.py ``` このプログラムでは,$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` を見つけることでデコードされます. このプログラムの出力は以下の通りです： ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` ## 変数の部分的固定 左上のセルに値2を割り当てた解を求めたいとします. one-hotエンコーディングでは,値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} $$ これらの固定された割り当ては残りのバイナリ変数の数を減らし,局所探索ベースのソルバーにとって有益な場合が多いです. ## 変数を部分的に固定した魔方陣のPyQBPPプログラム 上記のプログラムを以下のように修正します： ```{literalinclude} /../programFiles/pythonPrograms/example/puzzle/magic-program2.py :language: python :caption: magic-program2.py ``` このコードでは,固定された割り当てを含む辞書 `ml` を作成します. 次に,元の式 `f` に対する解オブジェクト `full_sol` を作成します. `replace(f, ml)` を呼び出すと,固定された値が `f` に代入され,`ml` に含まれる変数は `g` から消えます. その結果,ソルバーが返す解 `sol` にはそれらの固定された変数が含まれません. 最後に,`set()` を使って `sol` と `ml` を `full_sol` にマージすることで,完全な割り当てを再構築します. 再構築された解 `full_sol` は完全な魔方陣を表します. このプログラムの出力は以下の通りです： ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` 意図通り,左上のセルが2であることを確認できます. ## `einsum` を使った簡潔な制約構築 上のプログラムでは制約 `c2`,`c3`,`c4` を 3 重 for ループで構築していますが, これらはいずれも $3 \times 3 \times 9$ の二値配列 `x` に対するテンソル縮約 そのものなので,[`qbpp.einsum`](../advanced/Einstein-sum.md) を使えば各制約を 1 行で書けます: ```{include} /../programFiles/markDown/example/puzzle/magic.md :start-after: :end-before: ``` 各 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` オーバーヘッドを回避でき大幅に高速になります. :::