本記事は、当社オウンドメディア「Doors」に移転しました。

約5秒後に自動的にリダイレクトします。

ブレインパッドの社員が「数理最適化技術」に関して連載するこの企画。
第2回は、当社のデータサイエンティストが、有名問題「ナップサック問題」の様々な解法を紹介しながら、実際に筆者が設定した問題例を解く過程を紹介しています。

こんにちは。アナリティクスサービス部の内池です。この記事では【連載】ブレインパッドの数理最適化ブログの第2弾として、最適化手法入門データサイエンス入門シリーズを読んで学んだことを活かし、現実の問題を様々な方法で解いていきたいと思います。

今回のテーマは組合せ最適化の有名問題であるナップサック問題です。ナップサック問題といえば「動的計画法」を思い浮かべる方が多いと思いますが、動的計画法だけでも様々なバリエーションがある他、動的計画法以外の効率的な解法も存在します。これらについて、実際に問題を解きながら追っていきましょう。

ナップサック問題とは
- 注意
様々な解法
- 厳密解法
  - 動的計画法
  - 分枝限定法
    - 分枝限定法の処理ステップ
- 近似解法
解いてみた
おまけ：処理速度の比較
まとめ
参考文献

ナップサック問題とは

ナップサック問題とは、

「容量 C のナップサックが一つと、n 種類の品物（各々、価値 pi, 容積 ci）が与えられたとき、ナップサックの容量 C を超えない範囲でいくつかの品物をナップサックに詰め、ナップサックに入れた品物の価値の和を最大化するにはどの品物を選べばよいか」

という整数計画問題です。（wikipediaより引用）

噛み砕いて言うと、「大きさの決まったナップサックに、どの品物を選んで詰めると価値の合計が最大になるか？」となります。ナップサック問題は整数計画問題（解ベクトル $\boldsymbol{x}$ の要素に整数制約がある）の中でも0-1整数計画問題という問題に分類され、以下のように解が「1か0か」となります。*1

$x_i = \begin{cases}1 \cdots品物を詰める\\0 \cdots品物を詰めない\end{cases}$

以下、各品物が「重さ」を持っており、ナップサックに「合計の重さ（＝容量）」の制限がある場合を見ていきましょう。品物の個数を $N$ 個, 品物 $i$ の重さを $w_i$ kg, 価値を $v_i$ 円, ナップサックの容量を $W$ kgとすると、ナップサック問題は以下のように定式化できます。

$\begin{aligned} \mathrm{Maximize}&\quad \sum_{i=1}^{N}v_ix_i\\ \mathrm{subject\ to}&\quad \sum_{i=1}^{N}w_ix_i\leq W\\ &\quad x_i\in\{0, 1\}\ (i=1, \cdots, N)\\ \end{aligned}$

通常の線形計画問題であればこの問題はGreedy（貪欲）に求解できますが、ナップサック問題では解が離散量となるため、容易に解くことができません。直感的にも、重さあたりの価値が高い順に品物を選ぶ戦略では必ずしも最適解を得られないことがわかります。このように通常の線形計画問題と比べると求解の難易度が非常に高く、多項式時間で最適解を求めることも難しいため、ナップサック問題はNP困難という問題のクラスに属しています。*2

注意

以降、本記事では多くのプログラムが登場します。プログラム中では最初に以下を実行し、品物の個数・ナップサックの容量、各品物の重さ・価値を1行ずつ標準入力から受け取ります。

import sys
input = sys.stdin.readline

# 品物の個数とナップサックの容量を入力
N, W = map(int, input().split())

# 品物の重さと価値を入力
wv = [list(map(int, input().split())) for _ in range(N)]

例えば品物の個数が2個, ナップサックの容量が10kg, 品物1の重さと価値が3と8, 品物2の重さと価値が7と5の場合、入力後の各変数をprintすると以下となります。

print(N, W)
>>>
2 10

print(wv)
>>>
[[3, 8], [7, 5]]

様々な解法

それでは、ナップサック問題の様々な解法を見ていきましょう。

厳密解法

まずは、全組合せの中で最も良い解（厳密解）を求める方法を見ていきましょう。

ナップサック問題はナイーブに実装すると、全組合せを試す必要があるため、計算量は $O(2^N)$ となってしまいます。この計算量は指数時間であり、多くの場合に現実的ではありません。そのため効率的に厳密解を求めるための様々なアイデアが提案されています。

ここからは、ナップサック問題における厳密解法の代表的なアイデアを2つ見ていきましょう。

動的計画法

上で述べた通り、厳密解法をナイーブに実装すると非現実的な計算量となってしまいます。ここで、品物の「重さ」に注目し、以下のように元の問題を小問題に分割することを考えてみましょう。

$dp[i][j]$ = $i$ 番目までの品物から重さの合計が $j$ 以下になるように選んだ場合の合計価値の最大値

dpは動的計画法（dynamic programming）の略です。このように小問題に分割し、dpテーブルを逐次的に更新することで、全組合せを試すことなく効率的に元の問題の厳密解を得ることができます。詳しい説明は他のブログ等に譲りますが、以下のように実装すると計算量が $O(NW)$ となり、 $O(2^N)$ の実装と比べると比較的速く求解することができます。

import sys
input = sys.stdin.readline

# 品物の個数とナップサックの容量を入力
N, W = map(int, input().split())

# 品物の重さと価値を入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 最大化問題なのでdpテーブルを0で初期化
# dp[i][j]は、重量の合計j以下という制約下でi番目までの品物から任意の品物を選んだ場合の価値の最大値
dp = [[0] * (W + 1) for _ in range(N + 1)]

# 品物i+1を選ばない場合で仮更新してから、ナップサックにi+1個目を詰め込める場合は再更新
for i, (w, v) in enumerate(wv):
    for j in range(W + 1):
        dp[i + 1][j] = dp[i][j]
        if j - w >= 0:
            dp[i + 1][j] = max(dp[i + 1][j], dp[i][j - w] + v)

print(dp[-1][-1])

しかし、この解法ではナップサックの容量が非常に大きい場合に現実的な時間での求解が難しくなります。そこで、今度は品物の「価値」に注目してみましょう。

$dp[i][j]$ = $i$ 番目までの品物から任意の品物を選ぶ場合、最小いくつの重さで価値の合計 $j$ 以上となるか

このように「価値」に注目して考えると、品物の価値の合計を $V$ とすると計算量は $O(NV)$ となり、 $V$ が小さい場合に、重さに注目するよりも効率的に求解することができます。実装は以下となります。

import sys
input = sys.stdin.readline

# dpテーブルから特定の重さのindexを逆引きするため二分探索ライブラリをimport
import bisect

# 品物の個数とナップサックの容量を入力
N, W = map(int, input().split())

# 品物の重さと価値を入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 品物を全て選んだ場合の価値の合計を計算
V = sum(x[1] for x in wv)

# 最小化問題なのでdpテーブルをINFで初期化
# 1 << 60は2**60, pow(2, 60)と同義
INF = 1 << 60

# i行目への遷移はi-1行目のみから行われるため、工夫するとdpテーブルは1行で済む
dp = [INF] * (V + 1)
dp[0] = 0

# 価値jを達成するための最小の重さ = min(dp[j], dp[max(0, j - v)] + w)
for w, v in wv:
    dp_tmp = dp[:]
    for j in range(1, V + 1):
        # j - vが負の場合は0からの遷移としてよい
        idx_prev = max(0, j - v)
        dp_tmp[j] = min(dp[j], dp[idx_prev] + w)
    dp = dp_tmp

# 重さWで達成できる最大価値を求める
print(bisect.bisect_right(dp, W) - 1)

分枝限定法

ナップサック問題を効率的に解くためのアイデアは動的計画法だけではありません。元の整数計画問題を線形計画問題に「緩和」して解く、連続緩和問題の考え方を導入して解く方法があります。それが分枝限定法です。

分枝限定法では、元の整数計画問題の制約を以下のように変更（緩和）して求解を繰り返します。

$\begin{aligned} \mathrm{Maximize}&\quad\boldsymbol v^T\boldsymbol x\\ \mathrm{subject\ to}&\quad\boldsymbol w^T\boldsymbol x\leq W\\ &\quad0\leq\boldsymbol x\leq1 \end{aligned}$

元のナップサック問題の制約と比較すると、解の整数制約が「緩和」され、連続量をとれるようになっています。緩和後の問題はシンプルな線形計画問題のため、高速かつ容易に求解することができます。

また、分枝限定法の肝は「分枝操作」と「枝刈り」です。最適な組合せを探索する過程で、線形計画問題（緩和問題）の最適値が必ず整数計画問題の最適値以上となる性質を利用し、比較的容易に解ける線形計画問題の解が暫定値以下の場合、その枝を「刈る」操作をすることで、探索する組合せの数を減らすというのが分枝限定法の基本的な考え方です。

以下、分枝限定法の手続きを見ていきましょう。

分枝限定法の処理ステップ

品物を重さあたりの価値で降順ソートする。
整数計画問題を貪欲に解き、暫定値を更新する。
線形計画問題（緩和問題）を貪欲に解く。この時、線形計画問題の最適値が暫定値以下の場合、「枝刈り」を行う。（その分岐は探索対象から除外する） $\because$ 整数計画問題の最適値 $\leq$ 線形計画問題の最適値
線形計画問題の解ベクトル $\boldsymbol x$ の要素の1つが整数でない場合、その変数を0と1で固定した状態の問題に「分枝」する。これを「分枝操作」という。
以降、2〜4の処理を繰り返す。

試しに、5変数に絞って処理を追ってみましょう。

（例題）※重さあたりの価値でソート済

$\begin{aligned} \mathrm{Maximize}&\quad 40x_1+12x_2+23x_3+30x_4+25x_5\\ \mathrm{subject\ to}&\quad 3x_1+x_2+2x_3+4x_4+5x_5\leq5 \end{aligned}$

まずは整数計画問題を貪欲に解き、暫定値52を得る。

$\begin{aligned} &\boldsymbol x_0^*=(1,1,0,0,0),\ \boldsymbol v^T\boldsymbol x_0^*=52 \end{aligned}$

次に、線形計画問題（緩和問題）を解く。

$\begin{aligned} &\boldsymbol x_1^*=(1,1,\frac12,0,0),\ \boldsymbol v^T\boldsymbol x_1^*=63.5\\ \end{aligned}$

線形計画問題の最適値63.5を得る。これは暫定値の52を上回るため、 $x_3=1$ の場合と $x_3=0$ の場合に「分枝」し、それぞれの問題を解く。

（ $x_3=1$ の場合）

$\begin{aligned} \mathrm{Maximize}&\quad 40x_1+12x_2+23+30x_4+25x_5\\ \mathrm{subject\ to}&\quad 3x_1+x_2+4x_4+5x_5\leq3\\ &\quad 0\leq x_i\leq1 \ (i=1,2,4,5) \end{aligned}$

$\begin{aligned} &\boldsymbol x_2^*=(1,0,1,0,0),\ \boldsymbol v^T\boldsymbol x_2^*=63 \end{aligned}$

線形計画問題の最適値63を得る。これは暫定値の52を上回るため、この枝は「刈らない」。また、線形計画問題の解が整数解（実行可能解）であるため、暫定値を63に更新する。

（ $x_3=0$ の場合）

$\begin{aligned} \mathrm{Maximize}&\quad 40x_1+12x_2+30x_4+25x_5\\ \mathrm{subject\ to}&\quad 3x_1+x_2+4x_4+5x_5\leq5\\ &\quad 0\leq x_i\leq1 \ (i=1,2,4,5) \end{aligned}$

$\begin{aligned} &\boldsymbol x_3^*=(1,1,0,\frac14,0),\ \boldsymbol v^T\boldsymbol x_3^*=59.5 \end{aligned}$

線形計画問題の最適値59.5を得る。これは先程更新した暫定値の63以下のため、この枝は「刈る」。これで全ての枝の探索が終了したため、暫定値の63が最適値となる。処理を終了する。

近似解法

ここまでナップサック問題の厳密解法を紹介してきました。しかし、現実の問題を解く場合、必ずしも現実的な時間で求解できるとは限りません。そのため、より高速に解ける近似解法のアイデアが提案されています。

さて、それではナップサック問題を厳密に解かなくていい場合の戦略について考えてみましょう。すぐに思いつくのは、「重さあたりの価値が大きい順に品物を詰めていき、詰められなくなったら打ち切る」という貪欲な戦略（貪欲法）です。

試しに以下の例題について考えてみましょう。

（例題）

$\begin{aligned} \mathrm{Maximize}&\quad 40x_1+12x_2+23x_3+30x_4+25x_5\\ \mathrm{subject\ to}&\quad 3x_1+x_2+2x_3+4x_4+5x_5\leq5\\ &\quad x_i\in\{0, 1\}\ (i=1, \cdots, 5)\\ \end{aligned}$

上記の貪欲な戦略（貪欲法）で解くと、

\begin{aligned} &\boldsymbol x^*=(1,1,0,0,0),\ \boldsymbol v^T\boldsymbol x^*=52 \end{aligned}

となります。この問題の最適値は63ですので、そこそこいい解を得ることができました。

しかし、以下の場合はどうでしょうか？

（例題）

$\begin{aligned} \mathrm{Maximize}&\quad 14x_1+52x_2+40x_3+8x_4+4x_5\\ \mathrm{subject\ to}&\quad x_1+4x_2+4x_3+4x_4+4x_5\leq4\\ &\quad x_i\in\{0, 1\}\ (i=1, \cdots, 5)\\ \end{aligned}$

貪欲に解くと、

\begin{aligned} &\boldsymbol x^*=(1,0,0,0,0),\ \boldsymbol v^T\boldsymbol x^*=14 \end{aligned}

となり、この問題の最適値である52と比べるとあまり良くない解となってしまいます。これはナップサック問題が整数計画問題であることに起因しています。

そこで、上記の貪欲法を以下のように改良してみましょう。 $x^{G}$ を整数計画問題を貪欲（Greedy）に解いた場合の解, $v^{max}$ をvの最大値とすると、

$\begin{aligned} \boldsymbol x^G&=[1,\cdots,1,0,\cdots,0]\\ v^{max}&=\mathrm{max}\{v_1,\cdots\,v_n\} \end{aligned}$

ここで、貪欲法による解 $x^{G}$ と $v^{max}$ のみ選ぶ場合を比較して良い方を採用する。

上記の例題では、

$\boldsymbol v^T\boldsymbol x^G=14,\ v^{max}=52$

となり、良い方である52を採用することができました。

実は、改良した解法は厳密解の $\displaystyle \frac12$ 以上の精度（ $=$ 近似比 $=$ 近似解 / 厳密解）が得られることが保証されている、非常に強力なアイデアです。

それでは、なぜ近似比が $\displaystyle \frac12$ 以上になるかについて考えてみましょう。

＜証明＞

整数計画問題の最適解を $\boldsymbol x^{opt}$ , 線形計画問題（緩和問題）の最適解を $\boldsymbol x^{LP}$ *3 , 線形計画問題を解いて得られた解の端数分の価値を $v_q$ とすると、

$\begin{aligned} \boldsymbol v^T\boldsymbol x^{opt}&\leq\boldsymbol v^T\boldsymbol x^{LP}\quad(\because整数計画問題の最適値\leq線形計画問題の最適値)\\ &\leq\boldsymbol v^T\boldsymbol x^G+v_q\quad(\because線形計画問題の最適値\leq解\ \boldsymbol x^{G}\ における合計価値+端数)\\ &\leq\boldsymbol v^T\boldsymbol x^G+v^{max}\quad(\because端数\leq v^{max}) \end{aligned}$

よって、

$\displaystyle\frac{\mathrm{max}\{\boldsymbol v^T\boldsymbol x^G,v^{max}\}}{\boldsymbol v^T\boldsymbol x^{opt}}\geq\frac{\mathrm{max}\{\boldsymbol v^T\boldsymbol x^G,v^{max}\}}{\boldsymbol v^T\boldsymbol x^{G}+v^{max}}$

左辺は近似比を表しており、右辺は必ず $\displaystyle\frac12$ 以上となるため、近似比が $\displaystyle\frac12$ となることが示されました。

【参考】UTokyo OCWx 数理手法Ⅲ

解いてみた

それではいよいよ、ナップサック問題を様々な方法で解いてみましょう！

問題

あなたはメイドカフェ「ほぉ〜むしっくカフェ」に遊びに来ています。今日の予算は5,000円で、この範囲内で幸せさを最大化したいです。ほぉ〜むしっくカフェには以下のメニューがあり、各メニューに設定された料金を支払うと、一定の幸せさを得ることができます。但し、各メニューは1個までしかオーダーできません。最大でいくつの幸せさを得ることができるかを求めてください。なお、あなたの胃袋の容量は無限大で、幸せさは足すことができます。

メニュー	金額	幸せさ
メイドさんオリジナルカクテル	750円	7
メイドさんの手作りパフェ	700円	6
メイドさんのお絵かきオムライス	900円	10
メイドさんと記念撮影（チェキ）	600円	10
メイドさんと記念撮影（デカチェキ）	1,200円	22
ドリンクセット（ドリンク+チェキ）	1,250円	17
デザートセット（ドリンク+デザート+チェキ）	1,850円	23
フードセット（ドリンク+フード+チェキ）	2,050円	27
お土産セット（ドリンク+お土産+チェキ）	1,750円	20
フルセット（ドリンク+デザート+フード+チェキ）	2,700円	33
プレミアムセット（ドリンク+デザート+フード+お土産+チェキ）	3,150円	36

この問題は、メニューの個数を $N$ , メニュー $i$ の金額と幸せさを $w_i$ と $v_i$ , 予算を $W$ と解釈すると、ナップサック問題の式をそのまま適用することが可能です。

厳密解法

それでは、まずはこの問題の厳密解を求めてみましょう。

動的計画法1

まずは、上で解説したプログラムをそのまま実行します。

import sys
input = sys.stdin.readline

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 最大化問題なのでdpテーブルを0で初期化
# dp[i][j]は、料金の合計j以下という制約下でi番目までのメニューから任意のメニューを選んだ場合の幸せさの最大値
dp = [[0] * (W + 1) for _ in range(N + 1)]

# メニューi+1を選ばない場合で仮更新してから、メニューi+1個目を選べる場合は再更新
for i, (w, v) in enumerate(wv):
    for j in range(W + 1):
        dp[i + 1][j] = dp[i][j]
        if j - w >= 0:
            dp[i + 1][j] = max(dp[i + 1][j], dp[i][j - w] + v)

print(dp[-1][-1])
>>>
72

実行結果は72となりました。5,000円の予算内で最大72の幸せさを得られることがわかりました。

さて、それでは最適な組合せはどうなるでしょうか？実は、上記のプログラムでは最適な組合せを得ることはできません。そこで、以下のようにプログラムを修正してみましょう。

import sys
input = sys.stdin.readline

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 最大化問題なのでdpテーブルを0で初期化
# dp[i][j]は、料金の合計j以下という制約下でi番目までのメニューから任意のメニューを選んだ場合の幸せさの最大値
dp = [[0] * (W + 1) for _ in range(N + 1)]

# どのメニューを選んだかを管理する
choice = [[""] * (W + 1) for _ in range(N + 1)]

# メニューi+1を選ばない場合で仮更新してから、メニューi+1個目を選べる場合は再更新
for i, (w, v) in enumerate(wv):
    for j in range(W + 1):
        dp[i + 1][j] = dp[i][j]
        # メニューi+1を選ばない場合は遷移前の集合に0を付加する
        choice[i + 1][j] = choice[i][j] + "0"
        if j - w >= 0:
            # メニューi+1を選ぶ場合は遷移前の集合に1を付加する
            if dp[i + 1][j] < dp[i][j - w] + v:
                dp[i + 1][j] = dp[i][j - w] + v
                choice[i + 1][j] = choice[i][j - w] + "1"

print(choice[-1][-1])
>>>
00011110000

このようにプログラムを修正すると、各 $i, j$ の時にどの品物を選ぶ組合せが最適か？を管理することができ、最終的に元の問題の最適な組合せを得ることができます。*4

このプログラムを実行すると $(0,0,0,1,1,1,1,0,0,0,0)$ が得られ、フードやお土産に目移りせず、チェキやドリンクを重視すると幸せになれることがわかりました。

動的計画法2

今回の問題は上記の解法により現実的な時間で解くことができますが、より大きな問題を解く場合はどうすればいいでしょうか。上記のプログラムをよく見ると、実は $i+1$ 番目の品物について考える時、dpテーブルの $i$ 行目からしか遷移していないことがわかります。

ここまでわかれば、「dpテーブルはベクトルでよい」という結論に達することができます。

但し、dpテーブルをベクトルで管理するため、idxの先頭から愚直にループを回していくとテーブル更新処理が重複してしまいます。そこで、ループを「逆回し」することを考えてみましょう。

import sys
input = sys.stdin.readline

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = [list(map(int, input().split())) for _ in range(N)]

# dpテーブルをベクトルで管理し、効率よく処理することを考える
dp = [0] * (W + 1)
choice = [""] * (W + 1)

# 逆回しをしてdpテーブル更新の重複を避ける
for w, v in wv:
    for prev in range(W - w, -1, -1):
        new_v = dp[prev] + v
        if new_v > dp[prev + w]:
            dp[prev + w] = new_v
            choice[prev + w] = choice[prev] + "1"
        else:
            choice[prev + w] = choice[prev + w] + "0"
    # 遷移できない場合はそのメニューを選べないため、0を付加する
    for i in range(w):
        choice[i] += "0"
            
print(dp[-1], choice[-1])
>>>
72 00011110000

逆回しをすることで上手くテーブル更新処理の重複が避けられ、動的計画法1と同じ結果を得ることができました。

動的計画法3

ナップサック問題における動的計画法は、NumPyを使うことでも効率化を図れます。動的計画法2の逆回しの戦略では各idxにforループでアクセスしていましたが、NumPyを使うことでこの処理を一気に行うことができます。

import sys
input = sys.stdin.readline

import numpy as np

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = np.array([list(map(int, input().split())) for _ in range(N)])

# dpテーブルを1行のベクトルで管理し、一気に処理することを考える
dp = np.zeros(W + 1, dtype=int)

# 実は、w個分ずらしたベクトル同士を比較してmaxをとればよい
# np.maximumで書くことでより高速に動作する
for w, v in wv:
    dp[w:] = np.vstack([dp[w:], dp[:-w] + v]).max(axis=0)

print(dp[-1])
>>>
72

NumPyを使って一気に処理をすることで、更新処理の重複を気にすることなく最適値72を得ることができました。

分枝限定法

今度は、もう1つのアイデアである分枝限定法で最適値を求めてみましょう。以下プログラムは「分枝後の小問題を再帰処理で解く」というアプローチを採用し、分枝限定法の探索処理を実現しています。（簡単のため、最適値のみを求めるプログラムとしています）

import sys
input = sys.stdin.readline

# 最大再帰数の引き上げ
sys.setrecursionlimit(10**7)

import numpy as np

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 価格あたりの幸せさ（コスパ）で降順ソート
wv = np.array(sorted(wv, key=lambda x: x[1] / x[0], reverse=True))


# 分枝操作をする関数
def rec(N, W, wv, acc, ans):
    # 分枝操作が終わったら0を返す
    if N <= 0 or W <= 0:
        return 0

    # 整数計画問題を貪欲に解いた場合の値を計算
    idx = np.searchsorted(acc[:, 0], W, side="right") - 1
    ans_ip = acc[idx, 1]

    # 線形計画問題を貪欲に解いた場合の値を計算
    if idx < N:
        ans_lp = ans_ip + (W - acc[idx, 0]) / wv[idx, 0] * wv[idx, 1]
    else:
        ans_lp = ans_ip

    # 線形計画問題の最適値が暫定値以下の場合は分枝操作を終了する
    if ans_lp <= ans:
        return ans
    # そうでない場合は暫定値を更新
    ans = max(ans, ans_ip)

    # 線形計画問題の解が整数計画問題の解と一致しない場合は分枝操作を再度行う
    if ans_lp != ans_ip:
        new_wv = np.vstack([wv[:idx], wv[idx + 1:]])

        # 分枝後の累積和を事前計算
        nxt_acc = acc.copy()
        nxt_acc[idx + 1:] -= wv[idx]
        nxt_acc = np.vstack([nxt_acc[:idx], nxt_acc[idx + 1:]])

        # メニューidxを選ぶ場合
        if W - wv[idx, 0] >= 0:
            ans = max(ans, rec(
                N - 1, W - wv[idx, 0], new_wv, nxt_acc, ans - wv[idx, 1]) + wv[idx, 1])

        # メニューidxを選ばない場合
        if W >= 0:
            ans = max(ans, rec(N - 1, W, new_wv, nxt_acc, ans))

    return ans


acc = np.vstack([np.array([0, 0]), wv]).cumsum(axis=0)

print(rec(N, W, wv, acc, 0))
>>>
72

このプログラムを実行すると、動的計画法と同様に最適値72を得ることができます。ナップサック問題といえば、競技プログラミング等の文脈では動的計画法を用いるのが一般的ですが、このように分枝限定法でも最適値を求めることが可能です。実際、このプログラムは競技プログラミングのオンラインジャッジでも問題なくAcceptされます。

貪欲法

次は、近似比が保証されない貪欲法を試してみましょう。

import sys
input = sys.stdin.readline

import numpy as np

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 価格あたりの幸せさ（コスパ）で降順ソート
wv = np.array(sorted(wv, key=lambda x: x[1] / x[0], reverse=True))

# 整数計画問題を貪欲に解いた場合の値を計算
acc = np.vstack([np.array([0, 0]), wv]).cumsum(axis=0)
idx = np.searchsorted(acc[:, 0], W, side="right") - 1
ans_ip = acc[idx, 1]

# 貪欲に解いた結果を出力
print(ans_ip)
>>>
49

貪欲法による解は上記の分枝限定法の初期解と一致するため、得られる幸せさは49となります。これは最適値の72と比べて著しく悪い結果ではなく、近似比が保証されない貪欲法でもそれなりの結果を得られる場合があることがわかります。

近似解法

今度は近似比が保証された解法を試してみましょう。近似解は、上記の貪欲法による解 $x^{G}$ と $v^{max}$ のみ選ぶ場合のいずれか良い方でした。今回の問題では $v^{max}$ が36のため、得られる幸せさは貪欲法と同じく49となります。この解の近似比は $\displaystyle\frac{49}{72}$ であり、確かに $\displaystyle\frac12$ 以上となっていることが確認できました。また、この近似解は $(0,0,0,1,1,1,0,0,0,0,0)$ であり、やはりフードやお土産に目移りせず、チェキやドリンクを重視すると一定以上の幸せさが保証されることがわかります。このように、必ずしも厳密解を求めなくていい場合は、そこそこ良い解が保証されている解法を採用するのが得策といえます。

import sys
input = sys.stdin.readline

import numpy as np

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = [list(map(int, input().split())) for _ in range(N)]

# 価格あたりの幸せさ（コスパ）で降順ソート
wv = np.array(sorted(wv, key=lambda x: x[1] / x[0], reverse=True))

# 整数計画問題を貪欲に解いた場合の値を計算
acc = np.vstack([np.array([0, 0]), wv]).cumsum(axis=0)
idx = np.searchsorted(acc[:, 0], W, side="right") - 1
ans_ip = acc[idx, 1]

# 貪欲に解いた結果とvmaxの大きい方を出力
print(max(ans_ip, wv[:, 1].max()))
>>>
49

ソルバーによる解法

ナップサック問題は、各種ソルバーを利用することでも容易に求解することができます。試しにPuLPを利用して問題を解いてみましょう。

import sys
input = sys.stdin.readline

import numpy as np
from pulp import LpProblem, LpMaximize, LpVariable, LpBinary, lpDot, value

# メニューの個数と予算を入力
N, W = map(int, input().split())

# 各メニューの料金と幸せさを入力
wv = np.array([list(map(int, input().split())) for _ in range(N)])

obj = LpProblem(sense=LpMaximize)

# 変数の定義
x = [LpVariable("x_{}".format(i), cat=LpBinary) for i in range(N)]

# 制約
obj += lpDot(wv[:, 1], x)
obj += lpDot(wv[:, 0], x) <= W

# 求解
obj.solve()

# 出力
print("最適値 = {}".format(obj.objective.value()))
print("最適解 = {}".format(tuple(int(value(x[i])) for i in range(N))))
>>>
最適値 = 72.0
最適解 = (0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0)

これまでの実装と同様に、最適値72, 最適解 $(0,0,0,1,1,1,1,0,0,0,0)$ を得ることができました。このように、実際に問題を解く際は長年の知見が蓄積されたソルバーを用いるのが一般的です。

おまけ：処理速度の比較

今回の「解いてみた」では、メイドカフェ「ほぉ〜むしっくカフェ」におけるオーダー最適化問題を7種類の解法で解きました。せっかくですので、各解法の処理速度を比較してみました。

＜実行環境＞

MacBook Pro (13-inch, 2017, Two Thunderbolt 3 ports)
2.3 GHz Intel Core i5, 16 GB 2133 MHz LPDDR3
Python 3.8.2
Jupyter Notebook

なお、今回解いた問題は問題の規模が小さいため、以下のプログラムで新たに入力データを作成し、 %%timeit -n 10 -r 10 で処理速度を比較しました。また、同じ条件下で比較するため、最適値のみを求めるプログラム同士の比較としています。

import random

random.seed(3655)

N, W = 1000, 20000
wv = [[random.randint(1, 100), random.randint(1, 100)] for _ in range(N)]

＜比較結果＞

各解法の処理速度は以下となりました。

#	解法	処理速度
1	動的計画法（通常）	12 s ± 24.2 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)
2	動的計画法（逆回し）	3.52 s ± 21.8 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)
3	動的計画法（NumPy）	65.9 ms ± 2.78 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)
4	分枝限定法	96.7 ms ± 5.68 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)
5	貪欲法	3.2 ms ± 270 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)
6	近似解法	3.4 ms ± 451 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)
7	ソルバー	133 ms ± 9.84 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)

動的計画法による解法は、最も遅い #1 を基準にすると、逆回しによる解法で約3倍強、NumPyによる解法で約200倍弱改善されました。NumPyで一気に処理する解法で大幅に効率化できていることがわかります。

また、分枝限定法による解法も動的計画法（NumPy）による解法に匹敵する処理速度となっており、十分実用的であるといえそうです。（同様に、ソルバーによる解法も非常に速いことがわかりました。）

余談ですが、以下のように $W$ が大きい場合は動的計画法（NumPy）による解法を分枝限定法による解法が上回ります。特定の解法が常に優れているのではなく、問題に応じて適切に解法を選択するのが重要であるといえそうです。

N, W = 100, 2000000
wv = [[random.randint(1, 100000), random.randint(1, 100000)] for _ in range(N)]

なお、貪欲法・近似解法は厳密解を犠牲にしている分、圧倒的に速い処理速度となりました。

まとめ

長くなりましたが、今回の記事ではナップサック問題を応用して現実の問題を解く様々な方法についてご紹介させていただきました。「ナップサック問題といえば動的計画法」と思われがちですが、その動的計画法だけでも多様な考え方がある他、動的計画法以外にも様々な解法があることがおわかりいただけたのではないかと思います。

最後に…株式会社プレインパッドでは、今回ご紹介した組合せ最適化問題をはじめ、あらゆる数理最適化問題に真剣に取り組んでいます。現実の問題を数理最適化で解いていく仲間を広く募集しておりますので、興味を持っていただけた方はお気軽にお声掛けください！