組み合わせ最適化問題において、代表的な解法の概念を説明することができる #1 近似解法

# データサイエンティスト スキルチェックリストver.3

『[データサイエンティスト スキルチェックリストver.3](https://www.datascientist.or.jp/common/docs/skillcheck_ver3.00.pdf)』はデータサイエンティストに必要とされるスキルをチェックリスト化したもので、データサイエンティストに求められるスキルを「ビジネス力」「データサイエンス力」「データエンジニアリング力」の３分野に分けて定義したものです。

本記事では、データサイエンス力のチェックリストNo.269にある項目について、学んだことをまとめてみます。

| スキルカテゴリ | スキルレベル | サブカテゴリ | チェック項目                                                                                                                                                                                                                   |
|----------------|--------------|--------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| 最適化         | ★★           | 最適化       | 組み合わせ最適化問題において、代表的な解法の概念を説明することができる（厳密解法（分枝限定法、動的計画法、切除平面法）、近似解法（局所探索、貪欲法など）、メタヒューリスティック解法（遺伝的アルゴリズム、タブーサーチなど）） |

特に本記事は、この中でも近似解法について記述していきます。


# 組み合わせ最適化

参考：https://annealing-cloud.com/knowledge/1.html
> 組合せ最適化問題とは、様々な制約の下で多くの選択肢の中から、ある指標（価値）を最も良くする変数の値（組合せ）を求めることです。

## 問題

### ナップサック問題

参考：https://ja.wikipedia.org/wiki/%E3%83%8A%E3%83%83%E3%83%97%E3%82%B5%E3%83%83%E3%82%AF%E5%95%8F%E9%A1%8C
> 「容量 C のナップサックが一つと、n 種類の品物（各々、価値 pi, 容積 ci）が与えられたとき、ナップサックの容量 C を超えない範囲でいくつかの品物をナップサックに詰め、ナップサックに入れた品物の価値の和を最大化するにはどの品物を選べばよいか」という整数計画問題である。

例えば、ナップサックの重量制限が20[kg]の下で、以下のような表があるとします。

|商品|1|2|3|4|5|6|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|重さ[kg]|10|1|2|5|15|7|
|価値[万円]|20|2|5|10|25|15|


#### 定式化

今、ナップサックの重量制限が$C$[kg]で、商品$n$個に、重さ$c_i$と価値$p_i$があるとします。
そうすると、ナップサック問題は以下のように定義できます。

```math
\begin{align}
\max &\ f(x) = \sum_{i=1}^n c_i p_i \\
s.t. &\ \sum_{i=1}^n c_i x_i \leq C \\
&\ x_i \in \{0,1\}
\end{align}
```
$x_i$はその商品を入れるか(1)入れないか(0)を表した変数です。
上であげた例題を定式化してみると、以下のようになります。

```math
\begin{align}
\max &\ \ 20x_1 + 2x_2 + 5x_3 + 10x_4 + 25x_5 + 15x_6 \\
s.t. &\ \ 10x_1 + x_2 + 2x_3 + 5x_4 + 15x_5 + 7x_6 \leq 20 \\
&\ \ x_i \in \{0,1\} \ \ (i=1,\dots,6)
\end{align}
```

### 巡回セールスマン問題

参考：https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C
> 都市の集合と各2都市間の移動コスト（たとえば距離）が与えられたとき、全ての都市をちょうど一度ずつ巡り出発地に戻る巡回路のうちで総移動コストが最小のものを求める（セールスマンが所定の複数の都市を1回だけ巡回する場合の最短経路を求める）組合せ最適化問題である。
![Branchbound.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/531226/d8ec426d-5232-8019-697a-e5b3a03fa9a1.gif)
※巡回セールスマン問題を総当たりで解く場合のイメージ。左側で一つずつ探していき、より効率のいいルートが見つかった場合、右側のグラフが更新される。

#### 定式化

都市$i$から都市$j$までの距離を$c_{ij}$とし、都市$i$から$j$にセールスマンが行く場合$1$(それ以外は0)となる変数$x_{ij}$を用いると、以下のようになります。

```math
\begin{align}
\min &\ \sum_{i=1}^n \sum_{j=1}^n c_{ij} x_{ij}\\
s.t. &\ x_{ij} \in \{0,1\} \\
&\ \sum_{j:j \neq i} x_{ij} = 1 \ (i=1,\dots,n)\\
&\ \sum_{i:i \neq j} x_{ij} = 1 \ (j=1,\dots,n)\\
&\ \sum_{i\in S} \sum_{j\in S} x_{ij} \leq |S| - 1 \ (S \subset \{1,\dots,n\})
\end{align}
```

制約条件の2つ目と3つ目は都市$i$を出発する回数と都市$j$に到着する回数がそれぞれ1回ずつであることを示しています。
また、制約条件の4つ目は部分巡回路がないことを意味しています。部分巡回路は、例えば7つの都市がある場合に、4つの都市と3つの都市で、それぞれ四角形、三角形を作るような巡回路です。この問題では、7つの都市ならば、全てがつながって七角形を作ることが条件となっているため、この条件が加わっています。


# 解法
## 解法の分類

まずは、算出される最適解の厳密性で定義される分類として、厳密解法と近似解法があります。

* 厳密解法
    * 大域的最適解を保証する
    * 一般に計算量が大きくなる傾向にある
* 近似解法
    * 大域的最適解は保証しない
    * 良い解を高速に算出する

次に、問題の汎化性に関するものとして、ヒューリスティック解法とメタヒューリスティック解法があります。

* ヒューリスティック解法
    * 発見的手法と言われ、そこそこの近似解を算出する方法
    * 一般的にある問題に特化した方法
* メタヒューリスティック解法
    * ヒューリスティック解法の内、特定の問題に依存せず手法のみが独立した汎用的な近似アルゴリズム
    * 特定の問題のヒューリスティック解法に精度は劣ることが多いが、その分、汎化性を得ている

本記事では、近似解法について紹介していきます。

* 貪欲法
* 緩和法

(局所探索法は#3のメタヒューリスティック解法で記述予定)

## 貪欲法

目的関数への貢献度を示す局所的な評価に基づいて、ステップごとに、その時点で最も得をする選択をするアルゴリズムです。基本的に、試行錯誤は含まないものとなっています。
また、貪欲法で得られる解は、最大化(最小化)問題なら下界(上界)になっています。ここでは『最適解か、それよりも悪い解』であるということを意味しています。


### ナップサック問題の場合

1. 各商品の評価値を決定。価値/重さがよく使用されます
2. 評価値が大きい順に並び替えます
3. $i$番目の商品を入れても容量を超えないなら入れます
4. 手順3を繰り返します

例えば、上述した例題で実行すると

|商品|1|2|3|4|5|6|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|重さ[kg]|10|1|2|5|15|7|
|価値[万円]|20|2|5|10|25|15|
|効率|2|2|2.5|2|1.67|2.15|

実行手順を示していくと以下のようになります。

* 商品3を入れます(残り18[kg]、価値5[万円])
* 商品6を入れます(残り11[kg]、価値20[万円])
* 商品1と商品2が入ります(残り0[kg]、価値42[万円])

今回の問題では、最適解と一致します。

### 巡回セールスマン問題の場合

1. 初期設定
    * 最初の都市(現在地)を選ぶ
    * 距離の総和を0とする
2. 以下を繰り返す
    * まだ訪れていない都市の中で、現在地から最も近い都市を選ぶ
    * 最も近い都市までの距離を加算する
    * 現在地を最も近い都市に更新する


## 緩和法

難しい問題の制約条件を緩めて、解き易い問題に変換し、変換後の問題の解から元の問題の情報を得る解法のことを言います。
制約を緩めているため、ここで得られる緩和解は、最大化(最小化)問題なら、上界(下界)になっています。ここでは『最適解か、それよりも良い解』であるということを意味しています。

### ナップサック問題の場合

ナップサック問題の場合、線形計画緩和を行います。
線形計画緩和問題とは、整数計画問題の整数制約を取り除いてできる問題のことを指し、線形最適化問題になります。

例えば、上述した、ナップサック問題の例題は以下の定式化がされていました。

```math
\begin{align}
\max &\ \ 20x_1 + 2x_2 + 5x_3 + 10x_4 + 25x_5 + 15x_6 \\
s.t. &\ \ 10x_1 + x_2 + 2x_3 + 5x_4 + 15x_5 + 7x_6 \leq 20 \\
&\ \ x_i \in \{0,1\} \ \ (i=1,\dots,6)
\end{align}
```

この問題を、線形緩和すると、以下のような問題になります。

```math
\begin{align}
\max &\ \ 20x_1 + 2x_2 + 5x_3 + 10x_4 + 25x_5 + 15x_6 \\
s.t. &\ \ 10x_1 + x_2 + 2x_3 + 5x_4 + 15x_5 + 7x_6 \leq 20 \\
&\ \ 0 \leq x_i \leq 1 \ \ (i=1,\dots,6)\\
&\ \ x_i \in \mathbb{R}\ \ (i=1,\dots,6)
\end{align}
```

#### Pythonによる解

```python
# ライブラリ
from pulp import *
# モデルの定義
problem = LpProblem(sense=LpMaximize)
# 変数の定義
x1,x2,x3,x4,x5,x6 = [ LpVariable('x'+str(i), 0, 1) for i in [1,2,3,4,5,6] ]
# 目的関数の定義
problem += lpDot([20,2,5,10,25,15], [x1,x2,x3,x4,x5,x6])
# 制約条件の定義
problem += lpDot([10,1,2,5,15,7], [x1,x2,x3,x4,x5,x6]) <= 20

# 解く
status = problem.solve()
# 出力
print("目的関数値:")
print(value(problem.objective))
print("それぞれの値")
print([value(v) for v in [x1,x2,x3,x4,x5,x6]])
```

$(x_1,x_2,x_3,x_4,x_5,x_6)=(0.5,1.0,1.0,1.0,0.0,1.0)$のときに、重さ$19$[kg]、価値が$42$[万円]となります。
この時、$x_1=0.5$のように、整数計画問題としては成立しない解が出てきますが、線形緩和問題としては成立しています。

なお、元の整数計画問題について同様に解いてみると、

```python
# ライブラリ
from pulp import *
# 整数計画の定義
problem = LpProblem(sense=LpMaximize)
# 変数の定義
x1,x2,x3,x4,x5,x6 = [ LpVariable('x'+str(i), cat=LpBinary) for i in [1,2,3,4,5,6] ]
# 目的関数の定義
problem += lpDot([20,2,5,10,25,15], [x1,x2,x3,x4,x5,x6])
# 制約条件の定義
problem += lpDot([10,1,2,5,15,7], [x1,x2,x3,x4,x5,x6]) <= 20

# 解く
status = problem.solve()
# 出力
print("目的関数値:")
print(value(problem.objective))
print("それぞれの値")
print([value(v) for v in [x1,x2,x3,x4,x5,x6]])
```
$(x_1,x_2,x_3,x_4,x_5,x_6)=(1,1,1,0,0,1)$のときに、重さ$20$[kg]、価値が$42$[万円]となります。
今回の例題では、最適解における各変数値は違えど、最適解での価値は一緒になるという結果になりました。


### 巡回セールスマン問題の場合

巡回セールスマン問題の場合、2段階の緩和が考えられます。
その前に、元の定式化を再掲しておきます。

```math
\begin{align}
\min &\ \sum_{i=1}^n \sum_{j=1}^n c_{ij} x_{ij}\\
s.t. &\ x_{ij} \in \{0,1\} \\
&\ \sum_{j:j \neq i} x_{ij} = 1 \ (i=1,\dots,n)\\
&\ \sum_{i:i \neq j} x_{ij} = 1 \ (j=1,\dots,n)\\
&\ \sum_{i\in S} \sum_{j\in S} x_{ij} \leq |S| - 1 \ (S \subset \{1,\dots,n\})
\end{align}
```

制約条件の2つ目と3つ目はある都市への訪問回数の制限、制約条件の4つ目は部分巡回路がない制限を表していました。

まず、1個目の緩和では部分巡回路除去の制約式を除きます。
これにより、いわゆる割当問題になります。

```math
\begin{align}
\min &\ \sum_{i=1}^n \sum_{j=1}^n c_{ij} x_{ij}\\
s.t. &\ x_{ij} \in \{0,1\} \\
&\ \sum_{j:j \neq i} x_{ij} = 1 \ (i=1,\dots,n)\\
&\ \sum_{i:i \neq j} x_{ij} = 1 \ (j=1,\dots,n)
\end{align}
```

2個目の緩和は、ナップサック問題でも用いた線形緩和になります。

```math
\begin{align}
\min &\ \sum_{i=1}^n \sum_{j=1}^n c_{ij} x_{ij}\\
s.t. &\ 0 \leq x_{ij} \leq 1 \\
&\ \sum_{j:j \neq i} x_{ij} = 1 \ (i=1,\dots,n)\\
&\ \sum_{i:i \neq j} x_{ij} = 1 \ (j=1,\dots,n)
\end{align}
```

#### Pythonによる解

今回の地点ごとの距離行列$c_{ij}$は以下のように定義します。

|地点|0|1|2|3|
|:-:|:-:|:-:|:-:|:-:|
|0|0.0|3.0|4.0|8.5|
|1|3.0|0.0|5.0|5.0|
|2|4.0|5.0|0.0|6.0|
|3|8.5|5.0|6.0|0.0|

この線形緩和問題は以下のように解きます。

```python
from pulp import *
import numpy as np
#  設定の定義
N = 4
C = np.array([
               [0.0, 3.0, 4.0, 8.5],
               [3.0, 0.0, 5.0, 5.0],
               [4.0, 5.0, 0.0, 6.0],
               [8.5, 5.0, 6.0, 0.0],
               ])
# 最適化モデルの定義
problem = LpProblem(sense=LpMinimize)
# 変数の定義
X = []
for i in range(N):
  x = [ LpVariable('x{}_{}'.format(i, j), 0,1) for j in range(N) ]
  X.append(x)
# 目的関数の定義
problem += lpSum([lpDot(c, x) for c, x in zip(C, X)])
# 制約条件の定義
for i in range(N):
  problem += lpSum([X[j][i] for j in range(N)]) == 1
  problem += lpSum([X[i][j] for j in range(N)]) == 1
  # 自分自身へのルートは0に固定 (制約式にはない暗黙ルール)
  problem += X[i][i] == 0
# 最適化の実行
status = problem.solve()
print("距離:" + problem.objective.value())
```

これを実行すると、緩和問題の最小値(距離)は18であるとわかります。
なお、今回の問題では、簡単な問題のためか、元の問題の最小値と一致します。

## まとめ

今回は、ナップサック問題と巡回セールスマン問題を題材に、組み合わせ最適化問題の近似解法をいくつか紹介しました。(ただの勉強メモですが)
近いうちに、厳密解法、メタヒューリスティック解法など、他の解法についても同様にまとめていこうと思います。