Dynamic Programming

長江後浪催前浪,一替新人趲舊人。《張協狀元》

資之深,則取之左右逢其原。《孟子》

Dynamic Programming

先透過一個簡單的例子,感受一下「動態規劃」吧!

範例:階乘(Factorial)

1 × 2 × 3 × ⋯ × N。整數1到N的連乘積。N階乘。N!。

N!源自(N-1)!,如此就遞迴分割問題了。

陣列的每一格對應每一個問題。設定第一格的答案,再以迴圈依序計算其餘答案。

UVa 623 568 10220 10323

時間複雜度

總共N個問題,每個問題花費O(1)時間,總共花費O(N)時間。

空間複雜度

求1!到N!:總共N個問題,用一條N格陣列儲存全部問題的答案,空間複雜度O(N)。

只求N!:用一個變數累計乘積,空間複雜度O(1)。

Dynamic Programming - recurrence

Dynamic Programming
= Divide-and-Conquer Method + Memoization

動態規劃是分治法的延伸。當分治法分割出來的問題,一而再、再而三出現,就運用記憶法儲存這些問題的答案,避免重複求解,以空間換取時間。

動態規劃的過程,就是反覆地讀取數據、計算數據、儲存數據。

1. 把原問題遞迴分割成許多更小的問題。(recurrence)
   1-1. 子問題與原問題的求解方式皆類似。(optimal sub-structure)
   1-2. 子問題會一而再、再而三的出現。(overlapping sub-problems)
2. 設計計算過程:
   2-1. 確認每個問題需要哪些子問題來計算答案。(recurrence)
   2-2. 確認總共有哪些問題。(state space)
   2-3. 把問題一一對應到表格。(lookup table)
   2-4. 決定問題的計算順序。(computational sequence)
   2-5. 確認初始值、計算範圍。(initial states / boundary)
3. 實作,主要有兩種方式:
   3-1. Top-down
   3-2. Bottom-up

1. recurrence

遞迴分割問題時,當子問題與原問題完全相同,只有數值範圍不同,我們稱此現象為recurrence,再度出現、一再出現之意。

【註:recursion和recurrence,中文都翻譯為「遞迴」,然而兩者意義大不相同,讀者切莫混淆。】

此處以爬樓梯問題當作範例。先前於遞歸法章節,已經談過如何求踏法,而此處要談如何求踏法數目。

踏上第五階,只能從第四階或從第三階踏過去。因此「爬到五階」源自兩個子問題:「爬到四階」與「爬到三階」。

「爬到五階」的踏法數目,就是總合「爬到四階」與「爬到三階」的踏法數目。寫成數學式子是「f(5) = f(4) + f(3)」,其中「f(‧)」表示「爬到某階之踏法數目」。

依樣畫葫蘆,得到「f(4) = f(3) + f(2)」、「f(3) = f(2) + f(1)」。

「爬到兩階」與「爬到一階」無法再分割、沒有子問題,直接得到「f(2) = 2」、「f(1) = 1」。

整理成一道簡明扼要的遞迴公式:

f(n) =
 ⎧ 1                , if n = 1
 ⎨ 2                , if n = 2
 ⎩ f(n-1) + f(n-2)  , if n >= 3 and n <= 5

爬到任何一階的踏法數目,都可以藉由這道遞迴公式求得。n代入實際數值,遞迴計算即可。

為什麼分割問題之後,就容易計算答案呢?因為分割問題時,同時也分類了這個問題的所有可能答案。分類使得答案的規律變得單純,於是更容易求得答案。

2-1. recurrence

f(n) =
 ⎧ 1                , if n = 1
 ⎨ 2                , if n = 2
 ⎩ f(n-1) + f(n-2)  , if n >= 3

2-2. state space

想要計算第五階的踏法數目。

全部的問題是「爬到一階」、「爬到二階」、「爬到三階」、「爬到四階」、「爬到五階」。

至於「爬到零階」、「爬到負一階」、「爬到負二階」以及「爬到六階」、「爬到七階」沒有必要計算。

2-3. lookup table

建立六格的陣列,儲存五個問題的答案。

表格的第零格不使用,第一格是「爬到一階」的答案,第二格是「爬到二階」的答案,以此類推。

如果只計算「爬完五階」,也可以建立三個變數交替使用。

2-4. computational sequence

因為每個問題都依賴「階數少一階」、「階數少二階」這兩個問題,所以必須由階數小的問題開始計算。

計算順序是「爬到一階」、「爬到二階」、……、「爬到五階」。

2-5. initial states / boundary

最先計算的問題是「爬到一階」與「爬到二階」,必須預先將答案填入表格、寫入程式碼,才能繼續計算其他問題。心算求得「爬到一階」的答案是1,「爬到二階」的答案是2。

最後計算的問題是原問題「爬到五階」。

為了讓表格更順暢、為了讓程式碼更漂亮,可以加入「爬到零階」的答案,對應到表格的第零格。「爬到零階」的答案,可以運用「爬到一階」的答案與「爬到兩階」的答案,刻意逆推而得。

最後可以把初始值、尚待計算的部份、不需計算的部分,統整成一道遞迴公式:

f(n) =
 ⎧ 0                , if n < 0               [Exterior]
 ⎪ 1                , if n = 0               [Initial State]
 ⎨ 1                , if n = 1               [Initial State]
 ⎪ f(n-1) + f(n-2)  , if n >= 2 and n <= 5   [Computation]
 ⎩ 0                , if n > 5               [Exterior]

UVa 11069

3. 實作

直接用遞迴實作,而不使用記憶體儲存各個問題的答案,是最直接的方式,也是最慢的方式。時間複雜度O(f(n))。問題一而再、再而三的出現,不斷呼叫同樣的函式求解,效率不彰。剛接觸DP的新手常犯這種錯誤。

正確的DP,是一邊計算,一邊將計算出來的數值存入表格,以後便不必重算。這裡整理了兩種實作方式,各有優缺點:

1. Top-down
2. Bottom-up

3-1. Top-down

這個實作方式的好處是不必斤斤計較計算順序,因為程式碼中的遞迴結構會迫使最小的問題先被計算。這個實作方式的另一個好處是只計算必要的問題,而不必計算所有可能的問題。

這個實作方式的壞處是程式碼採用遞迴結構,不斷呼叫函式,執行效率較差。這個實作方式的另一個壞處是無法自由地控制計算順序,因而無法妥善運用記憶體,浪費了可回收再利用的記憶體。

UVa 10285 10446 10520

3-2. Bottom-up

訂定一個計算順序,然後由最小的問題開始計算。特色是程式碼通常只有幾個迴圈。這個實作方式的好處與壞處,與前一個方式恰好互補。

首先建立表格。

心算「爬到零階」的答案、「爬到一階」的答案,填入表格當中,作為初始值。分別填到表格的第零格、第一格。

尚待計算的部份就是「爬到兩階」的答案、……、「爬到五階」的答案。通常是使用迴圈,按照計算順序來計算。

計算過程的實作方式,有兩種迥異的風格。一種是「往回取值」,常見的實作方式。

另一種是「往後補值」,罕見的實作方式。

計算完畢之後,最後印出答案。

UVa 495 900 10334

總結

第一。先找到原問題和其子問題們之間的關係,寫出遞迴公式。如此一來,便可利用遞迴公式,用子子問題的答案,求出子問題的答案;用子問題的答案,求出原問題的答案。

第二。確認可能出現的問題全部總共有哪些,這樣才能知道要計算哪些問題,才能知道總共花多少時間、多少記憶體。

第三。有了遞迴公式之後,就必須安排出一套計算的順序。大問題的答案,總是以小問題的答案來求得的,所以,小問題的答案是必須先算的,否則大問題的答案從何而來呢?

一個好的安排方式,不但使程式碼容易撰寫,還可重複利用記憶體空間。

第四。記得先將最小、最先被計算的問題,心算出答案,儲存入表格,內建於程式碼之中。一道遞迴公式必須擁有初始值,才有辦法計算其他項。

第五。實作DP的程式時,會建立一個表格,在表格存入所有大小問題的答案。安排好每個問題的答案在表格的哪個位置,這樣計算時才能知道該在哪裡取值。

切勿存取超出表格的元素,產生溢位情形,導致答案算錯。計算過程當中,一旦某個問題的答案出錯,就會如骨牌效應般一個影響一個,造成很難除錯。

Dynamic Programming - counting / optimization

範例:樓梯路線(Staircase Walk),計數問題

一個方格棋盤,從左上角走到右下角,每次只能往右走一格或者往下走一格。請問有幾種走法?

對於任何一個方格來說,只可能「從上走來」或者「從左走來」,答案是兩者相加。

「從上走來」是一個規模更小的問題,「從左走來」是一個規模更小的問題,答案是兩者相加。

二維陣列的每一格對應每一個問題。設定第零行、第零列的答案,再以迴圈依序計算其餘答案。

時間複雜度分析:令H和W是棋盤長寬。計算一個問題需要O(1)時間(兩個子問題答案相加的時間)。總共HW個問題,所以計算所有問題需要O(HW)時間。

空間複雜度分析:總共HW個問題,所以需要O(HW)空間,簡單來說就是二維陣列啦!如果不需要儲存所有問題的答案,只想要得到其中一個特定問題的答案,那只需要一維陣列就夠了,也就是O(min(H,W))空間。

節省記憶體是動態規劃的重要課題!

如果只打算求出一個問題,那麼只需要儲存最近算出來的問題答案,讓計算過程可以順利進行就可以了。

兩條陣列輪替使用,就足夠儲存最近算出來的問題答案、避免f[i-1][j]超出陣列範圍。

事實上,一條陣列就夠了。也不能再少了。

如果某些格子上有障礙物呢?把此格設為零。

如果也可以往右下斜角走呢?添加來源f[i-1][j-1]。

如果可以往上下左右走呢?不斷繞圈子,永遠不會結束,走法無限多種。

UVa 10599 825 926

遞迴公式

若瞭解遞迴關係,就不必強記遞迴公式。若瞭解圖片意義,就不必強記數學符號。

f(i, j) =
 ⎧ 0                      , if i < 0 or j < 0     [Exterior]
 ⎪ 1                      , if i = 0 or j = 0     [Initial States]
 ⎨ f(i-1, j) + f(i, j-1)  , if i > 0 and i < 8    [Computation]
 ⎪                          and j > 0 and j < 8
 ⎩ 0                      , if i >= 8 or j >= 8   [Exterior]

f(i, j):從格子 (0, 0) 走到格子 (i, j) 的走法數目。

遞歸方向

這個問題雙向都可以遞歸。對於任何一個方格來說,只可能「向右走出」或者「向下走出」。

範例:樓梯路線(Staircase Walk),極值問題

動態規劃的問題,可以分為「計數問題」和「極值問題」。方才介紹「計數問題」,現在介紹「極值問題」。

一個方格棋盤,格子擁有數字。從左上角走到右下角,每次只能往右走一格或者往下走一格。請問總和最小的走法?(或者總和最大的走法?)

想要印出路線,另外用一個陣列,記錄從哪走來。

額外介紹一個技巧。為了避免減一超出邊界,需要添補許多程式碼。整個棋盤往右下移動一格,就能精簡許多程式碼。

範例:樓梯路線(Staircase Walk),極值問題

節省記憶體是動態規劃當中重要的課題!

方才的分割方式:分割最後一步,窮舉最後一步從哪走來;方才的實作方式:由小到大的迴圈。問題答案f[i][j],可以精簡成一維陣列。路線來源p[i][j],無法精簡成一維陣列。

想讓路線來源精簡成一維陣列,必須採用另一種分割方式:從地圖中線分割,窮舉穿過中線的所有地點;同時採用另一種實作方式:由大到小的遞迴。

Dynamic Programming - state / DAG

state / DAG

以State和DAG的觀點,重新看待動態規劃。

動態規劃得類比成「狀態State」:「問題」變「狀態」,「全部問題」變「狀態空間」,「遞迴關係」變「狀態轉移函式」。

動態規劃得類比成「有向無環圖DAG」:既然遞迴關係不能循環,顯然就是DAG。「問題」變「點」,「遞迴關係」變「邊」,「計算順序」變「拓撲順序」。

即便讀者不懂State和DAG也沒關係,只要抓住兩個要點:每個小問題各是一個狀態,只有數值範圍不同;狀態之間是單行道,依序求解,不能循環。

範例:爬樓梯(Climbing Stairs)

f(n) =
 ⎧ f(n-1) + f(n-2)  , if n > 2
 ⎨ 2                , if n = 2
 ⎩ 1                , if n = 1

一維狀態。兩個子問題。

時間複雜度O(N)。空間複雜度O(1)。

範例:河內塔(Tower of Hanoi)

f(n) =
 ⎰ f(n-1) + 1 + f(n-1)  , if n > 1
 ⎱ 1                    , if n = 1

除了子問題,還附帶其他成本。

時間複雜度O(N)。空間複雜度O(1)。

範例:樓梯路線(Staircase Walk)

f(h, w) =
 ⎰ f(h-1, w) + f(h, w-1)  , if h > 1 and w > 1
 ⎱ 1                      , if h = 1 or w = 1

二維狀態。兩個子問題。

時間複雜度O(HW)。空間複雜度O(HW)。

表格可以重複利用。空間複雜度降為O(W)甚至O(min(H,W))。

範例:組合(Combination)

f(n, m) =
 ⎰ f(n-1, m-1) + f(n-1, m)  , if n > 1 and m > 1 and m < n
 ⎱ 1                        , if n = 1 or n = m

二維狀態。兩個子問題。

時間複雜度O(N²)。空間複雜度O(N²)。

🡲

組合數表格,即是巴斯卡三角形。巴斯卡三角形靠左對齊,視覺上就可以一一對應至表格。巴斯卡三角形左右對稱,對稱部分可以精簡。

表格可以只用一半。空間複雜度仍然相同。但是沒有意義。個人電腦記憶體充足,節省一半差別不大。

表格可以重複利用。空間複雜度降為O(N)。但是沒有意義。與其重複利用表格,不如直接使用公式解。

範例:矩陣相乘次序(Matrix Chain Multiplication)

一連串矩陣相乘,無論從何處開始相乘,計算結果都一樣。

然而,計算時間卻有差異。兩個矩陣,大小為a x b及b x c,相乘需要O(abc)時間(當然還可以更快,但是此處不討論)。那麼一連串矩陣相乘,最少需要多少時間呢?

一連串矩陣,從最後一次相乘的地方分開,化作兩串矩陣相乘。考慮所有可能的分開方式。

f(i, k) =  min  { f(i, j) + f(j+1, k) + r[i] ⋅ c[j] ⋅ c[k] }
         j=i⋯k-1

f(i, k):從第i個矩陣乘到第k個矩陣,最少的相乘次數。
r[i]:第i個矩陣的row數目。
c[i]:第i個矩陣的column數目。

時間複雜度O(N³)。空間複雜度O(N²)。

計算順序可以調整成online版本。

同類型的動態規劃問題(只是順手整理,不是本節重點):

Matrix Chain Multiplication
Optimal Binary Search Tree
Hu–Tucker Compression
Minimum Weight Triangulation of Convex Polygon
Cocke–Younger–Kasami Algorithm

UVa 348 442

範例:Longest Increasing Subsequence

一連串數字,從左到右挑選數字,越來越大。最多幾個?

一連串數字,保留開頭幾個數字,隔離最後一個數字,化作一小串數字與一個新數字。考慮所有可能的保留方式。

f(j) = max ( 1, max { f(i) + ( a[i] < a[j] ? 1 : -∞ ) } )
              i=0⋯j-1

f(j):前j個數字可以挑出的數字數量,而且第j個數字被挑出。
a[j]:第j個數字。

遞迴公式內含判斷式:嘗試讓數字更大且更多。

正確答案必須搜尋表格:枚舉每個問題,其答案最小者。

時間複雜度O(N²)。空間複雜度O(N)。

同類型的動態規劃問題(只是順手整理,不是本節重點):

Longest Increasing Subsequence
Longest Common Subsequence
Longest Palindromic Subsequence

UVa 10900

設計狀態

設計狀態,目前沒有標準流程。變化多端、層出不窮。

正因如此,演算法競賽經常出現動態規劃問題,以便考驗選手思維。也因如此,選手經常探討如何設計狀態,以便應對各種動態規劃問題。關於設計狀態,也許大家可以請教競賽選手。競賽選手的網路文章、專家學者的學術論文,兩者相比,前者比較平易近人。

以下姑且提供一些個人經驗,即便我不擅長設計狀態。

設計狀態:狀態變數

簡易策略:所有未知變因,通通納入狀態變數。

一般來說,狀態變數越多,表格維度越大,計算時間越久。一般來說,狀態變數越少越好。

變數不一定是表面上看到的變因。例如「人潮最多的時段」,遞增的標的,是時刻,而不是訪客。

變數也有優先順序。例如「儲存座標」,變數有「點編號」和「座標軸編號」,優先順序導致表格造型改變、遞迴公式改變。

變數甚至可以是答案。答案納入狀態變數,等同於枚舉各種可能的答案。宛如試誤法。

變數甚至可以跟答案缺乏直接關聯。答案需要從表格抽取數字、重新計算而得。宛如預計算。

設計狀態:遞歸策略

變數減一。也可能減更多數值。

也可能需要枚舉許多個數值。

也可能需要複雜的公式。

也可能仰賴其他變數。

也可能包含判斷式。

Dynamic Programming - bitset

bitset

bitset是一個二進位數字,每一個bit分別代表每一件東西,1代表開啟,0代表關閉。例如現在有十個燈泡,編號設定為零到九,其中第零個、第三個、第九個燈泡是亮的,剩下來的燈泡是暗的。我們用一個10 bit的二進位數字1000001001,表示這十個燈泡的亮暗狀態。

建立一個大小為2¹⁰的陣列,便囊括了所有可能的狀態。陣列的每一格,就代表一種燈泡開關的狀態。

當狀態數量呈指數成長,可以利用bitset作為狀態。

UVa 10952

範例:Maximum Matching

以線相連的兩物,可以配對在一起。求最大配對數目暨配對方式。

Maximum Matching」有多項式時間演算法,可是很難實作;動態規劃雖然慢了些,是指數時間演算法,但是容易實作。移除匹配成一對的點,就得到遞迴公式。

f(S+{i}+{j}) = max { f(S) + adj(i, j) }
              i,j∉S
f(S) = max { f(S-{i}-{j}) + adj(i, j) }
      i,j∈S

使用bitset,已匹配標成1,未匹配標成0。時間複雜度O(2ᴺ N²),空間複雜度O(2ᴺ)。

這個演算法可以再修正,讓時間複雜度降為O(2ᴺ N)。各位可以試試看。

這個演算法需要大量記憶體,所以無法計算N很大的情況,何況編譯器也不准我們建立太大的陣列,N=28就是極限了。這個演算法同時也需要大量時間,以現在的個人電腦來說,N=17就已經要花上幾分鐘才能求出答案了。

UVa 10888 10911 11439 10296 11156

範例:Hamilton Path

找到一條路徑,剛好每一個點都去過一次。有可能找不到。

Hamilton Path」尚無多項式時間演算法。直覺的解法是backtracking,窮舉所有點的各種排列方式,一種排列方式當作一條路徑,判斷是不是Hamilton Path。

運用動態規劃,可以減少計算時間。拆掉一條路徑的最後一條邊,就得到遞迴公式。需要額外維度,記錄路徑終點。

f(S+{j}, j) = or { f(S, i) and adj(i, j) }   where j∉S
              i∈S
f(S, j) = or { f(S-{j}, i) and adj(i, j) }   where j∈S
          i∈S
          i≠j

時間複雜度O(2ᴺ N²),空間複雜度O(2ᴺ)。

UVa 216 10068 10496 10818 10937 10944 10605 10890 265

範例:骨牌平舖(Domino Tiling)

中文網路稱作「插头DP」或「轮廓线DP」。

UVa 11741

範例:子集子集和(Subset Sum of Subset)

英文網路稱作「Sum over Subsets DP」。

範例:不重複路線(Self-avoiding Walk),計數問題

先前介紹過樓梯路線(Staircase Walk)。樓梯路線問題,只能往兩個方向走,可以簡單的遞迴分割,得到多項式時間演算法。不重複路線問題,可以往四個方向走,無法簡單的遞迴分割,只有指數時間演算法。儘管如此,不重複路線還是可以使用動態規劃。

UVa 10572 10531

範例:不重複路線(Self-avoiding Walk),極值問題

極值問題反而存在多項式時間演算法!原問題化作圖:點有權重,邊無權重。原問題化作最短路徑:點權重挪至入邊權重。

請見本站文件「Shorest Path」。

Dynamic Programming Speedup

Dynamic Programming Speedup

極值問題有許多手法可以降低時間複雜度。根據recurrence的樣式,採取對應的手法。

詳細來說,一個問題(狀態)通常有許多個子問題。

計數問題:每一個子問題都納入統計。窮舉全部。

極值問題:找到最好的那一個子問題。挑選一個。

極值問題,只需要找到一個子問題,其答案最小(最大)。所有子問題,依照答案大小排序,就能快速找到那一個子問題,而不必逐一檢查。根據recurrence的樣式,採取對應的排序手法。

即時排序

詳細來說,極值問題當中,一個問題通常有許多個子問題(候選解),從中找到最好的那一個子問題(最佳解)。

正常做法:枚舉所有候選解,找出一個最佳解。

取巧做法:即時排序候選解,只枚舉一部分的候選解。

巧妙地安排計算順序,方便即時排序候選解。每個問題的候選解,彼此之間有關聯。下一個問題的候選解,繼承了上一個問題的候選解。不必重頭排序,只需微調順序。

候選解數值、候選解位置

詳細來說,即時排序的對象,可以分為兩種:候選解數值(樓梯路線的f陣列)、候選解位置(樓梯路線的p陣列)。

候選解數值:大家已經發現許多經典手法,例如stack/deque、envelope/convex hull、unimodal function、Monge matrix。

候選解位置:手法極少,例如Knuth–Yao speedup。

Dynamic Programming Speedup - stack / deque

概論

F[n] = min/max { F[i] ⋅ W[i] + C[i] }
       i=0⋯n-1

F[n]是未知數,F[i] W[i] C[i]是已知數。計算到F[n]時,F[i]早已計算完畢,因此F[i]是已知數。

這種形式的recurrence,直接計算是O(N²)。此處介紹更快的演算法,降為O(N)。

stack

括號配對極值。stack保持嚴格遞增(嚴格遞減),即時移除多餘數值,即時獲取過往最大值(最小值)。

          for finding maximum,
          keep monotone increasing  ---->
      ┌──────┬──────┬──────┬─────────┬───────╯
stack │ M[2] │ M[4] │ M[6] │   ...   │ M[10]
      └──────┴──────┴──────┴─────────┴───────╮
                                      ^^^^^^^ extract maximum
               clean tails until monotonicity
               then push new M[i] into tail

M[i] = F[i] ⋅ W[i] + C[i]塞入尾端。塞入前,先彈出尾端數值們,使得M[i]塞入尾端之後,stack嚴格遞增。最大值從尾端取得。

範例:Largest Empty Rectangle

詳見「Largest Empty Rectangle」。

範例:All Nearest Smaller Values

deque

滑動視窗極值。deque保持嚴格遞增(嚴格遞減),即時移除多餘數值,即時獲取過往最小值(最大值)。

           for finding minimum,
           keep monotone increasing  ---->
      ╰──────┬──────┬──────┬─────────┬───────╯
deque   M[2] │ M[4] │ M[6] │   ...   │ M[10]
      ╭──────┴──────┴──────┴─────────┴───────╮
       ^^^^^^                         ^^^^^^^
       extract minimum      clean tails until monotonicity
                            then push new M[i] into tail

M[i] = F[i] ⋅ W[i] + C[i]塞入尾端。塞入前,先彈出尾端數值們,使得M[i]塞入尾端之後,deque嚴格遞增。塞入後,再彈出頭端數值們,使得deque符合滑動視窗範圍。最小值從頭端取得。

中文網路稱作「单调队列优化」。

範例:Maxium Sum Subarray

詳見「Maxium Sum Subarray」。

範例:Maximum Average Subarray

詳見「Maxium Average Subarray」。

由於斜率是關鍵,因此中文網路稱作「斜率优化」。

Dynamic Programming Speedup - envelope / convex hull

概論

F[n] = min/max { F[i] ⋅ W[n] + C[i] }
       i=0⋯n-1

W[i]換成W[n]。處理方式大幅改變。

這種形式的recurrence,直接計算是O(N²)。此處介紹更快的演算法,降為O(Nlog²N)。

先備知識是「Convex Hull」、「Point-Line Duality」、「Fenchel Duality」。

envelope

直線y = F[i] x + C[i],鉛直線x = W[n],交點Y軸座標F[n]。

最小(大)值對應最低(高)交點。

最小(大)值位於下(上)包絡線。

convex hull

如果討厭包絡線,可透過點線對偶,從包絡線變成凸包。

直線穿過點(F[i], -C[i]),斜率W[n],Y軸截距-F[n]。

垂直方向翻轉,讓最小(大)值對應最低(高)截距。

直線穿過點(F[i], C[i]),斜率-W[n],Y軸截距F[n]。

最小(大)值位於下(上)凸包的切線,跟Y軸的交點。

如果不熟悉點線對偶,可透過移項推導,得到相同結論:

概論

根據F[i]與W[n]的性質,時間複雜度隨之變化:

一、W[n]皆相同:不必維護凸包,只需維護每個點在「斜率-W[n]的直線的垂直方向上」的先後順序,彷彿先前章節「deque加速」。一個常見的例子是W[n] = 1。O(N)。

中文網路稱作「单调队列优化」。

二、F[i]與W[n]皆單調:Andrew's monotone chain維護下凸包(求最小值時)。切線斜率-W[n]遞減/增。直接從上次切點開始,往左/右找到新切點。O(N)。

由於斜率遞增,中文網路也將此手法稱作「斜率优化」。

三、F[i]單調:Andrew's monotone chain 維護下凸包(求最小值時)。切線斜率-W[n]會變。三分搜尋找到切點,或者二元搜尋凸包斜率找到切點。O(NlogN)。

中文網路稱作「WQS二分」。

四、無特別性質:動態凸包資料結構。O(Nlog²N)。

中文網路稱作「李超线段树」。

UVa 12524

範例:1D p-Median Problem

詳見「p-Median Problem」。

範例:Bounded Knapsack Problem

詳見「Bounded Knapsack Problem」。

Dynamic Programming Speedup - unimodal function

概論

F[n] = min/max { M[i] }
       i=0⋯n-1

M[i] is monotone/unimodal

根據M[i]的性質,時間複雜度隨之變化:

一、F[n]的子問題們,其答案M[0]...M[n-1]恰好是單調函數:根本沒啥好算,最佳解顯然是第一個(最後一個)子問題。O(N)。

二、F[n]的子問題們,其答案M[0]...M[n-1]恰好是單峰函數:三分搜尋山峰,或者二元搜尋斜率。O(NlogN)。

三、每個問題的單峰函數,山峰位置恰好遞增(往右移動):用同一條掃描線尋找山峰。O(N)。

概論

什麼時候會是單峰函數呢?舉兩個例子。

F[n] =  min  { max(F[i], G[i]) + 5 }
      i=0⋯n-1

(1) F[i] is increasing
(2) G[i] is decreasing
(=>) max(F[i], G[i]) + 5 is unimodal

F遞增、G遞減,則max(F[i], G[i]) + 5是單峰函數。不過這個例子有點蠢,山谷恰好永遠不動。

F[n] =  min  { F[i] + G[i] }
      i=0⋯n-1

(1) F[i] is convex
(2) G[i] is convex
(=>) F[i] + G[i] is convex (also unimodal)

F G皆是凸函數,則F[i] + G[i]是凸函數,即是單峰函數。不過這個例子也有點蠢,山谷恰好永遠不動。

範例:Egg Drop

一堆蛋,已知耐力皆相同,不知耐力為多少。試求耐力。

耐力是以樓層衡量:大於某樓層,摔下必破,無法重複使用;小於等於某樓層,摔下必不破,完全不會折損,可以重複使用。

這個問題有許多變形,此處討論實驗場地是n層樓的情況:

一、最少摔破幾顆蛋?需要事先準備幾顆蛋?

答:一顆蛋。從一樓開始摔,逐步上樓,直到摔破為止。

二、無限多蛋,運氣不好時,最少摔幾次?

答:二元搜尋。F[0] = 0, F[1] = 1, F[n] = F[ceil((n-1)/2)] + 1。

三、一顆蛋,運氣不好時,最少摔幾次?

答:n次。耐力不幸是n樓,從一樓開始摔,要摔n次。

四、兩顆蛋,運氣不好時,最少摔幾次?

答:首發選在i樓。如果摔破了,剩下一顆蛋,只好從一樓開始測試;最差的情況,耐力是i-1樓,要摔i-1次。如果沒有摔破,仍是兩顆蛋,問題仍相同,範圍縮小成n-i樓。兩種情況取最大值。窮舉i,找到次數最少者。

F[n] = 1 +  min  { max(i-1, F[n-i]) }
           i=1⋯n
F[n] = 1 +  min  { max(i, F[n-1-i]) }   調整索引值
          i=0⋯n-1

i遞增,F[n-1-i]遞減,因此max(i, F[n-1-i])是單峰函數,而且山谷持續往右移動,可以使用掃描線解決。另外也有數學公式解:

五、k顆蛋,運氣不好時,最少摔幾次?

答:留給讀者。

UVa 10934 882

範例:Isotonic Regression

英文網路稱作「Slope Trick」。

範例:Word Wrap

詳見「Word Wrap」。

Dynamic Programming Speedup - Monge matrix

概論

F[n] = min/max { M[i][n] }
       i=0⋯n-1

M[i][n] is upper-triangular monotone/totally monotone/Monge

根據M[i][n]的性質,時間複雜度隨之變化:

一、M[i][n]是上三角monotone matrix:似乎沒有特別快的演算法,仍是O(N²)。

二、M[i][n]是上三角totally monotone matrix:兩個演算法,O(NlogN)與O(Nα(N))。此處只談第一個。

三、M[i][n]是上三角Monge matrix:至今仍然沒有專屬演算法。大家直接沿用totally monotone matrix的演算法。

三種矩陣皆細分成凹凸兩種情況。兩種情況的演算法稍微不同。此處均介紹。

三種矩陣當中,以Monge matrix最為常見。它也是性質最強的矩陣,若三成立則二成立,若二成立則一成立。

本文當中,「上三角Monge matrix」並非是「下三角恰為0的Monge matrix」,而是「僅上三角符合Monge matrix」。

Monge matrix的定義,在古代曾經稱作quadrangle inequality,導致中文網路稱作「四边形不等式优化」。

先備知識是「Pairwise Distance」。

概論

F[n] = min/max { F[i] ⋅ W[i] + C[i][n] }
       i=0⋯n-1

(1) C[i][n] is upper-triangular Monge
(=>) M[i][n] = F[i] ⋅ W[i] + C[i][n] is upper-triangular Monge

Monge matrix橫條加上同一數,仍是Monge matrix。藉由這個數學性質,我們得以處理這種形式的recurrence。

這種形式的recurrence,直接計算是 O(N²)。此處介紹更快的演算法,降為O(NlogN)。

convex totally monotone matrix

大意:直條最小值位置往上遞減。使用stack。

首先畫出C[i][n],一個上三角矩陣。圖片省略了實際數值。

計算F[1]:F[0] ⋅ W[0]加到C[0][1]即得。

計算F[2]:F[0] ⋅ W[0]加到C[0][2],F[1] ⋅ W[1]加到C[1][2],取直條C[:][2]最小值。

以此類推。時間複雜度O(N²)。

重新反省計算過程,改良演算法:

一、算出F[i]之後,F[i] ⋅ W[i]加到對應的橫條,比較省事。

二、計算F[i],即是取直條最小值。隨時記錄每個直條的最小值位置,每次F[i] ⋅ W[i]加到對應的橫條,順手更新最小值位置。

得到一個更容易解釋的演算法。時間複雜度仍是O(N²)。

凸全單調矩陣(直條版本),每個子矩陣都是凸單調矩陣,直條最小值位置總是往上遞減!

算出F[i]之後,從左往右更新最小值位置。最小值位置可能變i、可能不變。變與不變有著唯一一條分界線,左側一律變i,右側一律不變,以滿足遞減性質。

雖然可以提早break,但是時間複雜度仍是O(N²)。

位置相同者,合併成一個區間,最多N個區間。每當更新最小值位置,從左往右判斷每個區間的右端位置,直到遭遇不變的位置。然後二元搜尋該區間,找到分界線。

使用stack實作,一筆資料有兩個參數:區間右端的直條編號、最小值位置。每回合pop一些區間、做一次二元搜尋、push一個區間。時間複雜度O(NlogN)。

concave totally monotone matrix

大意:直條最小值位置往下遞增。使用deque。

改成從右往左更新最小值位置。時間複雜度O(NlogN)。

Dynamic Programming Speedup - interval

概論

最佳解位置,恰好已排序。

每個問題,依照大小排序。最佳解位置,恰好也依照大小排序。

一個問題不必枚舉所有候選解。參考先前問題的最佳解位置,一個問題只需要枚舉一部分的候選解。

中文網路稱作「决策单调性」。

Knuth–Yao speedup

F[i,j] =  min  { F[i,k] + F[k+1,j] + C[i,j] }
        k=i⋯j-1

(1) C[i,j] is Monge*:  C[a,b] + C[a+1,b+1] ≤ C[a,b+1] + C[a+1,b]
(2) C[i,j] is sorted:  C[a,b] ≤ C[a,b+1] and C[a,b] ≤ C[a-1,b]
                       (C[a,b] ≤ C[c,d] when [a,b] ⊆ [c,d])
(=>) F[i,j] is Monge*
(=>) P[i,j] is sorted: P[a,b] ≤ P[a,b+1] and P[a,b] ≤ P[a+1,b]
                       (P[a,b-1] ≤ P[a,b] ≤ P[a+1,b])

*non-negative upper-trianglar concave Monge
(1) C[i,j] is Monge:   ↖ + ↘ ≤ ↗ + ↙       (concave)
(2) C[i,j] is sorted:  ← ≤ → and ↓ ≤ ↑     (toward ↗)
                       (←↙↓ ≤ →↗↑)
(=>) F[i,j] is Monge
(=>) P[i,j] is sorted: ← ≤ → and ↑ ≤ ↓     (toward ↘)
                       (↖ ≤ ↗ ≤ ↘)

這種形式的recurrence,直接計算是O(N³)。此處介紹更快的演算法,降為O(N²)。

附帶成本C[i,j]既是Monge matrix又是sorted matrix,導致最佳解數值F[i,j]是Monge matrix,也導致最佳解位置P[i,j]是sorted matrix。【待補證明】

一個問題的最佳解位置,範圍限縮。原本是a ≤ P[a,b] ≤ b,現在是P[a,b-1] ≤ P[a,b] ≤ P[a+1,b]。Monge matrix夾住左邊界,sorted matrix夾住右邊界。候選解變少,時間複雜度變低。

過程當中使用了Monge matrix,導致中文網路也將此手法稱作「四边形不等式优化」。

英文網路稱作「Knuth's Optimization」。

範例:Optimal Binary Search Tree

N筆資料,欲建立成一棵「Binary Search Tree」。並且預測了每筆資料的搜尋次數。

請問Binary Search Tree是什麼形狀,才能讓拜訪到的節點數量最少呢?也就是說,每個節點的「深度」乘上「搜尋次數」,總和要最小。

遞迴公式類似於Matrix Chain Multiplication,都是記錄區間。窮舉樹根,分割成左右兩棵子樹遞迴下去。問題總共O(N²)個,一個問題要窮舉O(N)種分割點,時間複雜度O(N³)。

由於第二層迴圈的c(i,k)維持定值,不會影響最大值的判斷結果,所以可以挪到迴圈外面,減少加法次數,減少執行時間。

總共O(N²)個問題,每個問題必須窮舉O(N)個分割點,時間複雜度O(N³)。

到這裡都和普通的Dynamic Programming沒兩樣。接下來要更進一步。

每次計算一個問題,總是得窮舉所有的分割點。然而有些分割點顯然是錯誤的,尤其是靠近區間邊界的那些分割點,實在不太可能將兩棵子樹分割的夠均勻、令總和最小。

相近的問題,其分割點也很相近。問題[a,b],嘗試從從右端拿掉一筆資料,成為問題[a,b-1]。問題[a,b]、問題[a,b-1]的分割點很相近。

為了讓左右子樹均勻,[a,b]的分割點一定要大於等於[a,b-1]的分割點,才能降低總和。小於[a,b-1]的分割點,沒有窮舉的必要,樹只會越不平衡、總和只會更大不會更小!

問題[a+1,b]的情況也十分類似,不再贅述。

也就是說,問題[a,b]的分割點,必定位於更小的問題[a,b-1]和[a+1,b]的分割點之間。計算一個問題,大可不必窮舉所有的分割點。

觀察分割點表格,[a,b-1]是左方格子,[a+1,b]是下方格子。要計算一個分割點,窮舉範圍就是左方格子的值到下方格子的值。也就是說每一個格子都會大於等於左方格子、小於等於下方格子。

每一條左上右下斜線,左上最小值是0,右下最大值是N-1,每一條斜線最多窮舉2N = O(N)個分割點。

除了初始值之外,總共N-1條斜線,需要窮舉的分割點總共O(N²)個,所以時間複雜度降為O(N²)。

UVa 10304 10003 12057 12809

Hu–Shing speedup

F[i,j] =  min  { F[i,k] + F[k+1,j] + C[i,j,k] }
        k=i⋯j-1

這種形式的recurrence,直接計算是O(N³)。此處介紹更快的演算法,降為O(NlogN)。

範例:Matrix Chain Multiplication

我沒有研究。論文上百頁,我不想讀。

【尚無正式名稱】

F[i][j] =  min  { F[i-1][k-1] + C[k,j] }
          k=1⋯j

(1) C[i,j] is Monge: C[a,b] + C[a+1,b+1] ≤ C[a,b+1] + C[a+1,b]
(=>) P[i][j] is row-sorted: P[i,j] ≤ P[i,j+1]

(2) C(|i-j|) is convex
(=>) F[i][j] is row-sorted
(=>) F[i-1][k-1] + C[k,j] is unimodal

這種形式的recurrence,直接計算是O(NM²)。此處介紹更快的演算法,降為O(NM)。

一、附帶成本C[i,j]恰是Monge matrix,導致最佳解位置P[i][j]是row-sorted matrix。【尚待確認】

一個問題的最佳解位置,範圍限縮。Monge matrix夾住左邊界。因為沒有夾住右邊界,所以用二分法處理右邊界。O(NMlogM)。

英文網路稱作「Divide-and-Conquer Optimization」。

二、附帶成本C(|i - j|)恰是convex function,再導致最佳解數值F[i][j]是row-sorted matrix,也導致候選解數值F[i][k] + C[k,j]是unimodal function。【尚待確認】

一個問題的最佳解位置,範圍限縮。Monge matrix夾住左邊界。convex function夾住右邊界。O(NM)。

三、不管C[i,j]是什麼了。先前章節「convex hull加速」實施N次。W[n] = 1,彷彿「deque加速」。O(NM)。

範例:1D p-Center Problem

詳見「p-Center Problem」。

參考文獻

[envelope/convex hull] [unimodal function]
using geometric techniques to improve dynamic programming algorithms for the economic lot-sizing problem and extensions

[Monge matrix]
dynamic programming with convexity, concavity and sparsity
an almost linear time algorithm for generalized matrix searching

[Knuth–Yao speedup]
optimum binary search trees
efficient dynamic programming using quadrangle inequalities

[Hu–Shing speedup]
computation of matrix chain products, part I, part II
revisiting “computation of matrix chain products”