今回は Haskell で「動的計画法 (dynamic programming)」というアルゴリズムに挑戦してみましょう。難しそうな名前がついていますが、これに惑わされてはいけません。動的計画法は、大きな問題を小さな問題に分けて、それを一つ一つ解いていくことで大きな問題を解く方法です。
問題によっては、小さな問題に分割していくと同じ小問題が何回も現れる場合があります。この場合、同じ問題を何回も解くよりも、その解を表などに保存しておいて、必要なときにその表から答を求めた方が、効率良く問題を解くことができるはずです。
どうせ小問題を解かなければならないのであれば、はじめから必要になりそうな小問題を解いて表を埋めておいたほうが、プログラムを作りやすい場合もあります。このように、与えられた問題を解くために小問題の表を埋めてしまう、というのが「動的計画法」の基本的な考え方です。
簡単な例題として、組み合わせの数 \({}_n \mathrm{C}_r\) を求めるプログラムを作ってみましょう。\({}_n \mathrm{C}_r\) を求めるには、次の公式を使えば簡単です。
皆さんお馴染みの公式ですね。(1) と (2) の公式を使うと簡単に (高速に) 答えを求めることができます。ただし、(3) の公式をそのままプログラムすると二重再帰になるので、大きな値を求めると時間がかかってしまいます。実際にプログラムを作って確かめてみましょう。
リスト : 組み合わせの数 comb :: Integer -> Integer -> Integer comb n r | n == r || r == 0 = 1 | otherwise = comb (n - 1) r + comb (n - 1) (r - 1)
ghci> :set +s ghci> comb 20 10 184756 (0.35 secs, 133,101,200 bytes) ghci> comb 22 11 705432 (1.17 secs, 507,983,104 bytes) ghci> comb 24 12 2704156 (4.42 secs, 1,947,065,152 bytes) 実行環境 : Ubunts 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
このように、公式 (3) を使うと時間がかかります。公式 (1), (2) を使えば高速に答えを求めることができますが、今回は動的計画法の例題として、あえてこのプログラムの高速化に挑戦してみましょう。
公式からわかるように、\({}_n \mathrm{C}_r\) の値は \({}_{n-1} \mathrm{C}_r\) と \({}_{n-1} \mathrm{C}_{r-1}\) を足したものです。n = 0 から順番に組み合わせの数を求めて表に格納しておけば、n が大きな値でも簡単に求めることができるはずです。プログラムは次のようになります。
リスト : 組み合わせの数 (動的計画法[1])
combDP :: Int -> Int -> IO Integer
combDP n r = do
a <- newArray ((0, 0), (n, n)) 1 :: IO (IOArray (Int,Int) Integer)
mapM_ (\x -> mapM_ (\y -> do b <- readArray a (x - 1, y)
c <- readArray a (x - 1, y - 1)
writeArray a (x, y) (b + c))
[1 .. x-1])
[2 .. n]
readArray a (n, r)
変数 a に値を格納する二次元配列をセットします。組み合わせの数を (n, r) で表すことにすると、関数 mapM_ の中で (0, 0) から順番に、(1, 0), (1, 1), (2, 0), (2, 1), (2, 2) ... と組み合わせの数を求めて配列 a にセットします。ようするに「パスカルの三角形」を作っていくわけです。最後に配列 a から (n, r) の値を求めて返します。
それでは実際に試してみましょう。
ghci> combDP 100 50 100891344545564193334812497256 (0.02 secs, 14,023,912 bytes) ghci> combDP 200 100 90548514656103281165404177077484163874504589675413336841320 (0.05 secs, 56,201,024 bytes)
動的計画法の効果はとても高いですね。なお、表は二次元配列ではなく一次元配列で済ますこともできます。たとえば、(6, 3) を求めてみましょう。次の図を見てください。
0 1 2 3 4 5 6
-------------------------
0 [ 1 1 1 1 1 1 1 ]
1 [ 1 1 1 1 1 1 1 ]
＼|
2 [ 1 2 1 1 1 1 1 ]
＼|＼|
3 [ 1 3 3 1 1 1 1 ]
＼|＼|＼|
4 [ 1 4 6 4 1 1 1 ]
＼|＼|＼|＼
5 [ 1 5 10 10 5 1 1 ]
＼|＼|＼|＼|＼|
6 [ 1 6 15 20 15 6 1 ]
図 : パスカルの三角形
最初にベクタの内容を 1 に初期化します。n = 0, 1 の場合はこのままで大丈夫です。あとは図のように、隣の要素を足し算するだけです。3 番目の要素の値 20 が (6, 3) の値になります。プログラムは次のようになります。
リスト : 組み合わせの数 (動的計画法[2])
combDP' :: Int -> Int -> IO Integer
combDP' n r = do
a <- newArray (0, n) 1 :: IO (IOArray Int Integer)
mapM_ (\x -> mapM_ (\y -> do b <- readArray a y
c <- readArray a (y - 1)
writeArray a y (b + c))
[x-1, x-2 .. 1])
[2 .. n]
readArray a r
配列の値を書き換えていくので、配列の後方から計算していくことに注意してください。前方から計算すると、値がおかしくなります。
もちろん、リストを使ってプログラムすることもできます。パスカルの三角形は、両側がすべて 1 で内側の数はその左上と右上の和になっています。したがって、n 段の三角形を作る場合、一段前の値を記憶しておいて隣同士の要素を足し算すればいいわけです。一段前の値をリストに格納すると、プログラムは次のようになります。
リスト : 組み合わせの数 (動的計画法[3])
combDP'' :: Int -> Int -> Integer
combDP'' n r = iter1 n [1]
where iter1 0 xs = xs !! r
iter1 n xs = iter1 (n - 1) (iter2 (0:xs))
iter2 [_] = [1]
iter2 (x:y:zs) = (x + y) : iter2 (y:zs)
局所関数 iter2 は、リストの隣同士の要素を足した値をリストに格納して返します。この処理は再帰呼び出しを使えば簡単です。リストの先頭要素 x と次の要素 y を足し算し、その値を再帰呼び出しの返り値に追加すればいいわけです。リストの要素がひとつになったら [1] を返します。また、iter2 を呼び出すときはリストの先頭に 0 を追加します。これで、リストの先頭と最後尾を 1 にすることができます。あとは、関数 iter1 で iter2 を n 回呼び出すだけです。
それでは実行してみましょう。
ghci> combDP'' 100 50 100891344545564193334812497256 (0.01 secs, 1,275,680 bytes) ghci> combDP'' 200 100 90548514656103281165404177077484163874504589675413336841320 (0.01 secs, 4,885,928 bytes)
─┬─ 6 : 6
│
├─ 5 ─ 1 : 5 + 1
│
├─ 4 ┬ 2 : 4 + 2
│ │
│ └ 1 ─ 1 : 4 + 1 + 1
│
├─ 3 ┬ 3 : 3 + 3
│ │
│ ├ 2 ─ 1 : 3 + 2 + 1
│ │
│ └ 1 ─ 1 ─ 1 : 3 + 1 + 1 + 1
│
├─ 2 ┬ 2 ┬ 2 : 2 + 2 + 2
│ │ │
│ │ └ 1 ─ 1 : 2 + 2 + 1 + 1
│ │
│ └ 1 ─ 1 ─ 1 ─ 1 : 2 + 1 + 1 + 1 + 1
│
└─ 1 ─ 1 ─ 1 ─ 1 ─ 1 ─ 1 : 1 + 1 + 1 + 1 + 1 + 1
図 : 整数 6 の分割
もうひとつ、数値計算の例を取り上げましょう。整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。
簡単な例を示しましょう。上の図を見てください。6 の場合、分割の仕方は上図のように 11 通りあります。この数を「分割数」といいます。分割の仕方を列挙する場合、整数 n から k 以下の整数を選んでいくと考えてください。まず、6 から 6 を選びます。すると、残りは 0 になるので、これ以上整数を分割することはできません。次に、6 から 5 を選びます。残りは 1 になるので、1 を選ぶしか方法はありません。
次に、4 を選びます。残りは 2 になるので、2 から 2 以下の整数を分割する方法になります。2 から 2 を選ぶと残りは 0 になるので 2 が得られます。1 を選ぶと残りは 1 になるので、1 + 1 が得られます。したがって、4 + 2, 4 + 1 + 1 となります。同様に、6 から 3 を選ぶと、残りは 3 から 3 以下の整数を選ぶ方法になります。
6 から 2 以下の整数を選ぶ方法は、残り 4 から 2 以下の整数を選ぶ方法になり、そこで 2 を選ぶと 2 から 2 以下の整数を選ぶ方法になります。1 を選ぶと 4 から 1 以下の整数を選ぶ方法になりますが、これは 1 通りしかありません。最後に 6 から 1 を選びますが、これも 1 通りしかありません。これらをすべて足し合わせると 11 通りになります。
整数 n を k 以下の整数で分割する総数を求める関数を p(n, k) とすると、p(n, k) は次のように定義することができます。
たとえば、p(6, 6) は次のように計算することができます。
p(6, 6) => p(0, 6) + p(6, 5)
=> 1 + p(1, 5) + p(6, 4)
=> 1 + 1 + p(2, 4) + p(6, 3)
=> 1 + 1 + 2 + 7
=> 11
p(2, 4) => p(-2, 4) + p(2, 3)
=> 0 + p(-1, 3) + p(2, 2)
=> 0 + 0 + p(0, 2) + p(2, 1)
=> 0 + 0 + 1 + 1
=> 2
p(6, 3) => p(3, 3) + p(6, 2)
=> p(0, 3) + p(3, 2) + p(4, 2) + p(6, 1)
=> 1 + p(1, 2) + p(3, 1) + p(2, 2) + p(4, 1) + 1
=> 1 + 1 + 1 + p(0, 2) + p(2, 1) + 1 + 1
=> 1 + 1 + 1 + 1 + 1 + 1 + 1
=> 7
分割数を求める関数 partition_number は、関数 p(n, k) を使うと次のようにプログラムすることができます。
リスト : 分割数
partition_number :: Integer -> Integer
partition_number n = part_num n n
where
part_num 0 _ = 1
part_num 1 _ = 1
part_num _ 1 = 1
part_num n k
| n < 0 || k < 1 = 0
| otherwise = part_num (n - k) k + part_num n (k - 1)
ghci> partition_number 40 37338 (0.13 secs, 45,665,080 bytes) ghci> partition_number 50 204226 (0.67 secs, 241,764,216 bytes) ghci> partition_number 60 966467 (3.10 secs, 1,121,766,480 bytes)
局所関数 part_num は p(n, k) の定義をそのままプログラムしただけです。ただし、このプログラムは二重再帰で何度も同じ値を求めているため実行速度はとても遅くなります。
動的計画法を使うと、大きな値でも高速に計算することができます。次の図を見てください。
k 1 : [1, 1, 1, 1, 1, 1, 1] 2 : [1, 1, 1+1=2, 1+1=2, 2+1=3, 2+1=3, 3+1=4] => [1, 1, 2, 2, 3, 3, 4] 3: [1, 1, 2, 1+2=3, 1+3=4, 2+3=5, 3+4=7] => [1, 1, 2, 3, 4, 5, 7] 4: [1, 1, 2, 3, 1+4=4, 1+5=6, 2+7=9] => [1, 1, 2, 3, 5, 6, 9 5: [1, 1, 2, 3, 5, 1+6=7, 1+9=10] => [1, 1, 2, 3, 5, 7, 10] 6: [1, 1, 2, 3, 5, 7, 10+1=11] => [1, 1, 2, 3, 5, 7, 11]
大きさ n + 1 の配列を用意します。配列の添字が n を表していて、p(n, 1) から順番に値を求めていきます。p(n, 1) の値は 1 ですから、配列の要素は 1 に初期化します。次に、p(n, 2) の値を求めます。定義により p(n, 2) = p(n - 2, 2) + p(n, 1) なので、2 番目以降の要素に n - 2 番目の要素を加算すれば求めることができます。あとは、k の値をひとつずつ増やして同様の計算を行えば p(n, n) の値を求めることができます。
プログラムは次のようになります。
リスト : 分割数 (動的計画法)
partition_number' :: Int -> IO Integer
partition_number' n = do
a <- newArray (0, n) 1 :: IO (IOArray Int Integer)
mapM_ (\k -> mapM_ (\m -> do b <- readArray a m
c <- readArray a (m - k)
writeArray a m (b + c))
[k .. n])
[2 .. n]
readArray a n
説明をそのままプログラムしただけなので、とくに難しいところはないと思います。それでは実際に試してみましょう。
ghci> partition_number' 40 37338 (0.00 secs, 1,849,272 bytes) ghci> partition_number' 50 204226 (0.00 secs, 2,858,152 bytes) ghci> partition_number' 60 966467 (0.01 secs, 4,091,872 bytes) ghci> partition_number' 100 190569292 (0.01 secs, 11,285,040 bytes) ghci> partition_number' 1000 24061467864032622473692149727991 (1.00 secs, 1,148,937,616 bytes)
もちろん、リストを使ってプログラムすることもできます。
リスト : 別解 (動的計画法)
partition_number'' :: Int -> Integer
partition_number'' n = iter 2 (replicate (n + 1) 1)
where
iter k xs
| k > n = last xs
| otherwise = let ys = (take k xs ++ zipWith (+) ys (drop k xs))
in iter (k + 1) ys
実際の処理は局所関数 iter で行います。リストの値を更新するとき、新しいリストを ys とすると、先頭の k 個の要素は xs と同じで、それ以降の要素は xs の i 番目の要素と ys の i - k 番目の要素を加算します。この処理は、ys と (drop k xs) を zipWith で加算することと同じです。ys の先頭要素は take k xs で求められているので、zipWith で ys の要素を先頭から順番に取り出していくことができます。
それでは実行してみましょう。
ghci> partition_number'' 40 37338 (0.00 secs, 272,192 bytes) ghci> partition_number'' 50 204226 (0.00 secs, 379,072 bytes) ghci> partition_number'' 60 966467 (0.00 secs, 507,592 bytes) ghci> partition_number'' 100 190569292 (0.00 secs, 1,247,952 bytes) ghci> partition_number'' 1000 24061467864032622473692149727991 (0.24 secs, 163,995,424 bytes)
最後に、簡単な例題として「ナップザック問題」を取り上げます。ナップザック (knapsack) とは辞書を引いてみると、ランドセルのような背中にせおう四角形の袋や箱のことを意味します。ここでは物を入れる袋と簡単に考えてください。
ここで、ナップザックの中に品物を詰め込むことを考えてみます。一つのナップザックと複数の品物が与えられたとき、袋に詰めた品物の合計金額が最大になるような選び方を求めることが「ナップザック問題」です。ここでは、同じ品物を何個も選んでもいいのですが、ナップザックの大きさをオーバーしてはいけません。
実はこの「ナップザック問題」が「NP 問題」なのです。世の中にはさまざまな問題が山積していますが、スーパーコンピュータを使っても解くのに数億年かかる、というような難問が「NP 問題」です。これは、厳密に解を求めようとすると、全ての場合について総当たりで調べるしか方法がなく、データ数が多くなると時間がかかるため、現実的な時間では解答を出すことができないというものです。品物の詰め方が難問の一つ、といわれてもピンとこないと思いますが、ナップザック問題は品物の種類が増えるにしたがって、その組み合わせ方が爆発的に増えるのです。
ところが、幸いなことに「ナップザック問題」は実用的には解決済みの問題と考えられています。とくに有名なのが「動的計画法」を用いた解法です。ナップザックと品物の大きさを整数値に限定すれば、動的計画法を用いることで厳密解を求めることができるのです。
それでは具体的に、ナップザック問題に動的計画法を適用してみましょう。ナップザックの大きさは 10 で、次の 3 種類の品物を詰め込むことにします。
（Ａ）大きさ ４ ：金額 ６ （Ｂ）大きさ ３ ：金額 ４ （Ｃ）大きさ １ ：金額 １
まず、大きさが 0 から 10 までのナップザックを用意します。これらのナップザックに品物を順番に詰め込んで、その合計金額を配列に格納しておきます。この配列は品物を詰め込んでいない状態 (金額は全て 0) に初期化します。
最初に品物 A を詰め込みます。このとき、小さなナップザックから順番に詰め込んでいきます。
Ａを詰め込む
↓
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│／│／│／│／│／│／│／│／│／│／│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│０│０│０│０│０│０│０│０│０│０│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
↓ ↓
金額[0] + 6 = 6 > 0 金額[4] に 6 をセット
品物 A が入る大きさのナップザックから詰め込みます。品物を入れるときは、それより小さいナップザックには、その時点での最適な値が決定されていると考えます。ナップザック (4) に品物 A を詰め込む場合、A の大きさだけ空いているナップザック (0) の状態に詰め込めば、ナップザック (4) の最適な値を求めることができるはずです。このように、前に計算された値を使うところが動的計画法の特徴なのです。
具体的には、金額[0] に A の金額 6 を足した値を計算し、金額[4] より大きくなれば金額[4] をその値に更新します。もし、金額[4] より小さいのであれば金額[4] は更新しません。つまり、品物Ａはナップザック (4) には詰め込まないのです。ほかの組み合わせの方が正解だというわけです。
詰め込んだ品物を記憶しておくため、もう一つ配列を用意して、そこに追加した品物の種類を格納しておきます。ナップザックの中身全てを記憶しておく必要はありません。この配列を使って、あとからナップザックの中身を求めることができます。
次に、品物 B を詰め込んでいきます。
Ｂを追加する
↓
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│／│／│Ｂ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│０│０│４│６│６│６│６│12│12│12│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
↓ ↓
金額[1] + 4 = 4 < 6 金額[4] は更新しない
まず、ナップザック (3) に B が詰め込まれます。これは品物 A の場合と同じですね。次に、ナップザック (4) に B を詰めようとします。その値を計算すると 4 となり、金額[4] の値 6 より小さいので、B は詰め込みません。ナップザック (5) の場合も同様です。
Ｂを追加する
↓
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│／│／│Ｂ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│０│０│４│６│６│６│６│12│12│12│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
↓ ↓
金額[3] + 4 = 8 < 6
金額[6] を 8 に更新し、選択[6] にＢをセットする
次はナップザック (6) に B を詰めます。値を計算すると 8 になり、今度は金額[6] の値 6 より大きくなります。つまり、A を詰め込むよりも B を詰め込む方が金額が高くなるのです。金額[6] と選択[6] の値を更新します。ナップザック (7) の場合も同様ですね。
あとは、順番に同じことを繰り返して、配列の値を更新していきます。そして、品物 C を最後まで詰め込むと、次のようになります。
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│Ｃ│Ｃ│Ｂ│Ａ│Ｃ│Ｂ│Ｂ│Ａ│Ｃ│Ｂ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│１│２│４│６│７│８│10│12│13│14│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
このときの金額[10] の 14 が答となります。この状態からナップザックに詰め込まれた品物を求めます。
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│Ｃ│Ｃ│Ｂ│Ａ│Ｃ│Ｂ│Ｂ│Ａ│Ｃ│Ｂ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴┼┘
↑ │ │ │
└───────┴←────┴←────┘
まず、選択[10] にセットされた品物を取り出します。この場合は B ですね。次は、10 からＢの大きさ 3 を引いた 7 のときに追加された品物を取り出します。この場合も B ですね。同様に、7 から 3 を引いた 4 のときに追加された品物を求めます。これはＡですね。4 から A の大きさを引くと 0 になるので、これ以上品物は入っていません。したがって、ナップザックの中には A が 1 個、B が 2 個入っていることがわかります。
それでは、実際にプログラミングしましょう。まず、品物を表すデータ型 Item を定義します。
リスト : 品物の定義
data Item = None | Item {name :: String, price :: Integer, size :: Integer} deriving (Eq, Show)
item_list = [Item "A" 6 4, Item "B" 4 3, Item "C" 1 1]
name が品名、price が金額で、size が大きさを表します。None は品物を選んでいないことを表すために使います。item_list は品物 (Item) を格納したリストです。
次は、ナップザック問題の解を求める関数 knapsack を作ります。
リスト : ナップザック問題の解法
knapsack :: Integer -> IO ()
knapsack ksize = do
gain <- newArray (0, ksize) 0 :: IO (IOArray Integer Integer)
choice <- newArray (0, ksize) None :: IO (IOArray Integer Item)
mapM_ (\item -> mapM_ (\j -> do a <- readArray gain (j - size item)
b <- readArray gain j
when (b < a + price item) $ do
writeArray gain j (a + price item)
writeArray choice j item)
[size item .. ksize])
item_list
-- 結果の表示
print_answer ksize choice
where print_answer :: Integer -> IOArray Integer Item -> IO ()
print_answer i choice = do
item <- readArray choice i
when (item /= None) $ do
print item
print_answer (i - size item) choice
knapsack の引数 ksize はナップザックの大きさを表します。最初に、金額を表す配列 gain と選択した品物を格納する配列 choice を用意します。gain は 0 で、choice は None で初期化します。次の mapM_ で、品物を一つずつ item_list から取り出します。次の mapM_ で、item の size から ksize まで gain と choice を更新していきます。
まず、ナップザックから item size 分だけ空けた場合の金額を a にセットします。a に item を追加した値 (a + price item) が新しい値になります。この値が現在の値 b よりも大きければ、choice と gain の値を writeArray で更新します。
最後にナップザックの中身を表示します。まず choice の i 番目に格納された品物を item にセットします。それが None でなければ、print で Item を表示し、i から size item を引いて、次に詰め込まれた品物へ移ります。
それでは実行結果を示します。
ghci> knapsack 10
Item {name = "B", price = 4, size = 3}
Item {name = "B", price = 4, size = 3}
Item {name = "A", price = 6, size = 4}