演算法導論學習筆記——動態規劃

本文系轉載，原文地址：http://www.cppblog.com/Fox/archive/2008/05/07/Dynamic_programming.html

以前在學習非數值演算法的時候，曾經瞭解過動態規划算法（Dynamic programming），以下是對Wikipedia上動態規劃的翻譯，圖也是Wikipedia上的，倉促行文，不到之處，請方家指正。

這篇文章的術語實在是太多了，所以我在文中加入了少量注釋，一律以粗斜體註明。

本文的不足之處將隨時修正，MIT的《Introduction to Algorithms》第15章是專門講動態規劃的。

_____________________________________________________________

動態規劃

在數學與電腦科學領域，動態規劃用於解決那些可分解為重複子問題（overlapping subproblems，想想遞迴求階乘吧）並具有最優子結構（optimal
substructure，想想最短路徑演算法）（如下所述）的問題，動態規劃比通常演算法花費更少時間。

上世紀40年代，Richard Bellman最早使用動態規劃這一概念表述通過遍曆尋找最優決策解問題的求解過程。1953年，Richard Bellman將動態規劃賦予現代意義，該領域被IEEE納入系統分析和工程中。為紀念Bellman的貢獻，動態規劃的核心方程被命名為貝爾曼方程，該方程以遞迴形式重申了一個最佳化問題。

在“動態規劃”（dynamic programming）一詞中，programming與“電腦編程”（computer programming）中的programming並無關聯，而是來自“數學規劃”（mathematical programming），也稱最佳化。因此，規劃是指對組建活動的最佳化策略。舉個例子，編製一場展覽的議程可稱為規劃。
在此意義上，規劃意味著找到一個可行的活動計劃。

• 概述

 

圖1 使用最優子結構尋找最短路徑：直線表示邊，波狀線表示兩頂點間的最短路徑（路徑中其他節點未顯示）；粗線表示從起點到終點的最短路徑。

不難看出，start到goal的最短路徑由start的相鄰節點到goal的最短路徑及start到其相鄰節點的成本決定。

 

 

最優子結構即可用來尋找整個問題最優解的子問題的最優解。舉例來說，尋找圖上某頂點到終點的最短路徑，可先計算該頂點所有相鄰頂點至終點的最短路徑，然後以此來選擇最佳整體路徑，如圖1所示。

一般而言，最優子結構通過如下三個步驟解決問題：

a) 將問題分解成較小的子問題；

b) 通過遞迴使用這三個步驟求出子問題的最優解；

c) 使用這些最優解構造初始問題的最優解。

子問題的求解是通過不斷劃分為更小的子問題實現的，直至我們可以在常數時間內求解。

 

 

圖2 Fibonacci序列的子問題：使用有向非循環圖（DAG, directed acyclic graph）而非樹表示重複子問題的分解。

為什麼是DAG而不是樹呢？答案就是，如果是樹的話，會有很多重複計算，下面有相關的解釋。

 

 

一個問題可劃分為重複子問題是指通過相同的子問題可以解決不同的較大問題。例如，在Fibonacci序列中，F3 = F1 + F2和F4 = F2 + F3都包含計算F2。由於計算F5需要計算F3和F4，一個比較笨的計算F5的方法可能會重複計算F2兩次甚至兩次以上。這一點對所有重複子問題都適用：愚蠢的做法可能會為重複計算已經解決的最優子問題的解而浪費時間。

為避免重複計算，可將已經得到的子問題的解儲存起來，當我們要解決相同的子問題時，重用即可。該方法即所謂的緩衝（memoization，而不是儲存memorization，雖然這個詞亦適合，姑且這麼叫吧，這個單詞太難翻譯了，簡直就是可意會不可言傳，其意義是沒計算過則計算，計算過則儲存）。當我們確信將不會再需要某一解時，可以將其拋棄，以節省空間的。在某些情況下，我們甚至可以提前計算出那些將來會用到的子問題的解。

總括而言，動態規劃利用：

1)
重複子問題

2)
最優子結構

3) 緩衝

動態規劃通常採用以下兩種方式中的一種兩個辦法：

自頂向下：將問題劃分為若干子問題，求解這些子問題並儲存結果以免重複計算。該方法將遞迴和緩衝結合在一起。

自下而上：先行求解所有可能用到的子問題，然後用其構造更大問題的解。該方法在節省堆棧空間和減少函數調用數量上略有優勢，但有時想找出給定問題的所有子問題並不那麼直觀。

為了提高按名傳遞（call-by-name，這一機制與按需傳遞call-by-need相關，複習一下參數傳遞的各種規則吧，簡單說一下，按名傳遞允許改變實參值）的效率，一些程式設計語言將函數的傳回值“自動”緩衝在函數的特定參數集合中。一些語言將這一特性儘可能簡化（如Scheme、Common
Lisp和Perl），也有一些語言需要進行特殊擴充（如C++，C++中使用的是按值傳遞和按引用傳遞，因此C++中本無自動緩衝機制，需自行實現，具體實現的一個例子是Automated Memoization
in C++）。無論如何，只有指稱透明（referentially transparent，指稱透明是指在程式中使用運算式、函數本身或以其值替換對程式結果沒有任何影響）函數才具有這一特性。

• 例子

1. Fibonacci序列

尋找Fibonacci序列中第n個數，基於其數學定義的直接實現：

   function fib(n)
       if n = 0
           return 0
       else if n = 1
           return 1
       return fib(n-1) + fib(n-2)

如果我們調用fib(5)，將產生一棵對於同一值重複計算多次的調用樹：

1. fib(5)
2. fib(4) + fib(3)
3. (fib(3) + fib(2)) + (fib(2) + fib(1))
4. ((fib(2) + fib(1)) + (fib(1) + fib(0))) + ((fib(1) + fib(0)) + fib(1))
5. (((fib(1) + fib(0)) + fib(1)) + (fib(1) + fib(0))) + ((fib(1) + fib(0)) + fib(1))

特別是，fib(2)計算了3次。在更大規模的例子中，還有更多fib的值被重複計算，將消耗指數級時間。

現在，假設我們有一個簡單的映射（map）對象m，為每一個計算過的fib及其傳回值建立映射，修改上面的函數fib，使用並不斷更新m。新的函數將只需O(n)的時間，而非指數時間：

   var m := map(0 → 1, 1 → 1)
   function fib(n)
       if map m does not contain key n
           m[n] := fib(n-1) + fib(n-2)
       return m[n]

這一儲存已計算出的數值的技術即被稱為緩衝，這兒使用的是自頂向下的方法：先將問題劃分為若干子問題，然後計算和儲存值。

在自下而上的方法中，我們先計算較小的fib，然後基於其計算更大的fib。這種方法也只花費線性（O(n)）時間，因為它包含一個n-1次的迴圈。然而，這一方法只需要常數（O(1)）的空間，相反，自頂向下的方法則需要O(n)的空間來儲存映射關係。

   function fib(n)
       var previousFib := 0, currentFib := 1
       if n = 0
           return 0
       else if n = 1
           return 1
       repeat n-1 times
           var newFib := previousFib + currentFib
           previousFib := currentFib
           currentFib  := newFib
       return currentFib

在這兩個例子，我們都只計算fib(2)一次，然後用它來計算fib(3)和fib(4)，而不是每次都重新計算。

2. 一種平衡的0-1矩陣

考慮n*n矩陣的賦值問題：只能賦0和1，n為偶數，使每一行和列均含n/2個0及n/2個1。例如，當n=4時，兩種可能的方案是：

+ - - - - +             + - - - - +
| 0 1 0 1 |             | 0 0 1 1 |
| 1 0 1 0 |             | 0 0 1 1 |
| 0 1 0 1 |             | 1 1 0 0 |
| 1 0 1 0 |             | 1 1 0 0 |
+ - - - - +             + - - - - +

問：對於給定n，共有多少種不同的賦值方案。

至少有三種可能的演算法來解決這一問題：窮舉法（brute force）、回溯法（backtracking）及動態規劃（dynamic programming）。窮舉法列舉所有賦值方案，並逐一找出滿足平衡條件的方案。由於共有C(n,
n/2)^n種方案（在一行中，含n/2個0及n/2個1的組合數為C(n,n/2)，相當於從n個位置中選取n/2個位置置0，剩下的自然是1），當n=6時，窮舉法就已經幾乎不可行了。回溯法先將矩陣中部分元素置為0或1，然後檢查每一行和列中未被賦值的元素並賦值，使其滿足每一行和列中0和1的數量均為n/2。回溯法比窮舉法更加巧妙一些，但仍需遍曆所有解才能確定解的數目，可以看到，當n=8時，該題解的數目已經高達116963796250。動態規劃則無需遍曆所有解便可確定解的數目（意思是劃分子問題後，可有效避免若干子問題的重複計算）。

通過動態規劃求解該問題出乎意料的簡單。考慮每一行恰含n/2個0和n/2個1的k*n(1<=k<=n)的子矩陣，函數f根據每一行的可能的賦值對應為一個向量，每個向量由n個整數對構成。向量每一列對應的一個整數對中的兩個整數分別表示該列上該行以下已經放置的0和1的數量。該問題即轉化為尋找f((n/2,n/2),(n/2,n/2),...,(n/2,n/2))（具有n個參數或者說是一個含n個元素的向量）的值。其子問題的構造過程如下：

1) 最上面一行（第k行）具有C(n, n/2)種賦值；

2) 根據最上面一行中每一列的賦值情況（為0或1），將其對應整數對中相應的元素值減1；

3) 如果任一整數對中的任一元素為負，則該賦值非法，不能成為正確解；

4) 否則，完成對k*n的子矩陣中最上面一行的賦值，取k=k-1，計算剩餘的(k-1)*n的子矩陣的賦值；

5) 基本情況是一個1*n的細小的子問題，此時，該子問題的解的數量為0或1，取決於其向量是否是n/2個(0, 1)和n/2個(1, 0)的排列。

例如，在上面給出的兩種方案中，向量序列為：

((2, 2) (2, 2) (2, 2) (2, 2))       ((2, 2) (2, 2) (2, 2) (2, 2))     k = 4
  0      1      0      1              0       0       1      1

((1, 2) (2, 1) (1, 2) (2, 1))       ((1, 2) (1, 2) (2, 1) (2, 1))     k = 3
  1      0      1      0              0      0      1      1

((1, 1) (1, 1) (1, 1) (1, 1))       ((0, 2) (0, 2) (2, 0) (2, 0))     k = 2
  0      1      0      1              1      1      0      0

((0, 1) (1, 0) (0, 1) (1, 0))       ((0, 1) (0, 1) (1, 0) (1, 0))     k = 1
  1      0      1      0              1      1      0      0

((0, 0) (0, 0) (0, 0) (0, 0))       ((0, 0) (0, 0), (0, 0) (0, 0))

動態規劃在此的意義在於避免了相同f的重複計算，更進一步的，上面著色的兩個f，雖然對應向量不同，但f的值是相同的，想想為什麼吧 。

該問題解的數量（序列a058527在OEIS）是1, 2, 90, 297200, 116963796250, 6736218287430460752, ...

下面的外部連結中包含回溯法的Perl原始碼實現，以及動態規劃法的MAPLE和C語言的實現。

3. 棋盤

考慮n*n的棋盤及成本函數C(i,j)，該函數返回方格(i,j)相關的成本。以5*5的棋盤為例：

5 | 6 7 4 7 8
4 | 7 6 1 1 4
3 | 3 5 7 8 2
2 | 2 6 7 0 2
1 | 7 3 5 6 1
- + - - - - -
  | 1 2 3 4 5

可以看到：C(1,3)=5

從棋盤的任一方格的第一階（即行）開始，尋找到達最後一階的最短路徑（使所有經過的方格的成本之和最小），假定只允許向左對角、右對角或垂直移動一格。

5 |
4 |
3 |
2 |   x x x
1 |     o
- + - - - - -
  | 1 2 3 4 5

該問題展示了最優子結構。即整個問題的全域解依賴於子問題的解。定義函數q(i,j)，令：q(i,j)表示到達方格(i,j)的最低成本。

如果我們可以求出第n階所有方格的q(i,j)值，取其最小值並逆向該路徑即可得到最短路徑。

記q(i,j)為方格(i,j)至其下三個方格（(i-1,j-1)、(i-1,j)、(i-1,j+1)）最低成本與c(i,j)之和，例如：

5 |
4 |     A
3 |   B C D
2 |
1 |
- + - - - - -
  | 1 2 3 4 5

q(A) = min(q(B),q(C),q(D)) + c(A)

定義q(i,j)的一般形式：

            |-  inf.                                                  j<1 or
j>n
q(i,j) = -+-  c(i,j)                                                i=1
            |-  min(q(i-1,j-1),q(i-1,j),q(i-1,j+1))+c(i,j)   otherwise.

方程的第一行是為了保證遞迴可以退出（處理邊界時只需調用一次遞迴函式）。第二行是第一階的取值，作為計算的起點。第三行的遞迴是演算法的重要組成部分，與例子A、B、C、D類似。從該定義我們可以直接給出計算q(i,j)的簡單的遞迴代碼。在下面的虛擬碼中，n表示棋盤的維數，C(i,j)是成本函數，min()返回一組數的最小值：

function minCost(i, j)
    if j < 1 or j > n
        return infinity
    else if i = 1
        return c(i,j)
    else
        return min(minCost(i-1,j-1),minCost(i-1,j),minCost(i-1,j+1))+c(i,j)

需要指出的是，minCost只計算路徑成本，並不是最終的實際路徑，二者相去不遠。與Fibonacci數相似，由於花費大量時間重複計算相同的最短路徑，這一方式慢的恐怖。不過，如果採用自下而上法，使用二維數組q[i,j]代替函數minCost，將使計算過程快得多。我們為什麼要這樣做呢？選擇儲存值顯然比使用函數重複計算相同路徑要簡單的多。

我們還需要知道實際路徑。路徑問題，我們可以通過另一個前任數組p[i,j]解決。這個數組用於描述路徑，代碼如下：

function computeShortestPathArrays()
     for x from 1 to n
         q[1, x] := c(1, x)
     for y from 1 to n
         q[y, 0]     := infinity
         q[y, n + 1] := infinity
     for y from 2 to n
         for x from 1 to n
             m := min(q[y-1, x-1], q[y-1, x], q[y-1, x+1])
             q[y, x] := m + c(y, x)
             if m = q[y-1, x-1]
                 p[y, x] := -1
             else if m = q[y-1, x]
                 p[y, x] :=  0
             else
                 p[y, x] :=  1

剩下的求最小值和輸出就比較簡單了：

function computeShortestPath()
     computeShortestPathArrays()
     minIndex := 1
     min := q[n, 1]
     for i from 2 to n
         if q[n, i] < min
             minIndex := i
             min := q[n, i]
     printPath(n, minIndex)

function printPath(y, x)
     print(x)
     print("<-")
     if y = 2
         print(x + p[y, x])
     else
         printPath(y-1, x + p[y, x])

4. 序列比對

序列比對是動態規劃的一個重要應用。序列比對問題通常是使用編輯操作（替換、插入、刪除一個要素等）進行序列轉換。每次操作對應不同成本，目標是找到編輯序列的最低成本。

可以很自然地想到使用遞迴解決這個問題，序列A到B的最優編輯通過以下措施之一實現：

插入B的第一個字元，對A和B的剩餘序列進行最優比對；

刪去A的第一個字元，對A和B進行最優比對；

用B的第一個字元替換A的第一個字元，對A的剩餘序列和B進行最優比對。

局部比對可在矩陣中列表表示，單元(i,j)表示A[1..i]到b[1..j]最優比對的成本。單元(i,j)的成本計算可通過累加相鄰單元的操作成本並選擇最優解實現。至於序列比對的不同實現演算法，參見Smith-Waterman和Needleman-Wunsch。

對序列比對的話題並不熟悉，更多的話也無從談起，有熟悉的朋友倒是可以介紹一下。

• 應用動態規劃的演算法

1) 許多字串操作演算法如最長公用子列、最長遞增子列、最長公用字串；

2) 將動態規劃用於圖的樹分解，可以有效解決有界樹寬圖的產生樹等許多與圖相關的演算法問題；

3) 決定是否及如何可以通過某一特定上下文無關文法產生給定字串的Cocke-Younger-Kasami (CYK)演算法；

4) 電腦國際象棋中轉換表和駁斥表的使用；

5) Viterbi演算法（用於隱式馬爾可夫模型）；

6) Earley演算法（一類圖表分析器）；

7)
Needleman-Wunsch及其他生物資訊學中使用的演算法，包括序列比對、結構比對、RNA結構預測；

8)
Levenshtein距離（編輯距離）；

9)
弗洛伊德最短路徑演算法；

10)
連鎖矩陣乘法次序最佳化；

11)
子集求和、背包問題和分治問題的偽多項式時間演算法；

12) 計算兩個時間序列全域距離的動態時間規整演算法；

13) 關係型資料庫的查詢最佳化的Selinger（又名System R）演算法；

14) 評價B樣條曲線的De Boor演算法；

15) 用於解決板球運動中斷問題的Duckworth-Lewis方法；

16) 價值迭代法求解馬爾可夫決策過程；

17) 一些圖形映像邊緣以下的選擇方法，如“磁鐵”選擇工具在Photoshop；

18)
間隔調度；

19) 自動換行；

20)
巡迴旅行商問題（又稱郵差問題或貨擔郎問題）；

21)
分段最小二乘法；

22)
音樂資訊檢索跟蹤。

對於這些演算法應用，大多未曾接觸，甚至術語翻譯的都有問題，鑒於本文主要在於介紹動態規劃，所以倉促之中，未及查證。

• 相關

1) 貝爾曼方程

2)
馬爾可夫決策過程

3) 貪心演算法

• 參考

Adda, Jerome, and Cooper, Russell, 2003. Dynamic Economics. MIT Press. An accessible introduction to dynamic programming in economics. The link contains sample programs.

Richard Bellman, 1957, Dynamic Programming, Princeton University Press. Dover paperback edition (2003),
ISBN 0486428095.

Bertsekas, D. P., 2000. Dynamic Programming and Optimal Control, Vols. 1 & 2, 2nd ed. Athena Scientific.
ISBN 1-886529-09-4.

Thomas H. Cormen,
Charles E. Leiserson,
Ronald L. Rivest, and
Clifford Stein, 2001.
Introduction to Algorithms, 2nd ed. MIT Press & McGraw-Hill.
ISBN 0-262-03293-7. Especially pp. 323–69.

Giegerich, R., Meyer, C., and Steffen, P., 2004, "A Discipline of Dynamic Programming over Sequence Data,"
Science of Computer Programming 51: 215-263.

Nancy Stokey, and
Robert E. Lucas, with
Edward Prescott, 1989. Recursive Methods in Economic Dynamics. Harvard Univ. Press.

S. P. Meyn, 2007.
Control Techniques for Complex Networks, Cambridge University Press, 2007.

• 外部連結

• Dyna, a declarative programming language for dynamic programming algorithms
• Wagner, David B., 1995, "Dynamic Programming." An introductory article on dynamic programming in
  Mathematica.
• Ohio State University: CIS 680: class notes on dynamic programming, by Eitan M. Gurari
• A Tutorial on Dynamic programming
• MIT course on algorithms - Includes a video lecture on DP along with lecture notes -- See lecture 15.
• More DP Notes
• King, Ian, 2002 (1987), "A Simple Introduction to Dynamic Programming in Macroeconomic Models." An introduction to dynamic programming as an important tool
  in economic theory.
• Dynamic Programming: from novice to advanced A TopCoder.com article by Dumitru on Dynamic Programming
• Algebraic Dynamic Programming - a formalized framework for dynamic programming, including an
  entry-level course to DP, University of Bielefeld
• Dreyfus, Stuart, "Richard Bellman on the birth of Dynamic Programming."
• Dynamic programming tutorial
• An Introduction to Dynamic Programming

_____________________________________________________________

關於動態規劃，這隻是一篇譯文，後面將根據實際問題具體寫點動態規劃的應用。