動的計画法１ ナップザック問題

引き続き，動的計画法です．今回はフィボナッチ数の問題よりも少し難しくなりますが，動的計画法のありがたさがよくわかる例かと思います．取り上げる問題は，有名な "ナップザック問題" です．

問題は，

重さが wi, 価値が vi であるような n 個の品物がある．
これらの品物から，重さの総和が W を超えないように選んだ時の，価値の総和の最大値を求めなさい．

わかりやすく言うと，手持ちのナップザックが重さWまで入れてOKだとしたときに，入れたものの価値が最大になるような品物の組み合わせを見つけることです．
この問題を図示すると，下記のようになります．（とっても汚い図で恐縮です(ﾟДﾟ)ﾉ）

上記の図では，ナップザックに入れていい重さが10の時に，与えられた4つの品物（□，〇，△，☆）を選択して価値を高めるとすると，（□，〇，△）の組み合わせを選んでやれば重さが10，価値が14となって価値の総和が最大になります．

１．問題の解き方（統治分割法＆貪欲法）

結局はそれぞれの品物をナップザックに”入れる”か”入れない”かの選択になるので，品物がn個あった場合はそれぞれの品物について”入れる”か”入れない”を考えていくと，組み合わせ数は2ⁿになります．

上図では，ナップザックに"入れる"を1に，"入れない"を0として記載しましたが，全部で2⁴=16通りのパターンがあることがわかると思います．この方法で解く場合，下記のような手順になります．

１．組み合わせをもれなく洗い出す．
２．１で出てきた組み合わせのそれぞれに対して，重さと価値を計算する．
３．計算した組み合わせの中から，重さが１０以下で価値が最大となるものを選択する．

この手順を実装すると，下記のようになります．

class Product(object):
  def __init__(self, val, weight):
    self.value = val
    self.weight = weight

comb_cnt = 0

def find_maximum_knapsack_greedy(max_weight, product_list):
  return find_maximum_knapsack_greedy_int(0, 0, max_weight, product_list)

def find_maximum_knapsack_greedy_int(cur_idx, cur_weight, max_weight, plist):

  if (cur_idx >= len(product_list)):

    global comb_cnt
    comb_cnt = comb_cnt + 1

    if (cur_weight > max_weight):
      return -sys.maxsize
    return 0

  return max(find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist), \
    find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)

if __name__ == '__main__':

  print("Knapsack Problem")
  product_list = [Product(4, 2), Product(5, 2), Product(2, 1), Product(8, 3)]
  max_weight = 5
  print(find_maximum_knapsack_greedy(max_weight, product_list))
  print("Combination Cnt : " + str(comb_cnt))

上の実装では，全探索した時の組み合わせ数を計算したかったので，重さの条件チェックを最後までやっていません．この関数を実行すると組み合わせ数がたしかに16となっていることがわかります．

ただ，重さがナップサックの許容量を超えている時点で計算をやめてもいいはずなので，この条件を付け加えるとコードは下記のようになります．

class Product(object):
  def __init__(self, val, weight):
    self.value = val
    self.weight = weight

def find_maximum_knapsack_greedy_int(cur_idx, cur_weight, max_weight, plist):

  if (cur_weight > max_weight):
    return -sys.maxsize
  elif (cur_idx >= len(product_list)):
    global comb_cnt
    comb_cnt = comb_cnt + 1
    return 0

  return max(find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist), \
    find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)

重さを計算の打ち切り条件にしたことで，組み合わせ数が12となりました．この問題を図示すると，下記のような二分木になります．

結局は，「4つある品物の組み合わせを全部洗い出して，重さ条件を満たす組み合わせのものの中で価値が最大のものを探す」ということをやっていることがわかります．ここで，「品物nを選ぶ，選ばない」という選択が上記のコードのどこで行われているかというと，下記の記述になります．

  return max(find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist), \
    find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)

今考えている品物を選ばない：find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist)
今考えている品物を選らぶ　：find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)

となります．この分岐を書くこと自体は動的計画法というよりは統治分割法の話になるのかもしれませんが，この分岐がはっきりすれば動的計画法での実装に一歩近づくものと思われます．

２．問題の解き方（統治分割法＆メモ化再帰）

１の統治分割法＆貪欲法で書いた二分木ですが，もう少し詳しく書いて見ました．それぞれのノードが（i, j）という形で書かれていますが，ここでは下記のように解釈します．（※ここをどう解釈するかが動的計画法のキモです．）

検討すべき品物数がnであるとしたときに，

（i, j）： 許容重量が j [kg] のときの， i+1 ~ n 番目までの品物を組み合わせて作った価値の総和の最大値．

です．下記に具体的な図を描いてみましたが，二分木の上の決定（例えば，品物1を入れる）が許容重量の形でツリーの下側に伝播していく様子がわかります．また，注意深く見ていくと，同じ構造をしている木が複数あることに気づくと思います．
例えば，(2, 3) 「許容重量が 3[kg] の時の，3 ~ 4 番目までの品物を組み合わせて作った価値の総和の最大値」が2回出てきています．この計算自体は何度やっても同じ結果になるので，この結果を”メモ”することで不要な再帰計算を省略することができます．

ここからメモ化の出番です．計算したパスの結果を保存して再利用することで，再帰関数の呼び出しを抑えることが目的になります．

class Product(object):
  def __init__(self, val, weight):
    self.value = val
    self.weight = weight

def find_maximum_kanpsack_rec_memo(max_weight, plist):
  memo = np.zeros((len(plist), max_weight+1))
  memo.fill(-1)
  memo[:, 0] = 0
  return find_maximum_kanpsack_rec_memo_int(0, max_weight, plist, memo)

def find_maximum_kanpsack_rec_memo_int(cur_idx, cur_weight, plist, memo):

  if (cur_weight < 0):
    return -sys.maxsize
  elif (cur_idx == len(plist)):
    return 0
  elif (memo[cur_idx, cur_weight] != -1):
    return memo[cur_idx, cur_weight]
  
  inc_value = find_maximum_kanpsack_rec_memo_int(cur_idx+1, cur_weight-plist[cur_idx].weight, plist, memo) + plist[cur_idx].value
  exc_value = find_maximum_kanpsack_rec_memo_int(cur_idx+1, cur_weight, plist, memo)

  memo[cur_idx, cur_weight] = max(inc_value, exc_value)

  return memo[cur_idx, cur_weight]

if __name__ == '__main__':

  print("Knapsack Problem")
  product_list = [Product(4, 2), Product(5, 2), Product(2, 1), Product(8, 3)]
  max_weight = 5
  print(find_maximum_kanpsack_rec_memo(max_weight, product_list))

３．問題の解き方（動的計画法）

二分木による可視化とメモ化ができれば，あとはこれを再帰を使わずに for ループで実現するための表を作ることは簡単です．2で可視化したツリーを見るとわかるように，末端以外のノードは自分の下に2つの子供を持っています．当該ノードが持つ値は，下記のどちらか大きいほうの値になります．

１．自分を含めた時の価値の最大値
２．自分を含めなかった時の価値の最大値

(1, 3) のノードを例にして書くと，下記のようになります．

(1, 3) = max((2, 1) + 品物2の価値, (2, 3))

つまり，漸化式は下記のようになります．

val[n+1][j] = 0 
val[i][j] = max(val[i+1][j]) (j < w[i]) : 品物の重量が許容重量より大きく，考慮に入れられない時．
val[i][j] = max(val[i+1][j-w[i]]+value[i], val[i+1][j]) : それ以外

といった具合になります．また，ノードの値を求めるには当該ノードの子供の値が確定していないと求まらないので，for ループでテーブルを作成するには末端のノード，つまり4番目のノードの値から求めてあげる必要があります．この場合，テーブルは下記のようになります．
※ちなみに，5番目の品物というものは存在しないので，ダミーみたいなもんです．

これをコードに落とすと．．．

class Product(object):
  def __init__(self, val, weight):
    self.value = val
    self.weight = weight

def find_maximum_knapsack_dp(max_weight, plist):

  val_table = np.zeros((len(plist)+1, max_weight+1))
  val_table[4, :] = 0

  for i in range(len(plist)-1, -1, -1):
    for j in range(max_weight + 1):
      if (plist[i].weight > j):
        val_table[i, j] = val_table[i+1, j]
      else:
        val1 = val_table[i+1, j-plist[i].weight] + plist[i].value
        val2 = val_table[i+1, j]
        val_table[i, j] = max(val1, val2)

  return val_table[0, max_weight]


if __name__ == '__main__':

  print("Knapsack Problem")
  product_list = [Product(4, 2), Product(5, 2), Product(2, 1), Product(8, 3)]
  max_weight = 5
  print(find_maximum_knapsack_dp(max_weight, product_list))

ちなみに，上記のコードだとDPの表を求める時の一階層目のループが逆方向になってます．順方向のほうがすっきり感じる場合は（自分は少なくともそうでした），問題の見方を少し変えれば順方向のループでテーブルが作れます．

この二分木では，(i, j) を下記のように解釈しています．

（i, j）： 許容重量が j [kg] のときの， i 番目までの品物を組み合わせて作った価値の総和の最大値．

こっちのほうが自然ですよね？こう解釈すると，表が順方向に計算できます．

これをコードに落とすと，，，

def find_maximum_knapsack_dp_low2high(max_weight, plist):

  val_table = np.zeros((len(plist)+1, max_weight+1))
  val_table[0, :] = 0

  for i in range(1, len(plist)+1):
    for j in range(max_weight + 1):
      if (plist[i-1].weight > j):
        val_table[i, j] = val_table[i-1, j]
      else:
        val1 = val_table[i-1, j-plist[i-1].weight] + plist[i-1].value
        val2 = val_table[i-1, j]
        val_table[i, j] = max(val1, val2)

  return val_table[len(plist), max_weight]

となります．引き続き，もう少し動的計画法をやっていきます．