講義メモ: 離散最適化 2週(Knapsack)
離散最適化の授業の第2週でのテーマ。典型的なNP困難(NP-Hard)な問題。NP困難な中では簡単な問題で十分に実用的な解を得られる。この問題を題材にさまざまな手法やアプローチについて学んだ。
問題
ナップサック問題は次のような最適化問題で複雑性クラスはNP困難(NP-Hard)になる。
\[
\begin{aligned}
I = \{1,2,3,\dots\}& \\
\text{Objective function}:& \max(\sum_{i \in I}{x_i v_i}) \\
\text{s.t.}& \sum_{i \in I}{x_i w_i} \le C \\
\text{Decision variables}:& x_i \in N (\forall i \in I)
\end{aligned}
\]
決定変数(Decision variables)の条件が \(
x_i \in \{0,1\}\) になった問題を0-1ナップサック問題と呼ぶ。
これら最適化問題に対応した次の決定問題もナップサック問題と呼ぶ。
\[
\begin{aligned}
I = \{1, 2, 3, \dots\}& \\
\text{Constraint}:& \sum_{i \in I}{x_i v_i} \ge V \\
\text{Constraint}& \sum_{i \in I}{x_i w_i} \le C \\
\text{Variables and Domain}:& x_i \in \{0,1\} (\forall i \in I)
\end{aligned}
\]
ここで\(
V=C\)となる問題を部分和問題と呼び、これはNP完全(NP-Complete)に属する。
というわけでナップサック問題はとても難しい(NP困難)。でも完全多項式時間近似スキーム(FPTAS)が存在するので近似アルゴリズムの観点では簡単な問題に含まれる。
スキーナ下巻での補足。ナップサック問題に対する最良のアルゴリズムを選ぶときの観点はいくつかある。
- 全てのアイテムが同コスト、同価値であるか?
- 全てのアイテムで価値密度が同じか?
- コストが小さめの整数か?
- ナップサックは複数か?
1つめは貪欲法で最適解を求められる。2,3番目は動的計画法などで解く。
容量\(
C\)が極めて大きい時、基準値で割りスケールを小さくして近似することで高速に解くことができる。最後のナップサックが複数あるときはビンパッキングとして解く。
アルゴリズム
講義で説明された困難な問題に対するアプローチは2つ。
- 指数関数的に計算量が増えるのであれば指数の底を小さくする
- 最適(Optimal)を諦めて良い解を求める
1つめは最適解を求めるアプローチで、2週目の中で動的計画法と分枝限定法が説明される。
2つめは最適解を諦めるアプローチで、2週目の最後にLSに触れている。(LSの他にMIP,CPにも触れていたがこれらは最適解を求めるはず。)
最適解を求めるアプローチ
NP困難なので全探索では時間がかかりすぎる。これをどう改善するかがポイント。
理論的に問題のサイズに対して指数関数的に計算時間が増大する特性は変えられない。そこで指数の底が小さくなる工夫をする。
ここで紹介されたのは動的計画法と分枝限定法の2つ。
動的計画法
動的計画法は分割統治をボトムアップにおこなうアプローチだ。
部分問題の結果をメモリに残して利用することで探索対象を削減する。
まず問題を分割できて、さらに部分問題の参照関係に循環があっていはいけない。
また部分問題に関するデータを管理するため扱う状態が小さくならないと空間計算量が莫大になり実用的とは言えなくなる。
動的計画法では問題に内在する順序を用いて小さい順に解を構成していく。
ナップサック問題ではアイテム番号\(
i\)の順序を使うと自然に書ける。
アイディアは \(
i\) より小さいアイテムを使って\(
w\)までの重さ以下の最大の価値をdp[i][w]に残していくというもの。これにより次の\(
i+1\)番目のアイテムを含んだ最適解を求めるのは各\(
w\)ごとに使う使わないを選ぶだけで済む。ここでは重さの合計を部分解の状態として管理している。計算量は \(
O(C|I|)\)となる。
このように入力データサイズではなく入力値が多項式に含まれる計算量を擬多項式時間と呼ぶ。
データサイズとしてみると指数関数になるが、入力が非常に小さい場合には扱うことができる。
0-1ナップサック問題に対するPythonを用いた擬似コードを残しておく。
擬似コードなのでModel,Itemのフィールドは察して欲しい。ごめんなさい。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
|
def solve(problem: Model) -> Tuple[int, List[int]]:
dp = [ [0]*(problem.capacity+1) for i in range(problem.item_count+1)]
for idx, item in enumerate(problem.items):
for w in range(problem.capacity+1):
if w-item.weight >= 0:
dp[idx+1][w] = dp[idx][w-item.weight] + item.value
if dp[idx+1][w] < dp[idx][w]:
dp[idx+1][w] = dp[idx][w]
# 選択肢の復元
taken = [0]*problem.item_count
idx, w = problem.item_count, problem.capacity
value = dp[idx][w]
while idx > 0:
if dp[idx][w] == dp[idx-1][w]:
taken[idx-1] = 0
else:
taken[idx-1] = 1
w -= problem.items[idx-1].weight
idx -= 1
return (value, taken)
|
一方で、ナップサックのキャパシティ(\(
C\))が大きいと前述のDPは遅くなる。
そこで見方を変えると価値の合計を使ったDPも考えることもできる。
こっちでは\(
i\)より前のアイテムで\(
v\)以上の価値をだせる最小の\(
C\)をdp[i][v]に残していく。
こちらの計算量は\(
O(\sum_{i \in I}(v_i)|I|)\)となる。問題の特徴によってどちらを適用するか選ぶ必要がある。
価値を元にした方法もPythonの擬似コードを残しておく。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
def max_reachable_value(problem: Model) -> int:
return sum([ i.value for i in problem.items ])
def solve(problem: Model) -> Tuple[int, List[int]]:
max_value = max_reachable_value(problem)
import sys
inf = sys.maxsize
dp = [ [inf]*(max_value+1) for i in range(problem.item_count+1)]
dp[0][0] = 0
for idx, item in enumerate(problem.items):
for v in range(max_value+1):
dp[idx+1][v] = dp[idx][v]
if dp[idx+1][v] > dp[idx][v-item.value]+item.weight:
dp[idx+1][v] = dp[idx][v-item.value] + item.weight
# 選択肢の復元
taken = [0]*problem.item_count
idx, value = problem.item_count, max_value
weight = dp[idx][value]
while weight > problem.capacity:
value -= 1
weight = dp[idx][value]
v = value
while idx > 0:
if dp[idx][v] == dp[idx-1][v]:
taken[idx-1] = 0
else:
taken[idx-1] = 1
v -= problem.items[idx-1].value
idx -= 1
return (value, taken)
|
分枝限定法
一方で分枝限定法は全探索を改善する。探索中に目的関数の期待値を計算してすでに見つかった解より悪い場合に探索を止めることで探索空間を削減する。そのため目的関数の期待値を正しく計算する必要がある。
目的関数が最大化であれば、期待値は実際の値よりも大きくなるように計算する。
目的関数が最小化であれば、期待値は実際の値よりも小さくなるように計算する。
また期待値の計算は探索より効率的に行えないと意味がない。
効率的な期待値の計算は問題を単純化・簡素化することで得られることが多い。
この期待値としてナップサック問題には2つの計算方法が提示された。
- 重量の制約を取り払う
- 決定変数の変域を区間\(
[0,1]\)にする
後者は線形緩和と呼ばれる。
これは全探索を改善する。探索なので探索順序をどうするかの問題が残っている。
- Depth First
- Best First
- Least Discrepancy
の3戦略が紹介されていた。ここでは Depth-Firstでの実装例を書いておく。
DFSなので再帰で書いた。そして探索中の情報を維持するためにオブジェクトにした。
DFSで線形緩和で実装したがバギーな上にDPの方が早かった。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
|
def value(i: Item) -> float:
return float(i.value)
def simple_estimate(items: List[Item], idx: int, capacity: int) -> int:
return sum([i.value for i in items[idx:]])
def linear_relaxation(items: List[Item], idx: int, capacity: int) -> int:
cap = float(capacity)
total, total_weight = 0.0, 0.0
for i in items[idx:]:
if total_weight + i.weight <= capacity:
total += i.value
total_weight += i.weight
else:
rate = float(capacity)/float(i.weight)
total += rate * float(i.value)
break
return math.ceil(total)
def solve(problem: Model,
score_fun: Callable[[Item], float] = value,
estimater: Callable[[List[Item], int, int], int] = linear_relaxation) -> Tuple[int, int, List[int]]:
items = sorted(problem.items, key=score_fun, reverse=True)
value, weight = 0, 0
so = _solver(problem, items, estimater)
value, taken = so.solve()
return (value, taken)
class _solver():
def __init__(self,
problem: Model, items: List[Item], estimate: Callable[[List[Item], int, int], int]):
self.p = problem
self.items = items
self.estimate_score = estimate
self.max_score = -1
self.taken_at_highest_score = None
def solve(self) -> Tuple[int, List[int]]:
taken = [0]*self.p.item_count
self._solve(taken, 0, 0, 0)
taken = [0]*self.p.item_count
for idx, v in enumerate(self.taken_at_highest_score):
if v == 1:
taken[self.items[idx].index] = 1
return self.max_score, taken
def _solve(self, taken: List[int], value: int, weight: int, idx: int):
if idx == len(taken):
if value > self.max_score:
self.max_score = value
self.taken_at_highest_score = copy.copy(taken)
return
rest_capacity = self.p.capacity-weight
expected = value + self.estimate_score(self.items, idx, rest_capacity)
if expected <= self.max_score:
return
cur = self.items[idx]
if cur.weight + weight <= self.p.capacity:
taken[idx] = 1
self._solve(taken, value + cur.value, weight + cur.weight, idx+1)
taken[idx] = 0
self._solve(taken, value, weight, idx+1)
|
最適解を諦めるアプローチ
スキーナ S.S.によると最適解を諦めるアプローチには大きく分けて2つカテゴリがある。
最適解に対する精度保証の有無で呼び分けている。たとえば最適解に対して0.7倍のスコアが保証される計算手法は近似アルゴリズムに分類される。一方で多くの場合で十分な成果を出すものの精度の保証がなされない計算手法をヒューリスティックと呼ぶ。
似たような方法でも問題によりヒューリスティックになるか近似アルゴリズムになるかアルゴリズムになるかが変わる。たとえば、ある指標を元に局所的に最善策を選ぶ貪欲法というアプローチがある。
最小スパニング木問題でエッジの重みを元にした貪欲法を考えるとPrim’sアルゴリズムやKruskal’sアルゴリズムになり最適解が得られる。
一方で、ナップサック問題では価値密度で貪欲法を実装しても、近似率を保証できずにヒューリスティック(近似保証のない近似アルゴリズム)になる。
ただし、上記の貪欲法の最後の1つのアイテムを外し残った最大価値のアイテムを入れて比較するような貪欲法になると近似率2の近似アルゴリズムになるらしい。
貪欲法
ナップサック問題に対して素朴な貪欲法はヒューリスティック解法にすぎない。
それでも実装が単純でバグが入り込みにくく、他のヒューリスティックや近似アルゴリズムを評価するベンチマークになる。Pythonによる擬似コード。
書いてみると、他に比べて短いし単純。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
|
def value_density(i: Item) -> float:
return float(i.value)/float(i.weight)
def solve(problem: Model, score_fun=value_density) -> Tuple[int, List[int]]:
items = sorted(problem.items, key=score_fun, reverse=True)
value = 0
weight = 0
taken = [0]*problem.item_count
for item in items:
if weight + item.weight <= problem.capacity:
taken[item.index] = 1
value += item.value
weight += item.weight
return (value, taken)
|
