# 巡回セールスマン問題
:::{container} prog-cpp
巡回セールスマン問題（TSP）は,すべての頂点をちょうど1回ずつ訪問して出発点に戻る最短巡回路を求める問題です.
頂点は平面上に配置され,巡回路の長さはユークリッド距離で測られるものとします.
以下の図は,9頂点と最適巡回路の例を示しています:
## TSPのQUBO定式化
巡回路は頂点の順列として表現できます.
そこで,TSPの解を符号化するために[置換行列](../math/permutation.md)を使用します.
$X=(x_{i,j})$（$0\leq i,j\leq n-1$）を $n\times n$ のバイナリ値の行列とします.
行列 $X$ は**置換行列**であり,各行と各列にちょうど1つの1が含まれます.以下に例を示します.
$x_{k,i}$ を「巡回路の $k$ 番目の位置が頂点 $i$ である」と解釈します.
したがって,$X$ のすべての行とすべての列はone-hotでなければなりません.すなわち以下の制約が成り立つ必要があります:
$$
\begin{aligned}
{\rm row}:& \sum_{j=0}^{n-1}x_{i,j}=1 & (0\leq i\leq n-1)\\
{\rm column}:& \sum_{i=0}^{n-1}x_{i,j}=1 & (0\leq j\leq n-1)
\end{aligned}
$$
$d_{i,j}$ を頂点 $i$ と $j$ の間の距離とします.
置換行列 $X$ に対する巡回路の長さは次のように書けます:
$$
\begin{aligned}
{\rm objective}: &\sum_{k=0}^{n-1}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} d_{i,j}x_{k,i}x_{(k+1)\bmod n,j}
\end{aligned}
$$
この式は,頂点 $i$ が位置 $k$ で訪問され,頂点 $j$ が次の位置（$(k+1)\bmod n$）で訪問されるときにちょうど $d_{i,j}$ を加算するので,巡回路の総距離に等しくなります.
## TSPのHi-QUBOプログラム
上記の置換行列による定式化を用いて,TSPのHi-QUBOプログラムを以下のように記述できます:
```{literalinclude} /../programFiles/cppPrograms/example/real/traveling-salesman-program1.cpp
:language: cpp
:caption: traveling-salesman-program1.cpp
```
このプログラムでは,頂点 `0` から `8` の座標が `Nodes` オブジェクト `nodes` に格納されています.
バイナリ変数の2次元配列 `x` を作成し,one-hot制約と巡回路長の目的関数を構成します.
これらの項は,制約にペナルティ重み（ここでは `1000`）を付けて加算することで,1つのQUBO式 `f` にまとめられます.実行可能性が優先されます.
次に,1.0秒の制限時間で EasySolver を使って `f` を解きます.
得られた割り当て `sol(x)` は置換行列を形成します.
この行列は `qbpp::onehot_to_int()` を使って整数のリスト（順列）`tour` に変換され,出力されます.
最後に,計算された `tour` が有向グラフとして描画され,ファイル `tsp_solution.svg` に保存されます.
このプログラムは以下の出力を生成します:
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
### 正確な距離（double）で最適化する
上の `dist` 関数は,デフォルトの `coeff_t` が整数型のため,各ユークリッド距離を `std::llround` で最も近い整数に丸めています.丸めずに**正確な**ツアー長を最適化したい場合は,[実数（double）係数](../../basic/variables-and-expressions.md#real-double-coefficients) を使います.ヘッダのインクルード前に `DOUBLE_TYPE` を定義し,`dist` が `double` を返すようにするだけです:
```{literalinclude} /../programFiles/cppPrograms/example/real/traveling-salesman-program2.cpp
:language: cpp
:caption: traveling-salesman-program2.cpp
```
プログラムの残りは変更不要です.ソルバーは同じ最適巡回路を発見し,`sol.energy()` は丸めた距離の和ではなく,**正確なツアー長**を `double`（例: `960.443`）で返します.
## `slice`、`concat`、`einsum` を使った簡潔な目的関数
`objective` を構築する三重 for ループは数式
$\sum_{i,j,k} d_{j,k}\, x_{i,j}\, x_{(i+1) \bmod n, k}$ をそのまま書き下したものですが,
2 つの補助配列を用意すれば
[`qbpp::einsum`](../../advanced/Einstein-sum.md) の 1 行で書き換えられます.
1. 距離行列 `d`（形状 $n \times n$,`d[j][k] = nodes.dist(j, k)`）
2. `x` を軸 0 で巡回シフトした変数行列 `x_next`
（`x_next[i][k] = x[(i+1) % n][k]`,[slice と concat](../../advanced/slice-connection.md) で構築）
これらを使うと,目的関数は次の 1 行
`qbpp::einsum<0>("jk,ij,ik->", d, x, x_next)` に置き換わり,三重ループ全体が
不要になります.
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
`qbpp::array(nodes.size(), nodes.size())` は 0 で初期化された
2 次元の `Array<2, coeff_t>` を確保し,そこへ実行時に計算した距離値を
代入します.`slice` と `concat` で `x` の先頭行を末尾に移動し,巡回シフト
された行列 `x_next` を作っています.subscript `"jk,ij,ik->"` によって
`d` と `x` の間で `j` を,`d` と `x_next` の間で `k` を,`x` と `x_next`
の間で `i` を共有させ,`i, j, k` すべてを総和してスカラー目的関数を得ます.
ループ版にあった `if (k != j)` の対角項スキップは不要です.
`nodes.dist(j, j)` は常に 0 なので,対角項は自動的に消えます.
得られる QUBO 式の項集合はループ版と完全に同じですが,構築コードは大幅に
短くなり,しかも `einsum` は内部で複数 CPU スレッドを使って並列に式を構築
するため,$n$ が大きい場合は高速化も期待できます.
## 最初の頂点の固定
一般性を失うことなく,頂点0を巡回路の出発点と仮定できます.
TSPの巡回路は巡回シフトに対して不変であるため,出発位置を固定しても最適巡回路長は変わりません.
出発頂点を固定することで,QUBO式中のバイナリ変数の数を削減できます.
具体的には,置換行列において頂点0を位置0に割り当てることを強制します.
そのために,以下のバイナリ変数を固定します:
$$
\begin{aligned}
x_{0,0} &= 1\\
x_{i,0} &= 0& (i\geq 1)\\
x_{0,j} &= 0& (j\geq 1)
\end{aligned}
$$
これらの割り当てにより,頂点0は位置0にのみ現れ,他の頂点は位置0に割り当てられません.
結果として,実効的な問題サイズが削減され,局所探索ベースのソルバーにとってQUBOが解きやすくなります.
## 最初の頂点を固定するHi-QUBOプログラム
Hi-QUBOでは,固定された変数の割り当てを `qbpp::replace()` 関数を使って適用できます:
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
まず,変数の固定割り当てを格納する `qbpp::MapList` オブジェクト `ml` を作成します.
各割り当ては `push_back()` メンバ関数で追加されます.
次に,`qbpp::replace(f, ml)` を呼び出します.これは `ml` で指定された固定値を元のQUBO式 `f` に代入して得られる新しい式を返します.
結果の式は `g` に格納され,簡約化されます.
次に,`g` に対するソルバーを作成して解 `sol` を得ます.
`sol` は縮小された問題に対応するため,`f` に対する `qbpp::Sol` オブジェクトを作成し,ソルバーの出力 `sol` と固定割り当て `ml` の両方を設定します.
結果の `full_sol` は `x` のすべての変数に対する完全な割り当てを格納します.
最後に,`full_sol(x)` で表される置換行列が `qbpp::onehot_to_int()` を使って順列に変換され,出力されます.
このプログラムは頂点0から始まる以下の巡回路を生成します:
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
:::
:::{container} prog-python
巡回セールスマン問題（TSP）は,すべての頂点をちょうど1回ずつ訪問して出発点に戻る最短巡回路を求める問題です.
頂点は平面上に配置され,巡回路の長さはユークリッド距離で測られるものとします.
以下の図は,9頂点と最適巡回路の例を示しています:
## TSPのQUBO定式化
巡回路は頂点の順列として表現できます.
そこで,TSPの解を符号化するために[置換行列](../math/permutation.md)を使用します.
$X=(x_{i,j})$（$0\leq i,j\leq n-1$）を $n\times n$ のバイナリ値の行列とします.
行列 $X$ は**置換行列**であり,各行と各列にちょうど1つの1が含まれます.以下に例を示します.
$x_{k,i}$ を「巡回路の $k$ 番目の位置が頂点 $i$ である」と解釈します.
したがって,$X$ のすべての行とすべての列はone-hotでなければなりません.すなわち以下の制約が成り立つ必要があります:
$$
\begin{aligned}
{\rm row}:& \sum_{j=0}^{n-1}x_{i,j}=1 & (0\leq i\leq n-1)\\
{\rm column}:& \sum_{i=0}^{n-1}x_{i,j}=1 & (0\leq j\leq n-1)
\end{aligned}
$$
$d_{i,j}$ を頂点 $i$ と $j$ の間の距離とします.
置換行列 $X$ に対する巡回路の長さは次のように書けます:
$$
\begin{aligned}
{\rm objective}: &\sum_{k=0}^{n-1}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} d_{i,j}x_{k,i}x_{(k+1)\bmod n,j}
\end{aligned}
$$
この式は,頂点 $i$ が位置 $k$ で訪問され,頂点 $j$ が次の位置（$(k+1)\bmod n$）で訪問されるときにちょうど $d_{i,j}$ を加算するので,巡回路の総距離に等しくなります.
## TSPのPyQBPPプログラム
上記の置換行列による定式化を用いて,TSPのPyQBPPプログラムを以下のように記述できます:
```{literalinclude} /../programFiles/pythonPrograms/example/real/traveling-salesman-program1.py
:language: python
:caption: traveling-salesman-program1.py
```
このプログラムでは,頂点 `0` から `8` の座標がリスト `nodes` に格納されており,ヘルパー関数 `dist(i, j)` が2頂点間の丸めたユークリッド距離を計算します.
バイナリ変数の2次元配列 `x` を作成し,one-hot制約と巡回路長の目的関数を構成します.
これらの項は,制約にペナルティ重み（ここでは `1000`）を付けて加算することで,1つのQUBO式 `f` にまとめられます.実行可能性が優先されます.
次に,1.0秒の制限時間で `EasySolver` を使って `f` を解きます.
得られた割り当て `sol(x)` は置換行列を形成します.
この行列を,各行で値が1のエントリを探すことで整数のリスト（順列）`tour` に変換し,出力します.
このプログラムは以下の出力を生成します:
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
### 正確な距離（double）で最適化する
上の `dist` 関数は,デフォルトの係数型が整数のため,各ユークリッド距離を `round()` で最も近い整数に丸めています.丸めずに**正確な**ツアー長を最適化したい場合は,[実数（double）係数](../../basic/variables-and-expressions.md#real-double-coefficients) を使います.`pyqbpp.d` をインポートし,`round()` を外すだけです:
```{literalinclude} /../programFiles/pythonPrograms/example/real/traveling-salesman-program2.py
:language: python
:caption: traveling-salesman-program2.py
```
プログラムの残りは変更不要です.ソルバーは同じ最適巡回路を発見し,`sol.energy` は丸めた距離の和ではなく,**正確なツアー長**を `float`（例: `960.443`）で返します.
## スライス、`concat`、`einsum` を使った簡潔な目的関数
`objective` を構築する三重 for ループは数式
$\sum_{i,j,k} d_{j,k}\, x_{i,j}\, x_{(i+1) \bmod n, k}$ をそのまま書き下したものですが,
2 つの補助配列を用意すれば
[`qbpp.einsum`](../../advanced/Einstein-sum.md) の 1 行で書き換えられます.
1. 距離行列 `d`（形状 $n \times n$,`d[j, k] = dist(j, k)`）
2. `x` を軸 0 で巡回シフトした変数行列 `x_next`
（`x_next[i, k] = x[(i+1) % n, k]`,[スライスと `concat`](../../advanced/slice-connection.md) で構築）
これらを使うと,目的関数は次の 1 行
`qbpp.einsum("jk,ij,ik->", d, x, x_next)` に置き換わり,三重ループ全体が
不要になります.
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
`qbpp.array(..., shape=(n, n))` は平坦化したリストから $n \times n$ の距離行列を
作成します.`x[1:]` と `x[:1]` で `x` の行をスライスし,`qbpp.concat([...],
axis=0)` で連結することで,`x` を 1 行だけ巡回シフトした `x_next` が得られます.
subscript `"jk,ij,ik->"` によって `d` と `x` の間で `j` を,
`d` と `x_next` の間で `k` を,`x` と `x_next` の間で `i` を共有させ,
`i, j, k` すべてを総和してスカラー目的関数を得ます.
ループ版にあった `if k != j` の対角項スキップは不要です.
`dist(j, j)` は常に 0 なので,対角項は自動的に消えます.
得られる QUBO 式の項集合はループ版と完全に同じですが,構築コードは大幅に
短くなります.さらに $n$ が大きい場合は,einsum 版は縮約全体が C++
バックエンド内でマルチスレッド実行されるため,ループ版で 1 反復ごとに
発生する Python の `ctypes` オーバーヘッドを回避でき,大幅に高速になります.
## 最初の頂点の固定
一般性を失うことなく,頂点0を巡回路の出発点と仮定できます.
TSPの巡回路は巡回シフトに対して不変であるため,出発位置を固定しても最適巡回路長は変わりません.
出発頂点を固定することで,QUBO式中のバイナリ変数の数を削減できます.
具体的には,置換行列において頂点0を位置0に割り当てることを強制します.
そのために,以下のバイナリ変数を固定します:
$$
\begin{aligned}
x_{0,0} &= 1\\
x_{i,0} &= 0& (i\geq 1)\\
x_{0,j} &= 0& (j\geq 1)
\end{aligned}
$$
これらの割り当てにより,頂点0は位置0にのみ現れ,他の頂点は位置0に割り当てられません.
結果として,実効的な問題サイズが削減され,局所探索ベースのソルバーにとってQUBOが解きやすくなります.
## 最初の頂点を固定するPyQBPPプログラム
PyQBPPでは,固定された変数の割り当てを,変数から値へのPython辞書を `qbpp.replace()` 関数に渡すことで適用できます:
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
まず,変数の固定割り当てを格納するPython辞書 `ml` を作成します.
各キーはバイナリ変数,各値はその固定値（`0` または `1`）です.
次に,`qbpp.replace(f, ml)` を呼び出します.これは `ml` で指定された固定値を元のQUBO式 `f` に代入して得られる新しい式を返します.
結果の式は `g` に格納され,簡約化されます.
次に,`g` に対するソルバーを作成して解 `sol` を得ます.
`sol` は縮小された問題に対応するため,`f` に対する `qbpp.Sol` オブジェクトを作成し,ソルバーの出力 `sol` と固定割り当て `ml` の両方を `set(sol, ml)` で設定します.
結果の `full_sol` は `x` のすべての変数に対する完全な割り当てを格納します.
最後に,`full_sol(x)` で表される置換行列を各行の走査により順列に変換し,出力します.
このプログラムは頂点0から始まる以下の巡回路を生成します:
```{include} /../programFiles/markDown/example/real/traveling-salesman.md
:start-after:
:end-before:
```
## matplotlibによる可視化
以下のコードはTSPの解を可視化します.各頂点はラベル付きの点として描画され,巡回路の各辺は赤い有向矢印として表示されます:
```{literalinclude} /../programFiles/pythonPrograms/example/real/traveling-salesman-program3.py
:language: python
:caption: traveling-salesman-program3.py
```
巡回路は,訪問順に頂点を結ぶ赤い有向矢印で表示されます.
:::