數值分析：幂疊代和PageRank算法(Numpy實作)

計算矩陣特征值是一個經典數值計算問題，而很多相關方法都是基于幂疊代的基。如該思想的一個複雜版本被稱為QR算法，可以計算出典型矩陣所有的特征值。幂疊代的動機是與矩陣相乘可以将向量推向主特征值向量的方向，其在工業界的一大應用就是PageRank算法。PageRank算法作用在有向圖上的疊代算法，收斂後可以給每個節點賦一個表示重要性程度的值，該值越大表示節點在圖中顯得越重要。

已知方陣\(\bm{A} \in \R^{n \times n}\), \(\bm{A}\)的占優特征值是比\(\bm{A}\)的其他特征值（的絕對值）都大的特征值\(\lambda\)，若這樣的特征值存在，則與\(\lambda\)相關的特征向量我們稱為占優特征向量。

如果一個向量反複與同一個矩陣相乘，那麼該向量會被推向該矩陣的占優特征向量的方向。如下面這個例子所示：

該算法運作結果如下：

為什麼會出現這種情況呢？因為對\(\forall \bm{x} \in \R^{n}\)都可以表示為\(A\)所有特征向量的線性組合(假設所有特征向量張成\(\R^n\)空間)。我們設\(\bm{x}^{(0)} = (-5, 5)^T\)，則幂疊代的過程可以如下表示：

\[\begin{aligned}

& \bm{x}^{(1)} = \bm{A}\bm{x}^{(0)} = 4(1,1)^T - 2(-3, 2)^T\\

& \bm{x}^{(2)} = \bm{A}^2\bm{x}^{(0)} = 4^2(1, 1)^T + 2(-3, 2)^T\\

& ...\\

& \bm{x}^{(4)} = \bm{A}^4\bm{x}^{(0)} = 4^4(1, 1)^T + 2(-3, 2)^T = 256(1, 1)^T + 2(-3, 2)^T\\

\end{aligned} \tag{1}

\]

可以看出是和占優特征值對應的特征向量在多次計算後會占優。在這裡4是最大的特征值，而計算就越來越接近占優特征向量\((1, 1)^T\)的方向。

不過這樣重複與矩陣連乘和容易導緻數值上溢，我們必須要在每步中對向量進行歸一化。就這樣，歸一化和與矩陣\(\bm{A}\)的連乘構成了如下所示的幂疊代算法：

算法運作結果如下：

觀察上面的代碼，我們在第\(t\)輪疊代的第一行，得到的是歸一化後的第\(t-1\)輪疊代的特征向量近似值\(\bm{u}^{(t-1)}\)（想一想，為什麼），然後按照\(\bm{x}^{(t)} = \bm{A}\bm{u}^{(t-1)}\)計算出第\(t\)輪疊代未歸一化的特征向量近似值\(\bm{x}^{(t)}\)，需要計算出第\(t\)輪疊代對應的特征值。這裡我們我們直接直接運用了結論\(\lambda^{(t)} = (\bm{u}^{(t-1)})^T \bm{x^{(t)}}\)。該結論的推導如下：

證明

對于第\(t\)輪疊代，其中特征值的近似未\(\bm{\lambda}^{(t)}\)，我們想解特征方程

\[\bm{x^{(t-1)}} \cdot \lambda^{(t)} = \bm{A}\bm{x}^{(t-1)}

\tag{2}

以得到第\(t\)輪疊代時該特征向量對應的特征值\(\lambda^{(t)}\)。我們采用最小二乘法求方程\((2)\)的近似解。我們可以寫出該最小二乘問題的正規方程為

\[(\bm{x}^{(t-1)})^T\bm{x}^{(t-1)} \cdot \lambda^{(t-1)} = (\bm{x}^{(t-1)})^T (\bm{A}\bm{x}^{(t-1)})

\tag{3}

然而我們可以寫出該最小二乘問題的解為

\[\lambda^{(t)} = \frac{(\bm{x}^{(t-1)})^T\bm{A}\bm{x}^{(t-1)}}{(\bm{x}^{(t-1)})^T\bm{x}^{(t-1)}}

\tag{4}

式子\((4)\)就是瑞利（Rayleigh）商。給定特征向量的近似，瑞利商式特征值的最優近似。又由歸一化的定義有

\[\bm{u}^{(t-1)} = \frac{\bm{x}^{(t-1)}}{||\bm{x}^{(t-1)}||} = \frac{\bm{x}^{(t-1)}}{{[(\bm{x}^{(t-1)})^T\bm{x}^{(t-1)}]}^{\frac{1}{2}}}

\tag{5}

則我們可以将式\((4)\)寫作：

\[\lambda^{(t)} = (\bm{u}^{(t-1)})^T\bm{A}\bm{u}^{(t-1)}

\tag{6}

又因為前面已經計算出\(\bm{x}^{(t)} = \bm{A} \bm{u}^{(t-1)}\)，為了避免重複計算，我們将\(\bm{x}^{(t)}\)代入式\((5)\)得到：

\[\lambda^{(t)} = (\bm{u}^{(t-1)})^T\bm{x}^{(t)}

\tag{7}

證畢。

可以看出，幂疊代本質上每步進行歸一化的不動點疊代。

上面我們的幂疊代算法用于求解（絕對值）最大的特征值。那麼如何求最小的特征值呢？我們隻需要将幂疊代用于矩陣的逆即可。

我們有結論，矩陣\(\bm{A}^{-1}\)的最大特征值就是矩陣\(\bm{A}\)的最小特征值的倒數。事實上，對矩陣\(\bm{A} \in \R^{n \times n}\) ,令其特征值表示為\(\lambda_1, \lambda_2, ..., \lambda_n\)，如果其逆矩陣存在，則逆矩陣\(A\)的特征值為\(\lambda_1^{-1}, \lambda_2^{-1}, ..., \lambda_n^{-1}\), 特征向量和矩陣\(A\)相同。該定理證明如下：

有特征值和特征向量定義有

\[\bm{A}\bm{v} = \lambda \bm{v}

\tag{8}

這蘊含着

\[\bm{v} = \lambda \bm{A}^{-1}\bm{v}

\tag{9}

因而

\[\bm{A}^{-1}\bm{v} = \lambda^{-1}\bm{v}

\tag{10}

得證。

對逆矩陣\(\bm{A}^{-1}\)使用幂疊代，并對所得到的的\(\bm{A}^{-1}\)的特征值求倒數，就能得到矩陣\(\bm{A}\)的最小特征值。幂疊代式子如下：

\[\bm{x}^{(t+1)} = \bm{A}^{-1}\bm{x}^{(t)}

\tag{11}

但這要求我們對矩陣\(\bm{A}\)求逆，當矩陣\(\bm{A}\)過大時計算複雜度過高。于是我們需要稍作修改，對式\((11)\)的計算等價于

\[\bm{A}\bm{x}^{(t+1)} = \bm{x}^{(t)}

\tag{12}

這樣，我們就可以采用高斯消元對\(\bm{x}^{(t+1)}\)進行求解，

不過，我們現在這個算法用于找出矩陣最大和最小的特征值，如何找出其他特征值呢？

如果我們要找出矩陣\(\bm{A}\)在實數\(s\)附近的特征值，可以對矩陣做出接近特征值的移動。我們有定理：對于矩陣\(\bm{A} \in \R^{n \times n}\)，設其特征值為\(\lambda_1, \lambda_2, ..., \lambda_n\)，則其轉移矩陣\(\bm{A}-sI\)的特征值為\(\lambda_1 -s, \lambda_2 -s, ..., \lambda_n -s\)，而特征向量和矩陣\(A\)相同。該定理證明如下：

\tag{13}

從兩側減去\(sI\bm{v}\)，得到

\[(\bm{A} - sI)\bm{v} = (\lambda - s)\bm{v}

\tag{14}

因而矩陣\(\bm{A} - sI\)的特征值為\(\lambda - s\)，特征向量仍然為\(\bm{v}\)，得證。

這樣，我們想求矩陣\(\bm{A}\)在實數\(s\)附近的特征值，可以先對矩陣\((\bm{A}-sI)^{-1}\)使用幂疊代求出\((\bm{A}-sI)^{-1}\)的最大特征值\(b\)（因為我們知道轉移後的特征值為\((\lambda - s)^{-1}\)，要使\(\lambda\)盡可能接近\(s\)，就得取最大的特征值），其中每一步的\(x^{(t)}\)可以對\((\bm{A}-sI)\bm{x}^{(t)}=\bm{u}^{(t-1)}\)進行高斯消元得到。最後，我們計算出\(\lambda = b^{-1} + s\)即為矩陣\(A\)在\(s\)附近的特征值。該算法的實作如下：

幂疊代的一大應用就是PageRank算法。PageRank算法作用在有向圖上的疊代算法，收斂後可以給每個節點賦一個表示重要性程度的值，該值越大表示節點在圖中顯得越重要。

比如，給定以下有向圖：

其鄰接矩陣為：

\[\left(

\begin{matrix}

0 & 1 & 1 \\

0 & 0 & 1 \\

1 & 0 & 0 \\

\end{matrix}

\right)

\tag{15}

我們将鄰接矩陣沿着行歸一化，就得到了馬爾可夫機率轉移矩陣\(\bm{M}\)：

0 & \frac{1}{2} & \frac{1}{2} \\

\tag{16}

我們定義上網者從一個頁面轉移到另一個随機頁面的機率是\(q\)，停留在本頁面的機率是\(1-q\)。設圖的節點數為\(n\)，然後我們可以計算Google矩陣做為有向圖的一般轉移矩陣。對矩陣每個元素而言，我們有：

\[\bm{G}_{ij} = \frac{q}{n} + (1-q)\bm{M}_{ij}

\tag{17}

注意，Google矩陣每一列求和為1，這是一個随機矩陣，它滿足一個性質，即占優特征值為1.

我們采用矩陣表示形式，即：

\[\bm{G}_{ij} = \frac{q}{n}\bm{E} + (1-q)\bm{M}_{ij}

\tag{18}

其中\(\bm{E}\)為元素全為1的\(n \times n\)方陣。

然後我們定義向量\(\bm{p}\)，其元素\(\bm{p}_i\)是待在頁面\(i\)上的機率。我們由前面的幂疊代算法知道，矩陣與向量重複相乘後向量會被推到特征值為1的方向。而這裡，與特征值１對應的特征向量是一組頁面的穩态機率，根據定義這就是\(i\)個頁面的等級，即PageRank算法名字中的Rank的由來。（同時，這也是\(G^T\)定義的馬爾科夫過程的穩态解)。故我們定義疊代過程：

\[\bm{p}_{t+1} = \bm{G}^T\bm{p}_{t}

\tag{19}

注意，每輪疊代後我們要對\(\bm{p}\)向量歸一化（為了減少時間複雜度我們除以\(p\)向量所有次元元素中的最大值即可，以近似二範數歸一化）；而且，我們在所有輪次的疊代結束後也要對\(p\)向量進行歸一化（除以所有次元元素之和以保證所有次元之和為1）。

我們對該圖的PageRank算法代碼實作如下(其中移動到一個随機頁面的機率\(q\)按慣例取0.15)：

可以看到20步疊代結束後網頁的Rank向量\(\bm{R}=(0.38779177, 0.21480614, 0.39740209)^T\)，這也可以看做網頁的重要性程度。

PageRank算法有很多優秀的開源實作，這裡推薦幾個項目：

GraphX是一個優秀的分布式圖計算庫，從屬于Spark分布式計算架構，采用Scala語言實作了很多分布式的圖計算算法，也包括我們這裡所講的PageRank算法。

文檔位址：https://spark.apache.org/graphx

源碼位址：https://github.com/apache/spark

neo4j是一個采用Java實作的知名的圖資料庫，該資料庫也提供了PageRank算法的實作。

文檔位址：https://neo4j.com/

源碼位址：https://github.com/neo4j/neo4j.git

如果你有興趣深入研究搜尋引擎的實作，那麼向你推薦elastic-search項目，它是基于Java實作的一個搜尋引擎。

文檔位址：https://www.elastic.co/cn/

源碼位址：https://github.com/elastic/elasticsearch.git

[1] Timothy sauer. 數值分析(第2版)[M].機械工業出版社, 2018.

[2] 李航. 統計學習方法(第2版)[M]. 清華大學出版社, 2019.

數學是符号的藝術，音樂是上界的語言。