# 置換行列の生成
:::{container} prog-cpp
多くの組合せ最適化問題は,最適な置換を見つけることが目的であるという意味で置換ベースです.
このような最適化問題を定式化する基本的な手法として,QUBO定式化ではバイナリ変数の行列が使用されます.
## 置換行列
$X=(x_{i,j})$ ($0\leq i,j\leq n-1$) を $n\times n$ のバイナリ値の行列とします.
行列 $X$ は,以下に示すように,すべての行とすべての列にちょうど1つの1のエントリがある場合にのみ**置換行列**と呼ばれます.
**置換行列**は $n$ 個の数 $(0,1,\ldots,n-1)$ の置換を表し,$x_{i,j} = 1$ であるのは $i$ 番目の要素が $j$ である場合に限ります.
例えば,上記の置換行列は置換 $(1,3,0,2)$ を表しています.
## 置換行列のQUBO定式化
バイナリ変数行列 $X=(x_{i,j})$ ($0\leq i,j\leq n-1$) が置換行列を格納するのは,各行と各列の和が1である場合に限ります.
したがって,以下のQUBO関数は $X$ が置換行列を格納する場合にのみ最小値0を取ります：
$$
\begin{aligned}
f(X) &= \sum_{i=0}^{n-1}\left(1-\sum_{j=0}^{n-1}x_{i,j}\right)^2+\sum_{j=0}^{n-1}\left(1-\sum_{i=0}^{n-1}x_{i,j}\right)^2
\end{aligned}
$$
## 置換行列を生成するHi-QUBOプログラム
上記の式 $f(X)$ に基づいて,以下のようにHi-QUBOプログラムを設計できます：
```{literalinclude} /../programFiles/cppPrograms/example/math/permutation-program1.cpp
:language: cpp
:caption: permutation-program1.cpp
```
このプログラムでは,**`qbpp::var("x",4,4)`**は**`x`**という名前の形状 $4\times 4$ の `Var` の配列を返します.
`x` の厳密な型は **`qbpp::Array<2, qbpp::Var>`** で,`2` は次元数,`qbpp::Var` は要素の型を表します.
`qbpp::Expr` オブジェクト**`f`**に対して,2つの二重forループが $f(X)$ の式を追加します.
Exhaustive Solverを使用して,すべての最適解が計算され**`sols`**に格納されます.
`sols` 内のすべての解が1つずつ表示されます.
ここで,`sol(x)` は `sol` における `x` の値の行列を返します.
このプログラムは以下のように24個すべての置換を出力します：
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
> **注意**
> `qbpp::var("x",4,4)` は $4\times 4$ のバイナリ変数配列を生成します.
> 各 `qbpp::Var` オブジェクトは `x[i][j]` として表され,`sol` における `x[i][j]` の値は `sol(x[i][j])` または `x[i][j](sol)` のいずれかで取得できます.
## 配列関数と演算を使った置換行列のQUBO定式化
`qbpp::vector_sum()` を使用して,バイナリ変数の行列 `x` の行方向と列方向の和を計算できます：
- **`qbpp::vector_sum(x, 1)`**: `x` の各行の和を計算し,それらの和を含むサイズ `n` の配列を返します.
- **`qbpp::vector_sum(x, 0)`**: `x` の各列の和を計算し,それらの和を含むサイズ `n` の配列を返します.
> **注意**:
> 多次元配列 `x` と軸 `k` に対して,`qbpp::vector_sum(x, k)` は軸 `k` に沿った和を計算し,次元が1つ減った多次元配列を返します.
> 2次元配列（行列）`x` の場合,軸 `1` は行方向に,軸 `0` は列方向に対応します.
スカラー-配列演算を使用して,各要素から1を引くことができます：
- **`qbpp::vector_sum(x, 1) - 1`**: 行方向の各和から1を引きます.
- **`qbpp::vector_sum(x, 0) - 1`**: 列方向の各和から1を引きます.
これら2つのサイズ `n` の配列に対して,`qbpp::sqr()` は各要素を2乗し,`qbpp::sum()` はすべての要素の和を計算します.
以下のHi-QUBOプログラムは,これらの配列関数と演算を使用してQUBO定式化を実装しています：
```{literalinclude} /../programFiles/cppPrograms/example/math/permutation-program2.cpp
:language: cpp
:caption: permutation-program2.cpp
```
このプログラムでは,`x(sol)` は `sol` における `x` に割り当てられた値の行列を返します.これは `x` と同じ形状の整数値の行列です.
`qbpp::onehot_to_int()` は軸に沿ったone-hot配列を対応する整数に変換します.
- **`qbpp::onehot_to_int(x(sol), 1)`**: 各行に対応する整数を計算し,4つの整数の配列として返します.これが置換を表します.
- **`qbpp::onehot_to_int(x(sol), 0)`**: 各列に対応する整数を返し,4つの整数の配列として返します.これが置換の逆を表します.
このプログラムはすべての置換とその逆を整数ベクトルとして以下のように出力します：
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
## 割当問題とそのQUBO定式化
$C = (c_{i,j})$ をサイズ $n \times n$ のコスト行列とします.
$C$ に対する**割当問題**は,総コストを最小化する置換
$p:\lbrace 0,1,\ldots, n-1\rbrace \rightarrow \lbrace 0,1,\ldots, n-1\rbrace$
を見つけることです：
$$
\begin{aligned}
g(p) &= \sum_{i=0}^{n-1}c_{i,p(i)}
\end{aligned}
$$
この問題のQUBO定式化には,サイズ $n \times n$ の置換行列 $X = (x_{i,j})$ を使用して以下のように定義できます：
$$
\begin{aligned}
g(X) &= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j}
\end{aligned}
$$
明らかに,$X$ が置換 $p$ を表す場合にのみ $g(p) = g(X)$ が成り立ちます.
置換行列のQUBO定式化 $f(X)$ と総コスト $g(X)$ を組み合わせて,割当問題のQUBO定式化を得ます：
$$
\begin{aligned}
h(X) &= P\cdot f(X)+g(X) \\
&=P\left(\sum_{i=0}^{n-1}\left(1-\sum_{j=0}^{n-1}x_{i,j}\right)^2+\sum_{j=0}^{n-1}\left(1-\sum_{i=0}^{n-1}x_{i,j}\right)^2\right)+\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j}
\end{aligned}
$$
ここで,$P$ は $f(X)$ にエンコードされた置換制約を優先するための十分に大きな正の定数です.
## 割当問題のHi-QUBOプログラム
これで割当問題のHi-QUBOプログラムを設計する準備が整いました.
このプログラムでは,サイズ $4\times4$ の固定行列 $C$ が整数定数の配列として与えられます.
`qbpp::array({...})` は入れ子の初期化子リストから整数定数配列を構築します.形状はリストの入れ子の深さから自動判定され（リストのリストなら 2 次元配列になる）,ここでは $4\times4$ の `Array<2, coeff_t>` が生成されます.
$f(X)$ と $g(X)$ の式は配列関数と演算を使用して定義されます.
ここで,`qbpp::vector_sum(x, 1) == 1` は等式が満たされた場合に最小値0を取るQUBO式の配列を返します.
実際には,`qbpp::sqr(qbpp::vector_sum(x, 1) - 1)` と同じQUBO式を返します.
また,`c * x` は `c` と `x` の要素ごとの積を計算して得られる行列を返すため,`qbpp::sum(c * x)` は `g(X)` を返します.
```{literalinclude} /../programFiles/cppPrograms/example/math/permutation-program3.cpp
:language: cpp
:caption: permutation-program3.cpp
```
Easy Solverを使用して `h` の解を求めます.
`h` に対するEasy Solverオブジェクト `solver` について,`search()` に {% raw %}`{{"time_limit", 1.0}}`{% endraw %} を渡して解の探索の制限時間を1.0秒に設定します.
得られた置換は `result` に格納され,選択された `c[i][j]` の値が順に出力されます.
このプログラムの出力は以下のとおりです：
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
> **注意**
> 式 `f` と整数 `m` に対して,`f == m` は式 `qbpp::sqr(f - m)` を返します.
> これは等式 `f == m` が満たされる場合にのみ最小値0を取ります.
## ネイティブ制約 `cons()` による割当問題の記述
上のプログラムでは,置換制約 `f` を重み1000のペナルティ式として目的関数に加えました.
Hi-QUBO では,行と列の one-hot 制約を `qbpp::cons()` で囲んで**制約であることを明示**できます:
```{literalinclude} /../programFiles/cppPrograms/example/math/permutation-program4.cpp
:language: cpp
:caption: permutation-program4.cpp
```
配列の比較 `qbpp::vector_sum(x, 1) == 1` を `qbpp::cons()` で囲むと,
**要素ごと**（ここでは行ごと・列ごと）**に1本の制約**として宣言されます.
宣言された制約は**制約として特別に処理**され,バンドルされたソルバーは
制約を満たす置換行列を効率よく探索します.
`h.cons(sol)` は解 `sol` で違反している制約の本数を返します（0 なら `x` は置換行列）.
このプログラムの出力は以下のとおりです:
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
`cons()` の詳しい使い方は[ネイティブ制約](../basic/constraints.md)をご覧ください.
:::
:::{container} prog-python
多くの組合せ最適化問題は,最適な置換を見つけることが目的であるという意味で,置換に基づいています.
このような最適化問題を定式化するための基本的な手法として,QUBO定式化ではバイナリ変数の行列が使用されます.
## 置換行列
$X=(x_{i,j})$ ($0\leq i,j\leq n-1$) を $n\times n$ のバイナリ値の行列とします.
行列 $X$ が**置換行列**であるとは,以下に示すように,すべての行とすべての列にちょうど1つの1のエントリがある場合に限ります.
**置換行列**は $n$ 個の数 $(0,1,\ldots,n-1)$ の置換を表し,$x_{i,j} = 1$ であることと $i$ 番目の要素が $j$ であることが同値です.
例えば,上記の置換行列は置換 $(1,3,0,2)$ を表しています.
## 置換行列のQUBO定式化
バイナリ変数行列 $X=(x_{i,j})$ ($0\leq i,j\leq n-1$) が置換行列を格納しているのは,各行と各列の合計が1である場合に限ります.
したがって,以下のQUBO関数は $X$ が置換行列を格納している場合に限り最小値0をとります:
$$
\begin{aligned}
f(X) &= \sum_{i=0}^{n-1}\left(1-\sum_{j=0}^{n-1}x_{i,j}\right)^2+\sum_{j=0}^{n-1}\left(1-\sum_{i=0}^{n-1}x_{i,j}\right)^2
\end{aligned}
$$
## 置換行列を生成するPyQBPPプログラム
バイナリ変数の2次元配列は,`shape=` キーワード引数にタプルを渡して作成します.
例えば,**`qbpp.var("x", shape=(4, 4))`** は $4\times 4$ のバイナリ変数 array
を返し,各要素には `x[i][j]` でアクセスできます.
`shape=(2, 3, 4)` のようなタプルで高次元配列も同様に作成できます.
詳細は [多次元変数](../advanced/multi-dimensional-variables-and-expressions.md) を参照してください.
上記の式 $f(X)$ に基づいて,以下のようにPyQBPPプログラムを設計できます:
```{literalinclude} /../programFiles/pythonPrograms/example/math/permutation-program1.py
:language: python
:caption: permutation-program1.py
```
このプログラムでは,**`qbpp.var("x", shape=(4, 4))`** は **`x`** という名前の形状 $\{4, 4\}$ の array オブジェクトを返します.
`Expr` オブジェクト **`f`** に対して,2つの二重forループが $f(X)$ の式を追加します.
Exhaustive Solverを使用して,すべての最適解が計算され **`result.sols`** に格納されます.
`result.sols` 内のすべての解が1つずつ表示されます.
ここで,`sol(x)` は `sol` における `x` の値の行列 (array of int) を返します.
このプログラムは以下のように24個すべての置換を出力します:
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
> **注釈**
> バイナリ変数の行列は array クラスを使用した多次元配列として実装されています.
> 例えば,`qbpp.var("x", shape=(4, 4))` は形状 `(4, 4)` の array オブジェクトを返します.
> 各 `Var` オブジェクトは `x[i][j]` として表され,`sol` における `x[i][j]` の値は `sol(x[i][j])` または `x[i][j](sol)` のいずれかで取得できます.
## 配列関数と演算を使った置換行列のQUBO定式化
**`qbpp.vector_sum()`** を使用して,バイナリ変数の行列 `x` の行方向と列方向の和を計算できます:
- **`qbpp.vector_sum(x, axis=1)`**: `x` の各行の和を計算し,それらの和を含むサイズ `n` の配列を返します.
- **`qbpp.vector_sum(x, axis=0)`**: `x` の各列の和を計算し,それらの和を含むサイズ `n` の配列を返します.
> **注釈**:
> 多次元配列 `x` と軸 `k` に対して,`qbpp.vector_sum(x, axis=k)` は軸 `k` に沿った和を計算し,次元が1つ減った多次元配列を返します.
> 2次元配列（行列）`x` の場合,軸 `1` は行方向に,軸 `0` は列方向に対応します.
スカラー-配列演算を使用して,各要素から1を引くことができます:
- **`qbpp.vector_sum(x, axis=1) - 1`**: 行方向の各和から1を引きます.
- **`qbpp.vector_sum(x, axis=0) - 1`**: 列方向の各和から1を引きます.
これら2つのサイズ `n` の配列に対して,`qbpp.sqr()` は各要素を2乗し,`qbpp.sum()` はすべての要素の和を計算します.
以下のPyQBPPプログラムは,これらの配列関数と演算を使用してQUBO定式化を実装しています:
```{literalinclude} /../programFiles/pythonPrograms/example/math/permutation-program2.py
:language: python
:caption: permutation-program2.py
```
このプログラムでは,`sol(x)` は `sol` における `x` に割り当てられた値の行列を返します.これは `x` と同じ形状の整数値の行列です.
`qbpp.onehot_to_int()` は軸に沿ったone-hot配列を対応する整数に変換します.
- **`qbpp.onehot_to_int(sol(x), axis=1)`**: 各行に対応する整数を計算し,4つの整数の配列として返します.これが置換を表します.
- **`qbpp.onehot_to_int(sol(x), axis=0)`**: 各列に対応する整数を返し,4つの整数の配列として返します.これが置換の逆を表します.
このプログラムはすべての置換とその逆を整数ベクトルとして以下のように出力します:
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
## 割当問題とそのQUBO定式化
$C = (c_{i,j})$ をサイズ $n \times n$ のコスト行列とします.
$C$ に対する**割当問題**は,総コストを最小化する置換
$p:\lbrace 0,1,\ldots, n-1\rbrace \rightarrow \lbrace 0,1,\ldots, n-1\rbrace$
を見つける問題です:
$$
\begin{aligned}
g(p) &= \sum_{i=0}^{n-1}c_{i,p(i)}
\end{aligned}
$$
この問題のQUBO定式化には,サイズ $n \times n$ の置換行列 $X = (x_{i,j})$ を使用し,以下のように定義します:
$$
\begin{aligned}
g(X) &= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j}
\end{aligned}
$$
明らかに,$X$ が置換 $p$ を表す場合にのみ $g(p) = g(X)$ が成り立ちます.
置換行列のQUBO定式化 $f(X)$ と総コスト $g(X)$ を組み合わせて,割当問題のQUBO定式化を得ます:
$$
\begin{aligned}
h(X) &= P\cdot f(X)+g(X) \\
&=P\left(\sum_{i=0}^{n-1}\left(1-\sum_{j=0}^{n-1}x_{i,j}\right)^2+\sum_{j=0}^{n-1}\left(1-\sum_{i=0}^{n-1}x_{i,j}\right)^2\right)+\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j}
\end{aligned}
$$
ここで,$P$ は $f(X)$ にエンコードされた置換制約を優先するための十分に大きな正の定数です.
## 割当問題のPyQBPPプログラム
これで割当問題のPyQBPPプログラムを設計する準備が整いました.
このプログラムでは,サイズ $4\times4$ の固定行列 $C$ が配列として与えられます.
`qbpp.array()` はネストされたPythonリストを自動的にネストされた配列オブジェクトに変換するため,多次元配列を簡潔に作成できます.
$f(X)$ と $g(X)$ の式は配列関数と演算を使用して定義されます.
ここで,`qbpp.constrain(qbpp.vector_sum(x, axis=1), equal=1)` は等式 `vector_sum(x, axis=1) == 1` が満たされた場合に最小値0を取るQUBO式の配列を返します.
実際には,`qbpp.sqr(qbpp.vector_sum(x, axis=1) - 1)` と同じQUBO式を返します.
また,`c * x` は `c` と `x` の要素ごとの積を計算して得られる行列を返すため,`qbpp.sum(c * x)` は `g(X)` を返します.
```{literalinclude} /../programFiles/pythonPrograms/example/math/permutation-program3.py
:language: python
:caption: permutation-program3.py
```
Easy Solverを使用して `h` の解を求めます.
`h` に対するEasy Solverオブジェクト `solver` について,`search()` に `time_limit=1.0` を渡すことで解の探索の制限時間を1.0秒に設定しています.
得られた置換は `result` に格納され,選択された `c[i][j]` の値が順に出力されます.
このプログラムの出力は以下のとおりです:
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
> **注釈**
> 式 `f` と整数 `m` に対して,`qbpp.constrain(f, equal=m)` は式 `sqr(f - m)` を返します.
> これは等式 `f == m` が満たされる場合に限り最小値0をとります.
## ネイティブ制約 `cons()` による割当問題の記述
上のプログラムでは,置換制約 `f` を重み1000のペナルティ式として目的関数に加えました.
PyQBPP では,行と列の one-hot 制約を **`qbpp.cons()`** で囲んで**制約であることを明示**できます:
```{literalinclude} /../programFiles/pythonPrograms/example/math/permutation-program4.py
:language: python
:caption: permutation-program4.py
```
配列の比較 `qbpp.vector_sum(x, axis=1) == 1` を `qbpp.cons()` で囲むと,
**要素ごと**（ここでは行ごと・列ごと）**に1本の制約**として宣言されます.
宣言された制約は**制約として特別に処理**され,バンドルされたソルバーは
制約を満たす置換行列を効率よく探索します.
`h.cons(sol)` は解 `sol` で違反している制約の本数を返します（0 なら `x` は置換行列）.
このプログラムの出力は以下のとおりです:
```{include} /../programFiles/markDown/example/math/permutation.md
:start-after:
:end-before:
```
`cons()` の詳しい使い方は[ネイティブ制約](../basic/constraints.md)をご覧ください.
:::