推 hcsoso: Nice summary. 推! 03/18 03:32
推 cutekid: 推(Y)，法布施! 03/18 13:45
推 pigalan: 大推~~~ 03/20 22:58
推 wolfpig: 推你 03/29 15:09
跟大家分享最近研究 monotonicity、1D/1D DP、convex hull trick 的心得。
我想很多人高中的時候就搞懂了，這篇文章的主要貢獻大概就是文獻整理吧。
我每次看這種介紹都會被搞昏，因為有一大堆索引值。
所以我盡量的把索引值的大小關係寫清楚，
一般情況下的使用會是 i < j < k < k'。
當要比較兩個陣列元素的時候，我會盡量把索引值小的放左側，
也就是我會盡量寫 A[i] ?? A[j] 。
Dominatation
先從一個簡單的例子開始介紹 dominate 的概念。
給定一個長度為 n 的整數陣列 A ，對所有 k 計算 A[0] ~ A[k - 1]
中的最小值，亦即計算 B[k] = min_{i < k} A[i]
首先考慮兩個元素 A[i], A[j], i < j (A[i] 在 A[j] 之前)，
如果 A[i] > A[j] 的話，那對於所有 j < k, B[k] 不可能會是 A[i] 。
換句話說 A[i] 被 A[j] dominate 了，A[i] 對於 j 之後的元素是沒有
影響的，因為選擇 A[j] 必定是個更好的選擇，所以不需要保留關於
A[i] 的資訊。
反之，如果 A[i] < A[j] ，那 A[i] 就 dominate 了 A[j] 。
因此解法很單純，只要從 k = 0 ~ n-1 一次掃描，紀錄當前的最小值即可。
再考慮一個類似的問題 (all nearest smaller values)
https://en.wikipedia.org/wiki/All_nearest_smaller_values
給定一個長度為 n 的整數陣列 A ，對於所有 k 計算出現在 A[k]
之前，最靠近 A[k] 且比 A[k] 小的數字。
亦即計算 B[k] = A[ max_{i < k : A[i] < A[k]} ]
因為必須要找到最靠近的元素，所以問題就變得比較複雜。
考慮兩個元素 A[i], A[j], i < j (A[i] 在 A[j] 之前)，
A[i] > A[j] 的情況跟之前一樣，但是A[i] < A[j] 的情形就不同了，
因為 A[j] 還是有可能會是 B[k], j < k 的解答。
所以可以把所有沒被 dominate 的元素儲存起來，計算 B[k] 的時候，
只需要從這些沒被 dominate 的元素裡面找解答就好。
這概念跟 Pareto frontier 的概念一樣，只從比較好的候選人中挑解答。
但是 Pareto frontier 可能還是包含了許多元素，所以從 Pareto frontier
裡面挑解答還是很花時間，某些時候是可以把 Pareto frontier 裡面的元素
放在資料結構內，使得解答可以快速的找出來。
此時就需要引入另外一個概念 -- monotonicity。
Monotonicity
在 all nearest smaller values 問題中，令 i1, ..., it 表示 在 A[k] 前
的 Pareto frontier 中元素的索引值，我們可以把在 Pareto frontier 的元
素按照索引值排列
A[i1], A[i2], ..., A[it], i1 < i2 < ... it < k
那麼 這些元素的值必定是遞增的，也就是
A[i1] < A[i2] < .... A[it]
因為如果索引值大的元素有比較低的值，那索引值低的元素就
會被 dominate 了。
而且 all nearest smaller values 問題要求要找到最靠近且小於 A[k]
的元素，所以如果 A[it] < A[k] ，那麼 B[k] 必為 A[it] 。
如果 A[it] > A[k] ，那麼 A[it] 就被 A[k] dominate 了。
所以只要依照 A[it], A[it-1], ... 的順序去找，找到 A[ih] < A[k] ，
那麼 B[k] 必為 A[ih] 。
然後可以把 A[k] 放到 Pareto frontier 中，然後接著計算 B[k+1]。
因為所有元素的操作都是在 Pareto frontier 的尾端，所以可以使用
一個 stack 來儲存 Pareto frontier 。
上面這個例子是使用 stack ，這邊我再介紹一個使用 queue 來儲存
Pareto frontier 的例子：
給定一個長度為 n 的整數陣列 A 和一正整數 L ，計算所有
長度為 L 的子陣列中的最小值。
亦即計算 B[k] = min_{k - L + 1 <= i <= k} A[i]
簡單來說，這問題就是要找出 sliding window 內的最小值。
所以我們想要設計一個 queue ，同時支援最小值查詢。
(題外話， deque + 最小值查詢也是可以辦得到的
可以參考 Gajewska 和 Tarjan 的論文 Deques with heap order
http://www.sciencedirect.com/science/article/pii/0020019086900281 )
考慮兩個元素 A[i] 和 A[j], i < j, (A[i] 在 A[j] 之前)，
如果 A[i] > A[j] ，那麼 A[i] 就不可能會是 B[k], j < k 的答案了，
也就是 A[i] 被 A[j] dominate 了 。
如果是 A[i] < A[j] ，那麼 A[i] 還是有可能是 B[k], j < k 的答案。
所以跟 all nearest smaller values 問題類似，把 Pareto frontier
中的元素按照索引值排列，此時 Pareto frontier 中的元素必遞增。
所以最小值必定是索引值最小的元素。
跟 all nearest smaller values 類似，也會需要把索引值大的元素從
Pareto frontier 中移除，同時也會新增加索引值大的元素。
不同點在於這問題有長度的限制，計算 B[k] 時的 Pareto frontier 必不
能包含 A[i], i < k - L 。
也就是說在維護 Pareto frontier 時會需要把索引值小的元素從
Pareto frontier 中移除。
所以此時需要用一個 deque 來維護 Pareto frontier ，只是索引值小
的那端只有 deque 不會 enqueue ，而索引值大的那端則是 enqueue 和
deque 都會發生。
所以總歸來說，如果要使用 domination 的性質，首先必須要設計一個
domination 的關係，把不可能為答案的 candidate 刪除，這步驟雖然
不容易，但是並不是最困難的地方。
而如果要使用 monotonicity 性質，就必須要考慮以下三點：
1. 如何排序 Pareto frontier
2. 如何從有序的 Pareto frontier 中快速找到解答
3. 如何更新 Pareto frontier 且滿足有序的性質
有興趣的人可以研究 maximum density segment 的問題
給定一個長度為 n 的整數陣列 A 和一正整數 L ，找出長度至少
為 L 的子陣列中平均值最大的子陣列。也就是要找出 i, j 滿足
i + L - 1 <= j 使得 A[i..j] / (j - i + 1) 的值最大
解法可以參考這個 https://goo.gl/FUmnPt 。
這問題的 domination 是要考慮三個點 A[i], A[j], A[k] 的，只考慮
兩個點是找不到 domination 的。而且從 Pareto frontier 找答案的
方式也不太一樣。
如果除了有長度的下限 L 之外，又有長度上限 U 的限制的話，deque
的維護就更複雜了，可以參考
D. T. Lee, Tien-Ching Lin 和 Hsueh-I Lu 的論文
Fast Algorithms for the Density Finding Problem
http://link.springer.com/article/10.1007%2Fs00453-007-9023-8
1D/1D DP
在使用 DP 解決序列上的問題時，1D/1D 是很常見的形式，我接下來
都會用 line wrap 問題當例子。
https://en.wikipedia.org/wiki/Line_wrap_and_word_wrap#Minimum_raggedness
這是 Knuth 在設計 TeX 的時候遇到的問題，這邊只是簡化的版本。
http://onlinelibrary.wiley.com/doi/10.1002/spe.4380111102/abstract
給定一串 n 個單字和一個正整數 L，把這些單字按照分行，
使得每行的字母數與 L 的差距越小越好
正式的定義如下：
給定一正整數 L 和一個長度為 n 的整數陣列 A ，
找出 p (p 是算法自己決定的) 個索引值 k1, ..., kp = n-1
把 A 分成 p 個子陣列 A[0...k1], A[k1+1...k2], ...
且最小化 penalty 總和。
子陣列 A[i...j] 的 penalty 定義為 (L - sum(A[i...j]))^2
令 w[i][j] 表示 A[i...j] 的 penalty 。
令 B[k] 表示 A[0...k] 的最佳分法的 penalty 。
此問題可以利用下列遞迴關係式求解：
B[k] = min_{i < k} B[i] + w[i + 1][k]
時間複雜度是 O(n^2)
為了解釋方便 如果 B[k] = B[i] + w[i + 1][k] ，就稱之為 B[k] 選擇了 A[i] 。
因為 F 是一個一維函數，而 F[k] 考慮的 candidate 也只有一個維度，
所以這類的 DP 就叫做 1D/1D DP 了。
此外還有 2D/0D DP, 2D/1D DP, 2D/2D DP ，只是不在這篇文章討論範圍了。
首先讓我們分析這個問題。
考慮兩個元素 A[i] 和 A[j], i < j, (A[i] 在 A[j] 之前)，
不管是 A[i] > A[j] 還是 A[i] < A[j] ，似乎都沒有很明顯的 domainte
的關係，畢竟問題有本質上的不同。
但是直覺上，因為這問題是在分行，所以儘管在算 B[k] 時，讓 A[i...k]
在一行是好選擇，但是應該會存有一個 k' > k 時 A[j...k'] 會變成比較
好的選擇了。
這性質有人稱為決策單調性，只是這性質應該只對 concave 的 w 有用。
(concave 的定義我就省略了)。
Hirschberg 和 Larmore 首先對這類問題展開研究。
The Least Weight Subsequence Problem
http://epubs.siam.org/doi/abs/10.1137/0216043
他們發現當 w 是 concave 的時候，
對於任何兩個元素 A[i] 和 A[j], i < j, (A[i] 在 A[j] 之前)
如果只考慮選擇 A[i] 或是 A[j] 的話，
會存在一個界線 h ，使得 k <= h 時，B[k] 挑選 A[i]
當 k > h 時，B[k] 挑選 A[j]。
而這個界線是可以在 O(lg n) 的時間內用二分法計算出來的。
所以可以維護一個 deque 來放所有的 candidate ，及其是最佳解的區間
也就是在計算 B[k] 時，如果 deque 是 (A[i1], h1) (A[i2], h2) ...
那麼在只考慮從 A[0] ~ A[k-1] 中挑的話
A[i1] 會是 B[k'] 的最佳選擇，如果 k <= k' <= h1
A[i2] 會是 B[k'] 的最佳選擇，如果 h1+1 <= k' <= h2
等等。
如此一來這問題就可以在 O(n lg n) 的時間內被解決了。
如果想要知道詳細內容的，可以看這篇
Galil 和 Park 的論文
Dynamic programming with convexity, concavity and sparsity
http://www.sciencedirect.com/science/article/pii/0304397592901353
Convex hull trick
第一次出現這技巧應該是在 Wagelmans, Hoesel 和 Kolen 的論文
Economic lot-sizing- An O(n log n) algorithm that runs in linear time
in the Wagner-Whitin case
http://pubsonline.informs.org/doi/abs/10.1287/opre.40.1.S145
只是他們只拿來解特定的問題
在 Hoesel, Wagelmans 和 Moerman 的論文提到了比較一般化的方法
Using geometric techniques to improve dynamic programming algorithms
for the economic lot-sizing problem and extensions
http://www.sciencedirect.com/science/article/pii/0377221794900779
使用幾何的方法來加速 DP ，專門處理
B[k] = C[k] + min_{i < k} (B[i] + D[i] + E[k]F[i])
的 DP 問題， 如果 E 和 F 是 monotone 的話，就可以
在線性時間內找到最佳解。
雖然看起來跟 line wrap 不太像，但是實際上是可以轉化的。
在 line wrap 問題中，遞迴關係式是
B[k] = min_{i < k} B[i] + w[i + 1][k]
兩邊帶入 w 的定義
= min_{i < k} B[i] + (L - sum(A[i+1...k]))^2
令 P 為 prefix sum 陣列， P[i] = sum A[0..i], P[-1] = 0
B[k] = min_{i < k} B[i] + (L - (P[k] + P[i]))^2
展開 (L - (P[k] + P[i]))^2
= min_{i < k} B[i] + L^2 - 2L(P[k] + P[i]) + (P[k] + P[i])^2
展開 (P[k] + P[i])^2
= min_{i < k} B[i] + L^2 - 2L(P[k] + P[i]) + P[k]^2 + 2P[k]P[i] + P[i]^2
把跟 i 無關的提到外面來
= P[k]^2 + L^2 - 2LP[k] + min_{i < k} B[i] + P[i]^2 - 2LP[i] + 2P[k]P[i]
所以可以令 C[k] = P[k]^2 + L^2 - 2LP[k]
D[i] = P[i]^2 - 2LP[i]
E = F = P
就可以使用這個方法了。
這方法的關鍵就是可以把 1D/1D 的問題轉成一個幾何的問題，從而找出
一個 dominate 的關係，然後利用 monotonicity 來加速。
而在 Brucker 的論文中，則使用了純代數的方法。
Efficient algorithms for some path partitioning problems
http://www.sciencedirect.com/science/article/pii/0166218X94001465
(這篇論文的的內容後來出現在IOI 2002 的 Batch Scheduling 問題)
其中最關鍵的一步就是證明了
如果存在有一個 d 函數，和一個遞增的 f 函數
使得 對於所有 i < j < k,
計算 B[k] 時，A[i] 被 A[j] dominate 若且唯若 d(i, j) <= f(k)
那 B 就可以在線性時間內算出來。
套在 line wrap 上面， 考慮 A[i] 和 A[j], i < j, (A[i] 在 A[j] 之前)
在計算 B[k] 時, i < j < k, 如果選擇 A[i] 比選擇 A[j] 差的話代表
B[i] + w[i+1][k] >= B[j] + w[j+1][k]
兩邊帶入 w 的定義
iff B[i] + (L - sum(A[i+1...k]))^2 >= B[j] + (L - sum(A[j+1...k]))^2
利用之前的結論，把兩邊的 P[k]^2 + L^2 - 2LP[k] 都消掉
iff B[i] + P[i]^2 - 2LP[i] + 2P[k]P[i] <= B[j] + P[j]^2 - 2LP[j] + 2P[k]P[j]
把 所有跟 k 相關的移到右邊 跟 k 無關的移到左邊
iff (B[i] + P[i]^2) - (B[j] + P[j]^2) <= 2L(P[i] - P[j]) + 2P[k](P[j]-P[i])
兩邊同除 P[j] - P[i]
iff ((B[j] + P[j]^2) - (B[i] + P[i]^2)) / (P[j] - P[i]) <= 2(P[k] - L)
左邊就是 delta 右邊是 f 了。
其實也不難看得出來
Hoesel, Wagelmans 和 Moerman 的模型
B[k] = C[k] + min_{i < k} (B[i] + D[i] + E[k]F[i])
只要 E 和 F 都是 monotone 的，都可以套到 Brucker 的模型裡。
（但是有些 case <= 要換成 >=）。
簡單來說，看到一個 1D/1D 的問題，可以先試著套用 Brucker 的模型，
也就是在計算 B[k] 的時候，比較 A[i] 和 A[j] 兩種選擇, i < j < k。
如果可以化成 Brucker 要求的形式，那問題就可以線性時間解出來了。
我也嘗試過使用 Brucker 的模型來解釋 maximum density segment ，
不過似乎是沒辦法直接套用的，因為需要兩個點 i, k 才能 dominate 另一個
點 j ， Brucker 的模型都是只考慮用一個點就可以 dominate 另一個點的。
結語
其實還有一個 DP 的最佳化技巧叫做 Knuth-Yao quadrangle inequality，
可以拿來加速 optimal binary search tree 的建立。
https://en.wikipedia.org/wiki/Optimal_binary_search_tree
經過專家的研究發現，其實 Knuth-Yao quadrangle inequality 實際上
也可以用 total monotonicity 來解釋的。
Bein, Golin, Larmore and Zhang
Quadrangle-Inequality Speedup is a Consequence of Total Monotonicity
http://dl.acm.org/citation.cfm?id=1644032
至於 total monoonicity, Monge, convex/concave 之間的關係可以參考
http://www.csie.ntnu.edu.tw/~u91029/DynamicProgramming.html#7
參考文獻
想要了解各種 DP 的最佳化技巧可以參考
Bein 的論文
Advanced Techniques for Dynamic Programming
http://link.springer.com/referenceworkentry/10.1007/978-1-4419-7997-1_28
考古
當 w 是 concave 時 Wilber 設計出了線性演算法，但是是基於 SMAWK 的。
The concave least-weight subsequence problem revisited
http://www.sciencedirect.com/science/article/pii/0196677488900326
https://en.wikipedia.org/wiki/SMAWK_algorithm
Galil 和 Park 考慮 B[k] = min_{i < k} C[i] + w[i + 1][k] 的情形也
設計出了線性演算法，還是基於 SMAWK 。
A linear-time algorithm for concave one-dimensional dynamic programming
http://www.sciencedirect.com/science/article/pii/002001909090215J
然後 Larmore 和 Schieber 設計了不基於 SMAWK 的線性演算法
On-line dynamic programming with applications to the prediction of
RNA secondary structure
http://www.sciencedirect.com/science/article/pii/019667749190016R
Klawe 也設計了線性演算法，但是只有發表在 technical report 上。
可以在這篇論文裡面找到 pseudo code
Consecutive interval query and dynamic programming on intervals
http://www.sciencedirect.com/science/article/pii/S0166218X98000213
當 w 是 convex 的時候（定義我就省略了）
Miller and Myers 設計了一個 O(n lg n) 的方法。
Sequence comparison with concave weighting functions
http://www.sciencedirect.com/science/article/pii/S0092824088800168
Galil and Giancarlo 也設計了一個 O(n lg n) 的方法，
Speeding up dynamic programming with applications to molecular biology
http://www.sciencedirect.com/science/article/pii/0304397589901011
而且他的方法可以解決 B[k] = min_{i < k} C[i] + w[i + 1][k]
形式的問題。
Klawe 和 Kleitman 設計了 O(n α(n)) 的方法。
An Almost Linear Time Algorithm for Generalized Matrix Searching
http://epubs.siam.org/doi/abs/10.1137/0403009
Eppstein 考慮 convex 和 concave 可以 mix 的情況
Sequence comparison with mixed convex and concave costs
http://www.sciencedirect.com/science/article/pii/0196677490900319
Eppstein, Galil, Giancarlo, and Italiano 考慮了 sparse 的情況
Sparse dynamic programming II: convex and concave cost functions
http://dl.acm.org/citation.cfm?id=146656
對 FP 有興趣的人可以看這兩篇
Oege de Moor and Jeremy Gibbons
Bridging the Algorithm Gap: A Linear-Time Functional Program
for Paragraph Formatting
http://www.sciencedirect.com/science/article/pii/S0167642399000052
Sharon A. Curtis and Shin-Cheng Mu
Calculating a linear-time solution to the densest-segment problem
http://dx.doi.org/10.1017/S095679681500026X
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 65.96.6.117
※ 文章網址: https://www.ptt.cc/bbs/Prob_Solve/M.1458211168.A.907.html