"Mathematical Analysis of Algorithms" 閱讀心得

"Mathematical Analysis of Algorithms" 閱讀心得

“Mathematical Analysis of Algorithms”是著名計算機界大神Knuth在1971年發表的論文。以前只是聽說Knuth非常神，看了這篇論文才體會到Knuth到底有多神…Orz

此外，特別感謝 @solaaaaa 聚聚，沒有他的指導我可能根本看不懂這篇論文…...

這篇文章要解決什么問題？

作為算法分析這一領域的早期論文，這篇文章回答了以下兩個問題：

1. 算法分析的核心目的是什么？

   使用量化方法來衡量算法的“好壞”。

2. 算法分析解決的是哪些類型的問題？

   A. 對一個特定算法的時間復雜度和空間復雜度的分析。

   B. 對解決一類問題的所有算法進行分析，試圖找出“最好”的算法。

此外，本文特別指出，雖然B類的問題看上去比A類的問題更有價值去解決，但對(2)類問題的研究存在兩個重大的問題：1、只要對“最好”的定義稍有改變，一個B類問題就會變成另一個完全不同的(2)類問題；2、B類問題往往過於復雜，難以進行有效的數學建模，而如果建立過於簡化的模型，容易得到不合常理的結果。出於以上原因，只有少數B類問題（如基於比較的排序） 能得到有效的解決，絕大多數算法分析依然是對A類問題的研究。

之后本文的核心，就是通過兩個例子，來展現算法分析的基本思路。

例一：In Situ Permutation

1. 問題定義

給定一個一維數組\(x_1, x_2, …, x_n\)和一個對\(1,2,…,n\)的排列\(p(1),p(2),…,p(n)\)，我們需要將一維數組\(x\)根據函數\(p\)替換為一維數組$x_{p(1)}, x_{p(2)}, …,x_{p(n)} $，有以下要求：

1. 不能修改存儲排列\(p(1),p(2),…,p(n)\)的空間。
2. 整個算法的空間復雜度必須為\(O(1)\)。

2. 設計算法

對以下算法的設計基於以下一個重要的事實：在任意一個排列\(p(1),p(2),…,p(n)\)，我們總會存在若干個“環”，這個環形如\(p(i_1)=i_2,p(i_2)=i_3,…,p(i_k)=i_1\)。學過抽象代數的同學可以對此給出一個嚴格證明。例如對於一個排列：

\[(p(1), p(2),…,p(9))=(8,2,7,1,6,9,3,4,5) \]

通過觀察我們發現了以下環：

\[\begin{cases} p(1)=8, p(8)=4, p(4)=1 \\ p(2)=2 \\ p(3)=7 ,p(7 )=3 \\ p(5)=6,p(6)=9,p(9)=5 \end{cases} \]

我們就可以定義這個環中最小的值為這個環的頭元素。那么只要我們發現了一個環的頭元素\(k\)，就可以將原來的數組中的\(x_k\)空出來，將緊隨其后的元素\(x_{p(k)}\)填入，然后\(x_{p(k)}\)的位置就會被空出來可以填入\(x_{p(p(k))}\)……最后將頭元素\(x_k\)填入環尾部那個排列對應的\(x\)數組的位置即可。

對應的偽代碼描述如下：

for j = 1 to n: # 外循環
    # 判斷 j 是不是環的頭元素
    # 從p(j)開始遍歷這個環
    k = p(j)    
    # 如果j不是環的頭元素，那么就會存在一個環上點k<j
    while k > j: # 循環1 a
        k = p(k)
    if k == j:   # 語塊2 b
        # k就是環的頭元素
        y = x[j], l = p(k)
        while l != j: # 循環3 c 
            x[k] = x[l], k = l, l = p(k)
        x[k] = y
        # 到這一步為止一個環上的元素全部按照排列置換完畢

3. 算法分析

對一個特定算法分析（A類問題）而言，最重要的是證明算法的正確性。在正確的基礎上，我們討論它的最好情況復雜度，最壞情況復雜度和平均情況復雜度。對平均情況的分析往往有更重要的意義。但就算法分析而言，一般對平均情況的分析會復雜得多，最后常常歸結為某個很難的數學問題。

就這個問題而言，為了便於討論，我們將內層循環和語句塊如上面偽代碼注釋所示進行標號$a,b,。

正確性

事實上，算法的設計過程已經簡單說明了這個算法的正確性。要給出一個嚴謹的證明非常麻煩這里不再贅述。

最好情況與最壞情況

顯然，外層循環無論如何都會被執行\(n\)次，因此我們主要考慮內層\(a\)和\(b\)的執行次數。

我們先構造只有一個長\(n\)的環的情況。令

\[(p(1), p(2),…,p(n))=(2,3,...,n,1) \]

那么\(a\)就進入了最壞情況，這樣，對於第\(i\)個外層循環，\(a\)循環需要執行\(n-i\)次，總共執行次數為

\[\sum_{i=1}^n (n-i) = O(n^2) \]

巧妙的是，此時剛好是\(b\)的最好情況！因此\(b\)只用執行一次，這一次執行的復雜度為\(O(n)\)。

因此，這個最壞情況的執行時間為\(O(n^2)\)。

我們再構造有\(n​\)個環的情況。令

\[(p(1),p(2),...,p(n))=(1,2,...,n) \]

這時\(b\)進入了最壞情況，\(a\)進入了最好情況，類似分析可得這種情況的執行時間為\(O(n)\)。

平均情況

我們首先對這個問題進行重新建模。依然是對之前舉過的排列

\[(p(1), p(2),…,p(9))=(8,2,7,1,6,9,3,4,5) \]

我們可以將它的環寫成\((5,6,9)(3,7)(2)(1,8,4)\)，滿足：1、每個環開頭是其最小的元素。2、每個環的頭元素從大到小遞減，那么就可以直接表示成另一個排列\((q(1),q(2),…,q(9))=(5,6,9,3,7,2,1,8,4)\)，因為我們只要找到其中的數\(q_k\)使得\(q_k=min\{q_1,q_2,…,q_k\}\)，那么這個數就是一個環的頭元素，這個環從這個元素開始直到下一個具有相同性質的元素之前為止。比如，由於\(3=min\{5,6,9,3\}\)，我們就發現了\(q(4)=3\)是一個環的頭元素，這個環是\((3,7)\)，因為我們發現了\(7\)之后的\(2\)比前面所有數都小，所以\(2\)就是下一個環的元素了。這樣，我們就建立了從一個排列\(p\)到另一個排列\(q\)的映射。

我們對\(q\)這個排列，首先研究\(a\)循環在平均情況下的執行次數。

對於任意一個\(a\)循環，它從一個隨機的環上的元素開始向后遍歷，也就是說，如果外循環中\(j=q(i)\)，\(k\)就會一直往后進入\(q(i+1),q(i+2),…,\)直到找到一個位置\(q(i+r)\)使得\(q(i+r)<q(i)\)，或者\(q(i)\)就是環的頭元素。所以，當外循環\(j=q(i)\)時，就會從\(q(i)\)到\(q(i+r)\)執行這么多次。這樣我們可以用以下方式來表示循環\(a\)的執行次數。令

\[y_{ij} = \begin{cases} 1, if\ q(i)<q(k)\ for\ i < k \le j \\ 0, \ otherwise \end {cases} \]

那么

\[a=\sum_{1\le i<j \le n}y_{ij} \]

對於上面的例子而言第一行只有\(y_{12}=y_{13}=1\)，因為當外循環中的變量\(j=q(1)\)時，循環\(a\)會被執行\(2\)次，直到遍歷回到\(q(1)\)發現\(q(1)\)就是頭元素為止，由於我們構造\(q\)的時候讓頭元素遞減，所以如果考慮成因為發現\(q(1)<q(4)\)而退出，不會對研究循環\(a\)的執行次數造成影響。同樣地對於外循環變量\(j=q(2)\)的情況，只會往后執行\(1\)次就會發現\(q(2)<q(4)\)，因此第二行只有\(y_{23}=1\)。類似分析還能找到\(y_{45}=y_{78}=y_{79}=1\)。

為什么要定義這么一個\(y_{ij}\)呢？因為我們非常容易就可以發現\(y_{ij}\)的平均值\(\bar{y}_{ij}\)。這個平均值就是總共\(n!\)個排列中\(y_{ij}=1\)的排列個數。這也就是\(q(i)=min(q(k)|i \le k \le j)\)的概率，也就是\({1}/({j-i+1})\)（給K神跪了）因此

\[\begin{align} \bar{a} &= \sum_{1\le i<j\le n}\bar{y}_{ij}\\ &= \sum_{1\le i < j \le n}\frac{1}{j-i+1} \ \ {(我們固定i先對j求和，然后對i求和)}\\ &= (\frac{1}{2-1+1}+\frac{1}{3-1+1}+...+\frac{1}{n-1+1}) \\ &+(\frac{1}{3-2+1}+\frac{1}{4-2+1}+...+\frac{1}{n-2+1})+...+(\frac{1}{n-(n-1)+1}) \\ &= \frac{n-1}{2}+\frac{n-2}{3}+...+\frac{1}{n} \\ &= \sum_{2 \le r \le n}\frac{n+1-r}{r} \ \ (令r=j-i+1) \\ &= (n+1)(H_n-1)-(n-1) \\ &= (n+1)H_n - 2n \\ \end {align} \]

其中\(H_n\)是調和級數，使用積分法可以證明

\[H_n=\sum_{i=1}^{n}\frac{1}{i}=O(logn) \]

所以我們得到\(a\)循環的平均執行次數為\(O(logn)\)。

下面分析\(b\)的執行次數。顯然，環有多少個，\(b\)就會被進入多少次。這所有次數加起來的循環\(c\)的執行代價為\(O(n)\)固定不變。因此我們需要考慮\(b\)平均會被進入多少次，也就是所有長度為\(n\)排列中平均有多少個環。在Knuth的TAOCP, 1.2.10中已經對此有過分析，因此原文沒有寫出。（因為我沒有TAOCP有大概也看不懂）所以我們只寫一個結論：有\(k\)個環的排列共有\([n,k]\)個，這個\([n,k]\)是第一類Stirling數。最后可以得到\(b\)的平均值剛好是調和級數！（我推了推死活推不出來…）

\[\bar{b} = H_n = O(logn) \]

通過以上討論我們就嚴格證明了平均情況下這個算法的復雜度是\(O(n\log n)\)。

原文還有關於\(a,b\)的方差討論，由於計算極其繁雜在此就不寫了。重要的結論是，它的方差證明了\(O(n^2)\)的最壞情況是非常罕見的。

算法的改進空間

我們可以通過增加一個變量來避免全部元素移動完成后的不必要的遍歷，但這不會改善平均情況和最壞情況的復雜度。但是有一個方法可以將最壞情況復雜度降低到\(O(n\log n)\)，這個方法的關鍵是對於外循環遍歷到的一個\(j\)，我們同時搜索\(p(j), p^{-1}(j),p(p(j)),p^{-1}(p^{-1}(j)),…\),其中\(p^{-1}(j)\)是其反函數。這樣，我們重新考慮最壞情況，顯然就是整個排列只有一個長度為\(n\)的環，例如排列\((1,2,…,n)\)。這樣，最壞情況就是我們從\(j\)開始向環的兩邊搜索到頭元素的時間，再加上j在環前面和環后面的長度分別為\(k\)和\(n-k\)的部分都處在最壞情況的時間之和。設最壞情況為\(f(n)\)，就得到如下遞推式：

\[\begin{cases} f(n)=\max_{1\le k < n}\{\min\{k, n-k\} + f(k) + f(n-k)\} \\ f(1) = 0 \end{cases} \]

(內心OS：這種東西TM也能解？！)

K神拒絕回答你的提問並直接把答案甩在了你臉上：

\[f(n)=\sum_{0 \le n < n}v(k) \ \ (v(k)是k的二進制表示中1的個數) \\ f(2^{a_1}+2^{a_2}+...+2^{a_r})=\frac{1}{2}(a_12^{a_1}+(a_2+2)2^{a_2}+...+(a_r+2r-2)2^{a_r}) \\ 其中a_1>a_2>...>a_r \]

在其他幾篇由Bush、Mirsky、Drazin、和Griffith發表的論文里有對其復雜度的詳細分析。通過這些分析我們知道了這個解法在最壞情況下的復雜度為\(O(n\log n)\)。（我的內心已經崩潰了）

例二：Selecting the t-th largest

1. 問題定義

給定一個數組\(a_1,a_2,…,a_n\)，試用盡可能小的比較次數找出其中第\(t\)大的值。

2. 設計算法

這個問題相對比較常見一些，甚至曾經出現為數據結構作業題……但是其實這個問題想到算法容易，算法的分析並沒有那么容易。首先一個\(O(n\log n)\)的排序算法總能解決我們的問題，但有沒有復雜度更低的方法呢？通過快速排序思路的啟發我們可以想到這樣一個算法，假如我們對數組\(a_i,a_{i+1},…,a_j\)搜索，首先調用其中的Partition()方法將數組頭元素\(a_i\)的位置放到一個位置使得他左邊的元素都比他小，右邊的元素都比他大；然后根據他的位置\(k\)來縮小我們搜索第\(t\)大數的范圍：如果\(k=t\)我們就找到了這個元素；\(k<t\)則對\(a_{k+1},…,a_j\)搜索；\(k>t\)就對\(a_i,…,a_{k-1}\)搜索。這是一個標准的分治算法。可以寫出以下偽代碼：

FindtthNumber(a, i, j, t):
    key = a[i]
    # Partition的實現請參考快速排序相關資料
    # Partition返回的是分割后的數組下標
    # 減去數組開頭的位置得到a[k]是a[i]-a[j]里第幾大的數
    k = Partition(key, a, i, j) - i + 1 
    if k == t:
        return a[k]
    else if k < t:
        return FindtthNumber(a, k + 1, j, t - k)
    else:
        return FindtthNumber(a, i, k - 1, t)

3. 分析算法

通過偽代碼我們可以看出，影響一個子問題的變量有兩個：\(n\)和\(t\)，\(n\)是這個子問題中數組的長度，\(t\)是這個子問題中我們要找的第\(t\)大的數。不妨設這個子問題為\(C(n, t)\)。假設我們的分割到什么位置完全隨機，即這個子問題分割找到第\(1\)大、第\(2\)大、…、第\(n\)大的數的個數完全相等，均為\(1/n\)，我們就能分析得到以下遞推方程：

\[\begin{cases} C(1,1) = 0\\ C(n, t) = n - 1 + \frac{1}{n}(A(n,t)+B(n,t))\\ \end{cases}\\ \mbox{其中} \\ A(n, t) = C(n - 1,t-1)+C(n-2,t-2)+...+C(n-t+1,1)\\ B(n,t) = C(t,t)+C(t + 1,t) + ...+C(n-1, t) \]

顯然\(A(n,t)\)對應於遞歸函數中\(k<t\)的情況，\(B(n,t)\)對應於遞歸函數中\(k>t\)的情況。

之后我們任務就是求解這個遞推方程了。一般人看到這個遞推方程大概就會放棄了吧，然而Knuth發現這個遞推方程是可以求解的（跪了跪了）！依然是使用差消迭代法。首先通過觀察我們看到

\[A(n+1,t+1) = A(n,t)+C(n,t)\\ B(n+1,t) = B(n,t)+C(n,t) \]

這樣我們就可以把\(C(n+1,t+1), C(n,t+1), C(n,t),C(n-1,t)\)按如下方式相減：

\[(n+1)C(n+1,t+1)-nC(n,t+1)-nC(n,t)+(n-1)C(n-1,t) \\ = (n+1)n - n(n-1) - n(n-1) + (n-1)(n-2)\\ + A(n+1,t+1)-A(n,t+1)-A(n,t)+A(n-1,t) \\ +B(n+1,t+1)-B(n,t+1)-B(n,t)+B(n-1,t) \\ = 2 + C(n,t) - C(n-1,t)+C(n,t+1)-C(n-1,t) \]

推出來

\[C(n+1,t+1)-C(n,t+1)-C(n,t)+C(n-1,t)=\frac{2}{n+1}\\ (C(n+1,t+1)-C(n,t))-(C(n,t+1)-C(n-1,t))=\frac{2}{n+1} \]

之后需要列出邊界條件，類似之前的推導

\[\begin{align} C(n,1)&=n-1+\frac{1}{n}(C(1,1)+C(2,1)+...+C(n-1,1)) \\ (n+1)C(n+1,1)-nC(n,1)&=(n+1)n-n(n-1)+C(n,1)\\ C(n+1,1)-C(n,1)&=\frac{2n}{n+1}=2-\frac{2}{n+1}\\ C(n,1)&=2n-2H_n\\ \Rightarrow C(n,n)&=2n-2H_n\ \ (對稱性) \end{align} \]

由上述兩式可以計算得到

\[C(n+1,t+1)-C(n,t)=\frac{2}{n+1}+\frac{2}{n}+...+\frac{2}{t+2}+C(t+1,t+1)-C(t,t)\\ =2(H_{n+1}-H_{t+1})+2-\frac{2}{t+1} \]

不斷迭代

\[C(n,t)=2\sum_{k=2}^{t}(H_{n-t+k}-H_k+1-\frac{1}{k})+C(n+1-t,1)\\ =2((n+1)H_n-(n+3-t)H_{n+1-t}-(t+2)H_t+n+3) \]

由於調和級數\(H_n=O(\log n)\)，我們就嚴格證明了無論\(n,t\)取什么值，平均情況下算法的復雜度\(C(n,t)=O(n)\)。我們的算法比較次數已經足夠小了，那么如何尋找最少的比較次數？這就從一個A類問題變成了B類問題，並且非常難，感興趣的讀者可以參考相關論文了解對這方面的進展。

總結

被前面兩個算法搞暈了之后，相信大家已經了解算法分析是個什么樣的過程了（笑）。最后Knuth對算法分析做了如下總結：

1. 算法分析對計算機科學領域十分重要
2. 算法分析與離散數學密切相關
3. 算法分析正在形成科學方法
4. 算法分析領域還有很多問題沒有解決

參考文獻

[1] "Mathematical Analysis of Algorithms" Donald Knuth, 1971.