M.Hiroi's Home Page

Functional Programming

お気楽 Haskell プログラミング入門

[ PrevPage | Haskell | NextPage ]

動的計画法

今回は Haskell で「動的計画法 (dynamic programming)」というアルゴリズムに挑戦してみましょう。難しそうな名前がついていますが、これに惑わされてはいけません。動的計画法は、大きな問題を小さな問題に分けて、それを一つ一つ解いていくことで大きな問題を解く方法です。

問題によっては、小さな問題に分割していくと同じ小問題が何回も現れる場合があります。この場合、同じ問題を何回も解くよりも、その解を表などに保存しておいて、必要なときにその表から答を求めた方が、効率良く問題を解くことができるはずです。

どうせ小問題を解かなければならないのであれば、はじめから必要になりそうな小問題を解いて表を埋めておいたほうが、プログラムを作りやすい場合もあります。このように、与えられた問題を解くために小問題の表を埋めてしまう、というのが「動的計画法」の基本的な考え方です。

●組み合わせの数

簡単な例題として、組み合わせの数を求めるプログラムを作ってみましょう。組み合わせの数 nr を求めるには、次の公式を使えば簡単です。

(1)
nr = n * (n - 1) * (n - 2) * ... * (n - r + 1) / (1 * 2 * 3 * ... * r)
(2)
n0 = nn = 1
nr = nr-1 * (n - r + 1) / r
(3)
n0 = nn = 1
nr = n-1r-1 + n-1r

皆さんお馴染みの公式ですね。(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)
*Main> :set +s
*Main> comb 20 10
184756
(0.40 secs, 133,097,264 bytes)
*Main> comb 22 11
705432
(1.70 secs, 507,979,760 bytes)
*Main> comb 24 12
2704156
(5.80 secs, 1,947,061,776 bytes)

実行環境 : Ubunts 18.04 (WSL), Intel Core i5-6200U 2.30GHz

このように、公式 (3) を使うと時間がかかります。公式 (1), (2) を使えば高速に答えを求めることができますが、今回は動的計画法の例題として、あえてこのプログラムの高速化に挑戦してみましょう。

●動的計画法による高速化

公式からわかるように、nr の値は n-1rn-1r-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) の値を求めて返します。

それでは実際に試してみましょう。

*Main> combDP 100 50
100891344545564193334812497256
(0.02 secs, 14,029,456 bytes)
*Main> combDP 200 100
90548514656103281165404177077484163874504589675413336841320
(0.07 secs, 56,249,320 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 回呼び出すだけです。

それでは実行してみましょう。

*Main> combDP'' 100 50
100891344545564193334812497256
(0.00 secs, 1,275,736 bytes)
*Main> combDP'' 200 100
90548514656103281165404177077484163874504589675413336841320
(0.01 secs, 4,932,544 bytes)

●整数の分割

もうひとつ、数値計算の例を取り上げましょう。整数 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(n, k) = 0                          ; n < 0 または k < 1
p(n, k) = 1                          ; n = 0 または k = 1
p(n, k) = p(n - k, k) + p(n, k - 1)

たとえば、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)
*Main> partition_number 40
37338
(0.14 secs, 45,661,744 bytes)
*Main> partition_number 50
204226
(0.82 secs, 241,760,864 bytes)
*Main> partition_number 60
966467
(3.77 secs, 1,121,763,136 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

説明をそのままプログラムしただけなので、とくに難しいところはないと思います。それでは実際に試してみましょう。

*Main> partition_number' 40
37338
(0.00 secs, 1,845,976 bytes)
*Main> partition_number' 50
204226
(0.00 secs, 2,854,832 bytes)
*Main> partition_number' 60
966467
(0.01 secs, 4,088,552 bytes)
*Main> partition_number' 100
190569292
(0.02 secs, 11,281,640 bytes)
*Main> partition_number' 1000
24061467864032622473692149727991
(1.24 secs, 1,150,099,952 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 の要素を先頭から順番に取り出していくことができます。

それでは実行してみましょう。

*Main> partition_number'' 40
37338
(0.00 secs, 268,800 bytes)
*Main> partition_number'' 50
204226
(0.00 secs, 375,656 bytes)
*Main> partition_number'' 60
966467
(0.00 secs, 504,176 bytes)
*Main> partition_number'' 100
190569292
(0.00 secs, 1,244,464 bytes)
*Main> partition_number'' 1000
24061467864032622473692149727991
(0.19 secs, 165,321,496 bytes)

●ナップザック問題

最後に、簡単な例題として「ナップザック問題」を取り上げます。ナップザック (knapsack) とは辞書を引いてみると、ランドセルのような背中にせおう四角形の袋や箱のことを意味します。ここでは物を入れる袋と簡単に考えてください。

ここで、ナップザックの中に品物を詰め込むことを考えてみます。一つのナップザックと複数の品物が与えられたとき、袋に詰めた品物の合計金額が最大になるような選び方を求めることが「ナップザック問題」です。ここでは、同じ品物を何個も選んでもいいのですが、ナップザックの大きさをオーバーしてはいけません。

実はこの「ナップザック問題」が「NP 問題」なのです。世の中にはさまざまな問題が山積していますが、スーパーコンピュータを使っても解くのに数億年かかる、というような難問が「NP 問題」です。これは、厳密に解を求めようとすると、全ての場合について総当たりで調べるしか方法がなく、データ数が多くなると時間がかかるため、現実的な時間では解答を出すことができないというものです。品物の詰め方が難問の一つ、といわれてもピンとこないと思いますが、ナップザック問題は品物の種類が増えるにしたがって、その組み合わせ方が爆発的に増えるのです。

ところが、幸いなことに「ナップザック問題」は実用的には解決済みの問題と考えられています。とくに有名なのが「動的計画法」を用いた解法です。ナップザックと品物の大きさを整数値に限定すれば、動的計画法を用いることで厳密解を求めることができるのです。

それでは具体的に、ナップザック問題に動的計画法を適用してみましょう。ナップザックの大きさは 10 で、次の 3 種類の品物を詰め込むことにします。

(A)大きさ  4  :金額  6
(B)大きさ  3  :金額  4
(C)大きさ  1  :金額  1

まず、大きさが 0 から 10 までのナップザックを用意します。これらのナップザックに品物を順番に詰め込んで、その合計金額を配列に格納しておきます。この配列は品物を詰め込んでいない状態 (金額は全て 0) に初期化します。

最初に品物 A を詰め込みます。このとき、小さなナップザックから順番に詰め込んでいきます。

品物 A が入る大きさのナップザックから詰め込みます。品物を入れるときは、それより小さいナップザックには、その時点での最適な値が決定されていると考えます。ナップザック (4) に品物 A を詰め込む場合、A の大きさだけ空いているナップザック (0) の状態に詰め込めば、ナップザック (4) の最適な値を求めることができるはずです。このように、前に計算された値を使うところが動的計画法の特徴なのです。

具体的には、金額[0] に A の金額 6 を足した値を計算し、金額[4] より大きくなれば金額[4] をその値に更新します。もし、金額[4] より小さいのであれば金額[4] は更新しません。つまり、品物Aはナップザック (4) には詰め込まないのです。ほかの組み合わせの方が正解だというわけです。

詰め込んだ品物を記憶しておくため、もう一つ配列を用意して、そこに追加した品物の種類を格納しておきます。ナップザックの中身全てを記憶しておく必要はありません。この配列を使って、あとからナップザックの中身を求めることができます。

次に、品物 B を詰め込んでいきます。

まず、ナップザック (3) に B が詰め込まれます。これは品物 A の場合と同じですね。次に、ナップザック (4) に B を詰めようとします。その値を計算すると 4 となり、金額[4] の値 6 より小さいので、B は詰め込みません。ナップザック (5) の場合も同様です。

次はナップザック (6) に B を詰めます。値を計算すると 8 になり、今度は金額[6] の値 6 より大きくなります。つまり、A を詰め込むよりも B を詰め込む方が金額が高くなるのです。金額[6] と選択[6] の値を更新します。ナップザック (7) の場合も同様ですね。

あとは、順番に同じことを繰り返して、配列の値を更新していきます。そして、品物 C を最後まで詰め込むと、次のようになります。

このときの金額[10] の 14 が答となります。この状態からナップザックに詰め込まれた品物を求めます。

まず、選択[10] にセットされた品物を取り出します。この場合は B ですね。次は、10 からBの大きさ 3 を引いた 7 のときに追加された品物を取り出します。この場合も B ですね。同様に、7 から 3 を引いた 4 のときに追加された品物を求めます。これはAですね。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 を引いて、次に詰め込まれた品物へ移ります。

●実行結果

それでは実行結果を示します。

*Main> knapsack 10
Item {name = "B", price = 4, size = 3}
Item {name = "B", price = 4, size = 3}
Item {name = "A", price = 6, size = 4}

初版 2013 年 4 月 28 日
改訂 2021 年 8 月 8 日

部分和問題

部分和問題は、要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題です。たとえば、集合 {2, 3, 5, 8} の場合、総和が 10 となる部分集合は {2, 3, 5} と {2, 8} がありますが、14 となる部分集合はありません。部分集合の総数は、要素数を n とすると 2n 個になるので、n が大きくなるとナイーブな方法では時間がかかってしまいます。実際には、「分岐限定法」や「動的計画法」を使うことで、現実的な時間で部分和問題を解くことができます。

●ナイーブな方法

最初にナイーブな方法で部分和問題を解いてみましょう。今回は要素を正整数に限定します。部分和問題は「べき集合」を生成する関数 powerSet を作ると簡単です。次のリストを見てください。

リスト : べき集合

powerSet :: [a] -> [[a]]
powerSet []     = [[]]
powerSet (x:xs) =  powerSet xs ++ [x:ys | ys <- powerSet xs]

べき集合を求める関数 powerSet は簡単です。引数が空リストの場合は [ ] を格納したリストを返します。そうでなければ、引数を x:xs で分解します。 そして、powerSet を再帰呼び出しして xs のべき集合を求め、その集合に先頭要素 x を追加します。そして、その集合と xs のべき集合を演算子 ++ で連結します。

簡単な実行例を示します。

*Main> powerSet [2,3,5,8]
[[],[8],[5],[5,8],[3],[3,8],[3,5],[3,5,8],[2],[2,8],[2,5],[2,5,8],[2,3],[2,3,8],[2,3,5],[2,3,5,8]]

[2, 3, 5, 8] の部分集合は空集合 [ ] を含めて 16 通りあります。この powerSet を使うと部分和問題のプログラムは次のようになります。

リスト : 部分和問題

subsetSum :: Integer -> [Integer] -> [[Integer]]
subsetSum n xs = filter ((==n) . sum) $ powerSet xs

部分集合の総和を sum で求め、n と等しいものを filter で取り出します。それでは実行してみましょう。

*Main> subsetSum 10 [2,3,5,8]
[[2,8],[2,3,5]]
*Main> subsetSum 14 [2,3,5,8]
[]

とても簡単ですね。ただし、集合の要素数が多くなると、実行時間がかかるようになります。次のテストプログラムを見てください。

リスト : 簡単なテスト

nums :: [Integer]
nums = [  1,   2,   3,   5,   8,   13,   21,   34,   55,    89,
        144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]

test :: Int -> [[Integer]]
test n = subsetSum m xs
  where xs = take n nums
        m  = sum xs - 1

配列 nums はフィボナッチ数列になっています。要素の総和を M とすると、1 から M までの整数は、要素を組み合わせて必ず作ることができます。これはフィボナッチ数列の面白い特徴です。

テストは 総和 - 1 となる組み合わせを subsetSum で求め、その実行時間を計測します。結果は次のようになりました。

*Main> :set +s
*Main> test 16
[[2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597]]
(0.35 secs, 153,202,928 bytes)
*Main> test 17
[[2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584]]
(0.62 secs, 323,071,040 bytes)
*Main> test 18
[[2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]]
(1.26 secs, 679,590,728 bytes)
*Main> test 19
[[2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765]]
(2.63 secs, 1,426,180,688 bytes)
*Main> test 20
[[2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765,10946]]
(5.56 secs, 2,986,466,360 bytes)

実行環境 : Ubunts 18.04 (WSL), Intel Core i5-6200U 2.30GHz

要素がひとつ増えると実行時間は約 2 倍になっていることがわかります。要素数を n とすると、subsetSum の実行時間は 2n に比例する遅いプログラムなのです。

●分岐限定法による高速化

次は「分岐限定法」を使って部分和問題の高速化に挑戦してみましょう。分岐限定法を簡単に説明すると、深さ優先探索において不要な局面の生成を省くため枝刈りを行う方法です。思考ルーチンでよく使われる「アルファベータ法」や反復深化の高速化で用いられる「下限値枝刈り法」も分岐限定法の一種です。

今回の部分和問題は要素を正整数値に限定しているので、二種類の枝刈りを考えることができます。ひとつは部分集合の総和が求める値 N を超えた場合です。残りの要素は正整数なので、これ以上要素を追加しても解を得られないのは明白ですね。もうひとつは、部分集合の総和に残りの要素をすべて足しても N に満たない場合です。これも解を得られないのは明白です。

部分集合の総和を S, 残りの要素の総和を R, 求める値を N とすると、探索を行う条件は次の式で表すことができます。

S < N && S + R >= N

これをプログラムすると次のようになります。

リスト : 部分和問題 (分岐限定法)

subsetSum1 :: Integer -> [Integer] -> [[Integer]]
subsetSum1 n xs = iter [] 0 xs (sum xs) []
  where iter a m [] _ b
          | n == m    = a : b
          | otherwise = b
        iter a m (y:ys) r b
          | n == m = a : b
          | m < n && m + r >= n =
              iter a m ys (r - y) $ iter (y:a) (m + y) ys (r - y) b
          | otherwise = b

実際の処理は局所関数 iter で行います。引数 m は求めている部分集合の総和、r は残りの要素の総和を表します。引数 b は結果を格納する累積変数です。m < n かつ m + r >= n ならば条件を満たすので、iter を再帰呼び出しして探索を続行します。とても簡単ですね。

それでは試してみましょう。

リスト : テストプログラム

test1 :: Int -> [[Integer]]
test1 n = subsetSum1 m xs
  where xs = take n nums
        m  = sum xs - 1
*Main> test1 16
[[1597,987,610,377,233,144,89,55,34,21,13,8,5,3,2]]
(0.01 secs, 142,072 bytes)
*Main> test1 17
[[2584,1597,987,610,377,233,144,89,55,34,21,13,8,5,3,2]]
(0.00 secs, 148,304 bytes)
*Main> test1 18
[[4181,2584,1597,987,610,377,233,144,89,55,34,21,13,8,5,3,2]]
(0.00 secs, 154,536 bytes)
*Main> test1 19
[[6765,4181,2584,1597,987,610,377,233,144,89,55,34,21,13,8,5,3,2]]
(0.01 secs, 160,768 bytes)
*Main> test1 20
[[10946,6765,4181,2584,1597,987,610,377,233,144,89,55,34,21,13,8,5,3,2]]
(0.00 secs, 167,736 bytes)

とても速くなりましたね。ただし、これは枝刈りがうまくいった場合であり、データによっては枝刈りが機能しない場合もありえます。たとえば、xs の最後尾の要素から 1 を引いた値を判定してみましょう。結果は次のようになりました。

リスト : テストプログラム

test1' :: Int -> [[Integer]]
test1' n = subsetSum1 m xs
  where xs = take n nums
        m  = last xs - 1
*Main> test1' 16
[[987,377,144,55,21,8,3,1]]
(0.15 secs, 51,641,120 bytes)
*Main> test1' 17
[[1597,610,233,89,34,13,5,2]]
(0.29 secs, 103,197,408 bytes)
*Main> test1' 18
[[2584,987,377,144,55,21,8,3,1]]
(0.61 secs, 206,310,000 bytes)
*Main> test1' 19
[[4181,1597,610,233,89,34,13,5,2]]
(1.40 secs, 412,531,248 bytes)
*Main> test1' 20
[[6765,2584,987,377,144,55,21,8,3,1]]
(2.23 secs, 824,973,760 bytes)

実行時間はそれほど速くなりません。実をいうと、データの並び方によって実行時間は大きく左右されます。興味のある方はいろいろ試してみてください。

●動的計画法による高速化

今度は前回説明した「動的計画法」で部分和問題を解いてみましょう。部分和問題の場合、総和が N となる部分集合があるか判定するだけでよければ、動的計画法で解くことが可能です。

部分和問題の場合、要素をひとつずつ追加しながら、総和 N となる部分集合があるか判定します。簡単な例を示しましょう。次の図を見てください。

上図は xs = {2, 3, 5, 8} で N = 10 の部分集合があるか判定する場合です。最初に N + 1 の配列 Ai を用意します。空集合の総和は 0 なので A0[0] に○をセットします。次に、要素 2 を追加します。部分集合は { } と {2} になります。A1[0] と A1[2] に○をセットします。その次に要素 3 を追加します。追加される部分集合は {3} と {2, 3} になるので、A2[0], A2[2], A2[3] と A2[5] に○をセットします。

つまり、i 番目の要素 x を追加する場合、Ai-1 で○が付いている位置を y とすると、Ai[y] と Ai[x + y] に○をセットすればいいわけです。添字 y は部分集合の総和を表しています。Ai[y] に○をセットすることは、その部分集合に x を加えないことを意味し、Ai[x + y] に○をセットすることは、その部分集合に x を追加することを意味するわけです。

次に 5 を追加します。A2 の○の位置は 0, 2, 3, 5 なので、これに 5 を足した 5, 7, 8, 10 の位置に○をセットします。最後に 8 を追加します。A3 の○の位置は 0, 2, 3, 5, 7, 8, 10 なので、これに 8 を足した 8, 10 の位置に○をセットします。A4[10] の値が○になので、部分和が 10 となる部分集合があることがわかります。

もうひとつ簡単な例を示しましょう。今度は総和が 14 となる部分集合があるか判定します。

3 番目で○の位置は 0, 2, 3, 5, 7, 8, 10 です。次は 8 を追加しますが、総和 14 より大きい値は不要なので、8, 10, 11, 13 の位置に○を追加します。14 の位置は×なので、総和が 14 となる部分集合は無いことがわかります。

プログラムは次のようになります。

リスト : 部分和問題 (動的計画法)

subsetSum2 :: Integer -> [Integer] -> IO Bool
subsetSum2 n xs = do
  a <- newArray (0, n) False :: IO (IOArray Integer Bool)
  writeArray a 0 True
  mapM_ (\x -> mapM_ (\y -> do f <- readArray a y
                               when f $ do writeArray a (x + y) True)
               [n - x, n - x - 1 .. 0])
        xs
  readArray a n

○を True で、×を False で表します。配列をひとつで済ますため、配列の後ろから True の位置を検索していることに注意してください。また、検索の開始位置を n - x とすることで、True をセットするときの範囲チェックを省略しています。今回のプログラムでは xs の要素をすべてチェックしていますが、x + y が n と等しくなったら True を返えすようにプログラムしてもかまいません。あとは特に難しいところはないと思います。

それでは実際に試してみましょう。テストプログラムを示します。

リスト : テストプログラム

test2 :: Int -> IO Bool
test2 n = subsetSum2 m xs
  where xs = take n nums
        m  = sum xs - 1

実行結果は次のようになりました。

*Main> test2 16
True
(0.10 secs, 69,407,536 bytes)
*Main> test2 17
True
(0.16 secs, 119,313,936 bytes)
*Main> test2 18
True
(0.21 secs, 204,414,184 bytes)
*Main> test2 19
True
(0.51 secs, 349,146,136 bytes)
*Main> test2 20
True
(0.95 secs, 594,712,440 bytes)

ナイーブな方法よりも高速になりましたが、分岐限定法にはかないませんでした。集合の要素数を M, 総和を N とすると、今回のプログラムの実行速度は N * M に比例します。たとえば、N の値を 最後尾の要素 - 1 とすると、実行結果は次のようになります。

リスト : テストプログラム

test2' :: Int -> IO Bool
test2' n = subsetSum2 m xs
  where xs = take n nums
        m  = last xs - 1
*Main> test2' 16
True
(0.04 secs, 24,251,792 bytes)
*Main> test2' 17
True
(0.05 secs, 41,883,560 bytes)
*Main> test2' 18
True
(0.08 secs, 72,073,256 bytes)
*Main> test2' 19
True
(0.15 secs, 123,608,560 bytes)
*Main> test2' 20
True
(0.27 secs, 211,342,688 bytes)

N が小さくなったので、実行時間も速くなりました。このように、動的計画法では N が大きくなると、どうしても時間がかかるようになります。そこで、配列を使わずにプログラムを作ってみましょう。

●Data.Set を使う方法

配列を使う場合、N が大きくなると True を検索する処理に時間がかかるようになります。そこで、配列の代わりにモジュール Data.Set を使って部分集合の総和を管理することにします。プログラムは次のようになります。

リスト : 部分和問題

import qualified Data.Set as S

subsetSum3 :: Integer -> [Integer] -> Bool
subsetSum3 n xs = iter xs (S.singleton 0)
  where iter [] s = S.member n s
        iter (x:xs) s = iter xs $ foldl (\s' y -> if x + y <= n
                                                  then S.insert (x + y) s'
                                                  else s')
                                        s
                                        (S.elems s)

実際の処理は局所関数 iter で行います。最初に部分集合の総和を管理するセット s を 0 に初期化します。関数 S.elems で集合の要素を求め、それを foldl に渡します。ラムダ式の中で、リストの要素 x と集合の要素 y を足し算し、それが n 以下であれば S.insert で集合に追加します。最後に、S.member で集合 s に n が含まれているかチェックします。

それでは実行結果を示します。

リスト : テストプログラム

test3 :: Int -> Bool
test3 n = subsetSum3 m xs
  where xs = take n nums
        m  = sum xs - 1

test3' :: Int -> Bool
test3' n = subsetSum3 m xs
  where xs = take n nums
        m  = last xs - 1
*Main> test3 16
True
(0.03 secs, 6,501,392 bytes)
*Main> test3 17
True
(0.04 secs, 10,849,880 bytes)
*Main> test3 18
True
(0.03 secs, 18,217,672 bytes)
*Main> test3 19
True
(0.06 secs, 30,687,048 bytes)
*Main> test3 20
True
(0.08 secs, 51,411,680 bytes)
*Main> test3' 16
True
(0.01 secs, 3,353,728 bytes)
*Main> test3' 17
True
(0.02 secs, 5,570,240 bytes)
*Main> test3' 18
True
(0.01 secs, 9,308,936 bytes)
*Main> test3' 19
True
(0.04 secs, 15,544,152 bytes)
*Main> test3' 20
True
(0.04 secs, 25,952,104 bytes)

subsetSum2 よりも高速になりました。ちなみに、分岐限定法と同様の枝刈りを入れると、実行速度はもう少し速くなります。興味のある方は試してみてください。


初版 2013 年 4 月 28 日
改訂 2021 年 8 月 8 日

Copyright (C) 2013-2021 Makoto Hiroi
All rights reserved.

[ PrevPage | Haskell | NextPage ]