今回は「動的計画法 (dynamic programming) 」というアルゴリズムについて取り上げます。難しそうな名前がついていますが、これに惑わされてはいけません。動的計画法は、大きな問題を小さな問題に分けて、それを一つ一つ解いていくことで大きな問題を解く方法です。
問題によっては、小さな問題に分割していくと同じ小問題が何回も現れる場合があります。この場合、同じ問題を何回も解くよりも、その解を表などに保存しておいて、必要なときにその表から答を求めた方が、効率良く問題を解くことができるはずです。
どうせ小問題を解かなければならないのであれば、はじめから必要になりそうな小問題を解いて表を埋めておいたほうが、プログラムを作りやすい場合もあります。このように、与えられた問題を解くために小問題の表を埋めてしまう、というのが「動的計画法」の基本的な考え方です。
簡単な例題として、組み合わせの数を求めるプログラムを作ってみましょう。組み合わせの数 nＣr を求めるには、次の公式を使えば簡単です。
(1) nＣr = n * (n - 1) * (n - 2) * ... * (n - r + 1) / (1 * 2 * 3 * ... * r)
(2) nＣ0 = nＣn = 1 nＣr = nＣr-1 * (n - r + 1) / r
(3) nＣ0 = nＣn = 1 nＣr = n-1Ｃr-1 + n-1Ｃr
皆さんお馴染みの公式ですね。(1) と (2) の公式を使うと簡単に (高速に) 答えを求めることができます。ただし、(3) の公式をそのままプログラムすると二重再帰になるので、大きな値を求めると時間がかかってしまいます。実際にプログラムを作って確かめてみましょう。
リスト : 組み合わせの数
(defun comb (n r)
(if (or (zerop r) (= n r))
1
(+ (comb (1- n) r) (comb (1- n) (1- r)))))
* (time (comb 26 13)) Evaluation took: 0.828 seconds of real time 0.828125 seconds of total run time (0.828125 user, 0.000000 system) 100.00% CPU 1,155,687,460 processor cycles 0 bytes consed 10400600 * (time (comb 28 14)) Evaluation took: 3.172 seconds of real time 3.171875 seconds of total run time (3.171875 user, 0.000000 system) 100.00% CPU 4,450,983,445 processor cycles 0 bytes consed 40116600 * (time (comb 30 15)) Evaluation took: 12.313 seconds of real time 12.312500 seconds of total run time (12.312500 user, 0.000000 system) 99.99% CPU 17,229,207,758 processor cycles 0 bytes consed 155117520 (Windows XP, celeron 1.40 GHz, SBCL ver 1.0.37)
このように SBCL (ver 1.0.37) でも時間がかかります。公式 (1), (2) を使えば高速に答えを求めることができますが、今回は動的計画法の例題として、あえてこのプログラムの高速化に挑戦してみましょう。
公式からわかるように、nＣr の値は n-1Ｃr と n-1Ｃr-1 を足したものです。n = 0 から順番に組み合わせの数を求めて表に格納しておけば、n が大きな値でも簡単に求めることができるはずです。プログラムは次のようになります。
リスト : 組み合わせの数 (動的計画法[1])
(defvar *comb-table* nil)
(defun comb-dp-sub (n r)
(if (or (zerop r) (= n r))
1
(+ (aref *comb-table* (1- n) r) (aref *comb-table* (1- n) (1- r)))))
(defun comb-dp (n r)
(setq *comb-table* (make-array (list (1+ n) (1+ n))))
(do ((i 0 (1+ i)))
((< n i) (aref *comb-table* n r))
(do ((j 0 (1+ j)))
((< i j))
(setf (aref *comb-table* i j) (comb-dp-sub i j)))))
大域変数 *comb-table* に値を格納する 2 次元配列をセットします。関数 comb-dp-sub は *comb-table* から組み合わせの数を求めます。組み合わせの数を (n, r) で表すことにすると、関数 comb-dp は (0, 0) から順番に、(1, 0), (1, 1), (2, 0), (2, 1), (2, 2) ... と組み合わせの数を求めて *comb-table* にセットします。ようするに、「パスカルの三角形」を作っていくわけです。最後に *comb-table* から (n, r) の値を求めて返します。
それでは実際に試してみましょう。
* (time (comb-dp 30 15)) Evaluation took: 0.000 seconds of real time 0.000000 seconds of total run time (0.000000 user, 0.000000 system) 100.00% CPU 119,330 processor cycles 4,096 bytes consed 155117520 * (time (comb-dp 100 50)) Evaluation took: 0.000 seconds of real time 0.000000 seconds of total run time (0.000000 user, 0.000000 system) 100.00% CPU 2,413,489 processor cycles 110,328 bytes consed 100891344545564193334812497256 (Windows XP, celeron 1.40 GHz, SBCL ver 1.0.37)
動的計画法の効果はとても高いですね。なお、表は二次元配列ではなくベクタで済ますこともできます。たとえば、(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])
(defun comb-dp1 (n r)
(let ((table (make-array (1+ n) :initial-element 1)))
(do ((i 1 (1+ i)))
((< n i) (aref table r))
(do ((j (1- i) (1- j)))
((zerop j))
(incf (aref table j) (aref table (1- j)))))))
ベクタの値を書き換えていくので、ベクタの後方から計算していくことに注意してください。前方から計算すると、値がおかしくなります。
なお、関数 comb はメモ化関数を使って高速化することもできます。次のリストを見てください。
リスト : 組み合わせの数 (メモ化)
; メモ化関数
(defun memoize (func)
(let ((table (make-hash-table :test #'equal)))
#'(lambda (&rest args)
(let ((value (gethash args table nil)))
(unless value
(setf value (apply func args))
(setf (gethash args table) value))
value))))
; メモ化
(setf (symbol-function 'comb) (memoize #'comb))
プログラムの説明は拙作のページ メモ化と遅延評価 をお読みください。
動的計画法とメモ化関数は、どちらも表を使って計算を高速化する方法です。動的計画法は小さな問題の解を積み上げていく、つまりボトムアップな方法なのに対し、メモ化関数は大きな問題を小さな問題に分割していく、トップダウンな方法 (分割統治法) といえます。
もうひとつ、数値計算の例を取り上げましょう。整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。簡単な例を示しましょう。次の図を見てください。
─┬─ 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 の分割
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) を使うと次のようにプログラムすることができます。
リスト : 分割数
(defun part-num (n k)
(cond ((or (zerop n) (= n 1) (= k 1)) 1)
((or (< n 0) (< k 1)) 0)
(t (+ (part-num (- n k) k) (part-num n (- k 1))))))
(defun partition-number (n)
(part-num n n))
関数 part-num は p(n, k) の定義をそのままプログラムしただけです。ただし、このプログラムは二重再帰で何度も同じ値を求めているため実行速度はとても遅くなります。次のように関数 part-num をメモ化することで高速化することができます。
リスト : 分割数 (メモ化による高速化) ; メモ化 (setf (symbol-function 'part-num) (memoize #'part-num)) (defun partition-number1 (n) (part-num n n))
動的計画法を使うと、もっと速くなります。次の図を見てください。
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) の値を求めることができます。
プログラムは次のようになります。
リスト : 分割数 (動的計画法)
(defun partition-number2 (n)
(let ((a (make-array (+ n 1) :initial-element 1)))
(do ((k 2 (1+ k)))
((< n k) (aref a n))
(do ((m k (1+ m)))
((< n m))
(incf (aref a m) (aref a (- m k)))))))
説明をそのままプログラムしただけなので、とくに難しいところはないと思います。実行例を示します。
* (time (partition-number1 1000)) Evaluation took: 6.875 seconds of real time 6.875000 seconds of total run time (6.578125 user, 0.296875 system) [ Run times consist of 0.187 seconds GC time, and 6.688 seconds non-GC time. ] 100.00% CPU 9,623,093,846 processor cycles 67,090,904 bytes consed 24061467864032622473692149727991 * (time (partition-number2 1000)) Evaluation took: 0.219 seconds of real time 0.218750 seconds of total run time (0.140625 user, 0.078125 system) [ Run times consist of 0.125 seconds GC time, and 0.094 seconds non-GC time. ] 100.00% CPU 307,823,617 processor cycles 11,854,264 bytes consed 24061467864032622473692149727991 (Windows XP, celeron 1.40 GHz, SBCL ver 1.0.37)
部分和問題は、要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題です。たとえば、集合 {2, 3, 5, 8} の場合、総和が 10 となる部分集合は {2, 3, 5} と {2, 8} がありますが、14 となる部分集合はありません。部分集合の総数は、要素数を n とすると 2n 個になるので、n が大きくなるとナイーブな方法では時間がかかってしまいます。実際には、「分岐限定法」や「動的計画法」を使うことで、現実的な時間で部分和問題を解くことができます。
最初にナイーブな方法で部分和問題を解いてみましょう。今回は要素を正整数に限定します。部分和問題は「べき集合」を生成する高階関数 power-set を用意すると簡単です。次のリストを見てください。
リスト : べき集合 (defun power-set (func xs &optional (a nil)) (cond ((null xs) (funcall func (reverse a))) (t (power-set func (cdr xs) a) (power-set func (cdr xs) (cons (car xs) a)))))
引数 func が関数で、xs が集合を表すリスト、a が累積変数 (リスト) です。処理内容は簡単で、xs が空リストの場合、部分集合がひとつできたので func に (reverse a) を渡して評価します。あとは、xs の先頭要素を選ばない場合は power-set を再帰呼び出しし、選ぶ場合は a に (car xs) を追加してから power-set を再帰呼び出しします。これですべての部分集合を生成することができます。
簡単な実行例を示します。
* (power-set #'print '(1 2 3 4)) NIL (4) (3) (3 4) (2) (2 4) (2 3) (2 3 4) (1) (1 4) (1 3) (1 3 4) (1 2) (1 2 4) (1 2 3) (1 2 3 4) (1 2 3 4) *
最後の (1 2 3 4) は power-set の返り値です。部分集合は空集合 NIL を含めて 16 通りあります。この power-set を使うと部分和問題のプログラムは次のようになります。
リスト : 部分和問題 (defun subset-sum0 (xs n) (power-set #'(lambda (ys) (when (= (apply #'+ ys) n) (print ys))) xs))
部分集合 ys の総和を (apply #'+ ys) で求め、n と等しい場合は print で出力します。それでは実行してみましょう。
* (subset-sum0 '(2 3 5 8) 10) (2 8) (2 3 5) NIL * (subset-sum0 '(2 3 5 8) 14) NIL *
とても簡単ですね。ただし、集合の要素数が多くなると、実行時間がかかるようになります。次のテストプログラムを見てください。
リスト : 簡単なテスト
(defvar *fibo* '(1 2 3 5 8 13 21 34 55 89 144 233 377 610 987
1597 2584 4181 6765 10946 17711 28657 46368
75025 121393 196418 317811 514229 832040 1346269))
;
(defun test (func)
(dolist (x '(20 21 22 23 24))
(let ((xs (subseq *fibo* 0 x))
(s (get-internal-real-time)))
(funcall func xs (1- (apply #'+ xs)))
(print (/ (- (get-internal-real-time) s) 1000.0)))))
リスト *fibo* はフィボナッチ数列になっています。要素の総和を M とすると、1 から M までの整数は、要素を組み合わせて必ず作ることができます。これはフィボナッチ数列の面白い特徴です。
テストは数列の長さを 20 から一つずつ増やしながら、総和 - 1 となる組み合わせを subset-sum0 で求め、その実行時間を計測します。結果は次のようになりました。
* (test #'subset-sum0) (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946) 0.859 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711) 1.719 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657) 3.656 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368) 7.547 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025) 15.625 NIL (Windows XP, celeron 1.40 GHz, SBCL ver 1.0.37, 単位 : 秒)
要素がひとつ増えると実行時間は約 2 倍になっていることがわかります。要素数を n とすると、subset-sum0 の実行時間は 2n に比例する遅いプログラムなのです。
次は「分岐限定法」を使って部分和問題の高速化に挑戦してみましょう。分岐限定法を簡単に説明すると、深さ優先探索において不要な局面の生成を省くため枝刈りを行う方法です。思考ルーチンでよく使われる「アルファベータ法」や反復深化の高速化で用いられる「下限値枝刈り法」も分岐限定法の一種です。
今回の部分和問題は要素を正整数値に限定しているので、二種類の枝刈りを考えることができます。ひとつは部分集合の総和が求める値 N を超えた場合です。残りの要素は正整数なので、これ以上要素を追加しても解を得られないのは明白ですね。もうひとつは、部分集合の総和に残りの要素をすべて足しても N に満たない場合です。これも解を得られないのは明白です。
これをプログラムすると次のようになります。
リスト : 部分和問題 (分岐限定法)
(defun subset-sum1 (xs n)
(labels ((iter (xs s r a)
(cond ((= s n)
(print (reverse a)))
((<= s n (+ s r))
(let ((x (car xs)))
(iter (cdr xs) s (- r x) a)
(iter (cdr xs) (+ s x) (- r x) (cons x a)))))))
(iter xs 0 (apply #'+ xs) nil)))
実際の処理は局所関数 iter で行います。引数 s は求めている部分集合の総和、r は残りの要素の総和を表します。s < n かつ s + r >= n ならば条件を満たすので、iter を再帰呼び出しして探索を続行します。とても簡単ですね。
それでは試してみましょう。
* (test #'subset-sum1) (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946) 0.0 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711) 0.0 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657) 0.0 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368) 0.0 (2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025) 0.0 NIL
実行時間は計測できないほど速くなりました。ただし、これは枝刈りがうまくいった場合であり、データによっては枝刈りが機能しない場合もありえます。たとえば、xs の要素数を 26, 27, 28, 29, 30 に増やして、xs の最後尾の要素から 1 を引いた値 (1- (car (last xs))) を判定してみましょう。結果は次のようになりました。
* (test #'subset-sum1) (1 3 8 21 55 144 377 987 2584 6765 17711 46368 121393) 6.891 (2 5 13 34 89 233 610 1597 4181 10946 28657 75025 196418) 13.922 (1 3 8 21 55 144 377 987 2584 6765 17711 46368 121393 317811) 27.734 (2 5 13 34 89 233 610 1597 4181 10946 28657 75025 196418 514229) 55.438 (1 3 8 21 55 144 377 987 2584 6765 17711 46368 121393 317811 832040) 110.75
実行時間は速くなりません。ところが、xs を降順に並べておくと、実行時間はとても高速になります。
* (test #'subset-sum1) (121393 46368 17711 6765 2584 987 377 144 55 21 8 3 1) 0.0 (196418 75025 28657 10946 4181 1597 610 233 89 34 13 5 2) 0.0 (317811 121393 46368 17711 6765 2584 987 377 144 55 21 8 3 1) 0.0 (514229 196418 75025 28657 10946 4181 1597 610 233 89 34 13 5 2) 0.0 (832040 317811 121393 46368 17711 6765 2584 987 377 144 55 21 8 3 1) 0.0
ちなみに、xs を乱数でシャッフルしてみたところ、実行時間は次のようになりました。
* (test #'subset-sum1) (1 3 17711 55 121393 6765 144 987 21 2584 377 46368 8) 0.453 (5 610 1597 233 2 75025 4181 89 28657 34 196418 13 10946) 11.328 (317811 6765 1 3 55 46368 121393 21 17711 2584 144 377 987 8) 0.172 (610 10946 1597 2 28657 5 75025 4181 196418 89 34 233 514229 13) 2.063 (377 832040 144 8 317811 46368 21 3 6765 2584 55 1 121393 17711 987) 1.531
今回の問題では、大きな値から試した方が枝刈りの効果は高いようです。興味のある方はいろいろ試してみてください。
次は「動的計画法」で部分和問題を解いてみましょう。総和が N となる部分集合があるか判定するだけでよければ、動的計画法で簡単に解くことができます。
部分和問題の場合、要素をひとつずつ追加しながら、総和 N となる部分集合があるか判定します。簡単な例を示しましょう。次の図を見てください。
xs = {2, 3, 5, 8}, N = 10
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
-------------------------------------------------
0: {} ○ × × × × × × × × × ×
1: {2} ○ × ○ × × × × × × × ×
2: {2, 3} ○ × ○ ○ × ○ × × × × ×
3: {2, 3, 5} ○ × ○ ○ × ○ × ○ ○ × ○
4: {2, 3, 5, 8} ○ × ○ ○ × ○ × ○ ○ × ○
上図は 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 となる部分集合があるか判定します。
xs = {2,3,5,8}, N = 14
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10 11 12 13 14
----------------------------------------------------------
0: {} ○ × × × × × × × × × × × × × ×
1: {2} ○ × ○ × × × × × × × × × × × ×
2: {2,3} ○ × ○ ○ × ○ × × × × × × × × ×
3: {2,3,5} ○ × ○ ○ × ○ × ○ ○ × ○ × × × ×
4: {2,3,5,8} ○ × ○ ○ × ○ × ○ ○ × ○ ○ × ○ ×
3 番目で○の位置は 0, 2, 3, 5, 7, 8, 10 です。次は 8 を追加しますが、総和 14 より大きい値は不要なので、8, 10, 11, 13 の位置に○を追加します。14 の位置は×なので、総和が 14 となる部分集合は無いことがわかります。
プログラムは次のようになります。
リスト : 部分和問題 (動的計画法)
(defun subset-sum-dp (xs n)
(let ((table (make-array (1+ n) :initial-element nil)))
(setf (aref table 0) t)
(dolist (x xs (aref table n))
(do ((y (- n x) (1- y)))
((< y 0))
(when (aref table y)
(setf (aref table (+ x y)) t))))))
○を t で、×を nil で表します。配列をひとつで済ますため、配列の後ろから t の位置を検索していることに注意してください。また、検索の開始位置を n - x とすることで、t をセットするときの範囲チェックを省略しています。今回のプログラムでは xs の要素をすべてチェックしていますが、x + y が n と等しくなったら return-from で t を返してもかまいません。あとは特に難しいところはないと思います。
それでは実際に試してみましょう。数列 xs の要素数は 26, 27, 28, 29, 30 で、(総和 - 1) があるか判定します。
* (test #'subset-sum-dp) T 0.172 T 0.265 T 0.454 T 0.796 T 1.344 NIL
集合の要素数を M, 総和を N とすると、今回のプログラムの実行速度は N * M に比例します。たとえば、N の値を (1- (car (last xs))) とすると、実行結果は次のようになります。
* (test #'subset-sum-dp) T 0.062 T 0.094 T 0.172 T 0.281 T 0.453 NIL
N が小さくなったので、実行時間も速くなりました。このように、動的計画法では N が大きくなると、どうしても時間がかかるようになります。ご注意くださいませ。
もうひとつ簡単な例題として「ナップザック問題」を取り上げます。ナップザック (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 := (name size price)
リスト : 品物のアクセス関数 (defun get-name (item) (car item)) (defun get-size (item) (cadr item)) (defun get-price (item) (caddr item))
リストの先頭要素が名前、 2 番目の要素が大きさで、3 番目の要素が金額を表します。そして、アクセス関数 get-name, get-size, get-price を用意します。
次は、ナップザック問題の解を求める関数 knapsack を作ります。
リスト : ナップザック問題の解法
(defun knapsack (item-list knap-size)
(let ((gain (make-array (1+ knap-size) :initial-element 0))
(choice (make-array (1+ knap-size) :initial-element nil)))
(dolist (item item-list)
(do ((i (get-size item) (1+ i))
(j 0 (1+ j)))
((< knap-size i))
(let ((new-price (+ (get-price item) (aref gain j))))
(when (< (aref gain i) new-price)
(setf (aref gain i) new-price
(aref choice i) item)))))
(print-answer choice knap-size)))
引数 item-list は品物を格納したリスト、knap-size はナップザックの大きさです。最初に、金額を表す配列 gain と選択した品物を格納する配列 choice を用意します。gain は 0 で、choice は nil で初期化します。次の dolist で、品物を一つずつ item-list から取り出します。その次の do ループで、品物のサイズから knap-size まで gain と choice を更新していきます。
変数 j は変数 i から item の大きさを引いた値になります。品物を追加した場合の金額は gain の j 番目の要素に item の金額を足したものになります。この値を new-value にセットします。new-value が gain の i 番目の要素より大きくなれば、gain と choice を更新します。
最後にナップザックの中身を表示する関数 print-answer を作ります。
リスト : 解の表示
; 選択した品物を集める
(defun collect-item (choice knap-size)
(do ((a nil)
(i knap-size (- i (get-size (aref choice i)))))
((null (aref choice i)) a)
(push (aref choice i) a)))
; 品物の個数をカウントする
(defun count-item (xs)
(let ((ys nil))
(dolist (x xs ys)
(let ((y (assoc x ys)))
(if y
(incf (cdr y))
(push (cons x 1) ys))))))
; 解を表示する
(defun print-answer (choice knap-size)
(let ((size 0) (price 0))
(dolist (x (count-item (collect-item choice knap-size)))
(format t "~A, ~D~%" (car x) (cdr x))
(incf size (* (get-size (car x)) (cdr x)))
(incf price (* (get-price (car x)) (cdr x))))
(format t "size : ~D, price : ~D~%" size price)))
関数 collect-item は選んだ品物を集めてリストに格納して返します。choice の knap-size 番目の要素から順番にたどっていき、値が nil になったら終了です。なお、サイズの合計は kanp-size に一致するとは限りません。knap-size よりも小さな値になる場合もあります。関数 count-item は品物の個数を数えてリストに格納して返します。関数 print-answer は品物の個数、サイズと金額の合計値を表示します。
それでは実行結果を示します。
* (knapsack '((a 4 6) (b 3 4) (c 1 1)) 10) (B 3 4), 2 (A 4 6), 1 size : 10, price : 14 NIL
正常に動作していますね。興味のある方はいろいろなデータで試してみてください。