PTT數位生活區 / Prob_Solve (計算數學 Problem Solving)

[心得] 1D/1D DP and convex hull trick

看板Prob_Solve (計算數學 Problem Solving)作者FRAXIS (喔喔)時間8年前 (2016/03/17 18:39)推噓5(5推 0噓 0→)

留言5則, 5人參與討論串1/1

跟大家分享最近研究 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