算法的完整實現代碼我已經上傳到了GitHub倉庫：NumericalAnalysis-Python（包括其它數值分析算法），感興趣的童鞋可以前往查看。
1 冪迭代算法（簡稱冪法）
1.1 占優特征值和占優特征向量
已知方陣\(\boldsymbol{A} \in \R^{n \times n}\), \(\boldsymbol{A}\)的占優特征值是比\(\boldsymbol{A}\)的其他特征值（的絕對值）都大的特征值\(\lambda\)，若這樣的特征值存在，則與\(\lambda\)相關的特征向量我們稱為占優特征向量。
1.2 占優特征值和占優特征向量的性質
如果一個向量反復與同一個矩陣相乘，那么該向量會被推向該矩陣的占優特征向量的方向。如下面這個例子所示：
import numpy as np
def prime_eigen(A, x, k):
x_t = x.copy()
for j in range(k):
x_t = A.dot(x_t)
return x_t
if __name__ == '__main__':
A = np.array(
[
[1, 3],
[2, 2]
]
)
x = np.array([-5, 5])
k = 4
r = prime_eigen(A, x, k)
print(r)
該算法運行結果如下：
250, 260
為什么會出現這種情況呢？因為對\(\forall \boldsymbol{x} \in \R^{n}\)都可以表示為\(A\)所有特征向量的線性組合(假設所有特征向量張成\(\R^n\)空間)。我們設\(\boldsymbol{x}^{(0)} = (-5, 5)^T\)，則冪迭代的過程可以如下表示：
可以看出是和占優特征值對應的特征向量在多次計算后會占優。在這里4是最大的特征值，而計算就越來越接近占優特征向量\((1, 1)^T\)的方向。
不過這樣重復與矩陣連乘和容易導致數值上溢，我們必須要在每步中對向量進行歸一化。就這樣，歸一化和與矩陣\(\boldsymbol{A}\)的連乘構成了如下所示的冪迭代算法：
import numpy as np
def powerit(A, x, k):
for j in range(k):
# 每次迭代前先對上一輪的x進行歸一化
u = x/np.linalg.norm(x)
# 計算本輪迭代未歸一化的x
x = A.dot(u)
# 計算出本輪對應的特征值
lam = u.dot(x)
# 最后一次迭代得到的特征向量x需要歸一化為u
u = x / np.linalg.norm(x)
return u, lam
if __name__ == '__main__':
A = np.array(
[
[1, 3],
[2, 2]
]
)
x = np.array([-5, 5])
k = 10
# 返回占優特征值和對應的特征向量
u, lam = powerit(A, x, k)
print("占優的特征向量:\n", u)
print("占優的特征值:\n", lam)
算法運行結果如下：
占優的特征向量:
[0.70710341 0.70711015]
占優的特征值:
3.9999809247674625
觀察上面的代碼，我們在第\(t\)輪迭代的第一行，得到的是歸一化后的第\(t-1\)輪迭代的特征向量近似值\(\boldsymbol{u}^{(t-1)}\)（想一想，為什么），然后按照\(\boldsymbol{x}^{(t)} = \boldsymbol{A}\boldsymbol{u}^{(t-1)}\)計算出第\(t\)輪迭代未歸一化的特征向量近似值\(\boldsymbol{x}^{(t)}\)，需要計算出第\(t\)輪迭代對應的特征值。這里我們我們直接直接運用了結論\(\lambda^{(t)} = (\boldsymbol{u}^{(t-1)})^T \boldsymbol{x^{(t)}}\)。該結論的推導如下：
證明
對於第\(t\)輪迭代，其中特征值的近似未\(\boldsymbol{\lambda}^{(t)}\)，我們想解特征方程
以得到第\(t\)輪迭代時該特征向量對應的特征值\(\lambda^{(t)}\)。我們采用最小二乘法求方程\((2)\)的近似解。我們可以寫出該最小二乘問題的正規方程為
然而我們可以寫出該最小二乘問題的解為
式子\((4)\)就是瑞利（Rayleigh）商。給定特征向量的近似，瑞利商式特征值的最優近似。又由歸一化的定義有
則我們可以將式\((4)\)寫作：
又因為前面已經計算出\(\boldsymbol{x}^{(t)} = \boldsymbol{A} \boldsymbol{u}^{(t-1)}\)，為了避免重復計算，我們將\(\boldsymbol{x}^{(t)}\)代入式\((5)\)得到：
證畢。
可以看出，冪迭代本質上每步進行歸一化的不動點迭代。
2 逆向冪迭代
上面我們的冪迭代算法用於求解（絕對值）最大的特征值。那么如何求最小的特征值呢？我們只需要將冪迭代用於矩陣的逆即可。
我們有結論，矩陣\(\boldsymbol{A}^{-1}\)的最大特征值就是矩陣\(\boldsymbol{A}\)的最小特征值的倒數。事實上，對矩陣\(\boldsymbol{A} \in \R^{n \times n}\) ,令其特征值表示為\(\lambda_1, \lambda_2, ..., \lambda_n\)，如果其逆矩陣存在，則逆矩陣\(A\)的特征值為\(\lambda_1^{-1}, \lambda_2^{-1}, ..., \lambda_n^{-1}\), 特征向量和矩陣\(A\)相同。該定理證明如下：
證明
有特征值和特征向量定義有
這蘊含着
因而
得證。
對逆矩陣\(\boldsymbol{A}^{-1}\)使用冪迭代，並對所得到的的\(\boldsymbol{A}^{-1}\)的特征值求倒數，就能得到矩陣\(\boldsymbol{A}\)的最小特征值。冪迭代式子如下：
但這要求我們對矩陣\(\boldsymbol{A}\)求逆，當矩陣\(\boldsymbol{A}\)過大時計算復雜度過高。於是我們需要稍作修改，對式\((11)\)的計算等價於
這樣，我們就可以采用高斯消元對\(\boldsymbol{x}^{(t+1)}\)進行求解，
不過，我們現在這個算法用於找出矩陣最大和最小的特征值，如何找出其他特征值呢？
如果我們要找出矩陣\(\boldsymbol{A}\)在實數\(s\)附近的特征值，可以對矩陣做出接近特征值的移動。我們有定理：對於矩陣\(\boldsymbol{A} \in \R^{n \times n}\)，設其特征值為\(\lambda_1, \lambda_2, ..., \lambda_n\)，則其轉移矩陣\(\boldsymbol{A}-sI\)的特征值為\(\lambda_1 -s, \lambda_2 -s, ..., \lambda_n -s\)，而特征向量和矩陣\(A\)相同。該定理證明如下：
證明
有特征值和特征向量定義有
從兩側減去\(sI\boldsymbol{v}\)，得到
因而矩陣\(\boldsymbol{A} - sI\)的特征值為\(\lambda - s\)，特征向量仍然為\(\boldsymbol{v}\)，得證。
這樣，我們想求矩陣\(\boldsymbol{A}\)在實數\(s\)附近的特征值，可以先對矩陣\((\boldsymbol{A}-sI)^{-1}\)使用冪迭代求出\((\boldsymbol{A}-sI)^{-1}\)的最大特征值\(b\)（因為我們知道轉移后的特征值為\((\lambda - s)^{-1}\)，要使\(\lambda\)盡可能接近\(s\)，就得取最大的特征值），其中每一步的\(x^{(t)}\)可以對\((\boldsymbol{A}-sI)\boldsymbol{x}^{(t)}=\boldsymbol{u}^{(t-1)}\)進行高斯消元得到。最后，我們計算出\(\lambda = b^{-1} + s\)即為矩陣\(A\)在\(s\)附近的特征值。該算法的實現如下：
import numpy as np
def powerit(A, x, s, k):
As = A-s*np.eye(A.shape[0])
for j in range(k):
# 為了讓數據不失去控制
# 每次迭代前先對x進行歸一化
u = x/np.linalg.norm(x)
# 求解(A-sI)xj = uj-1
x = np.linalg.solve(As, u)
lam = u.dot(x)
lam = 1/lam + s
# 最后一次迭代得到的特征向量x需要歸一化為u
u = x / np.linalg.norm(x)
return u, lam
if __name__ == '__main__':
A = np.array(
[
[1, 3],
[2, 2]
]
)
x = np.array([-5, 5])
k = 10
# 逆向冪迭代的平移值，可以通過平移值收斂到不同的特征值
s = 2
# 返回占優特征值和對應的特征值
u, lam = powerit(A, x, s, k)
# u為 [0.70710341 0.70711015]，指向特征向量[1, 1]的方向
print("占優的特征向量:\n", u)
print("占優的特征值:\n", lam)
算法運行結果如下：
占優的特征向量:
[0.64221793 0.7665221 ]
占優的特征值:
4.145795530352381
3 冪迭代的應用：PageRank算法
冪迭代的一大應用就是PageRank算法。PageRank算法作用在有向圖上的迭代算法，收斂后可以給每個節點賦一個表示重要性程度的值，該值越大表示節點在圖中顯得越重要。
比如，給定以下有向圖：
其鄰接矩陣為：
我們將鄰接矩陣轉置后，再沿着列歸一化，就得到了馬爾可夫概率轉移矩陣\(\boldsymbol{M}\)：
我們定義上網者從一個頁面轉移到另一個隨機頁面的概率是\(q\)，點擊本頁面上鏈接的概率是\(1-q\)。設圖的節點數為\(n\)，然后我們可以計算Google矩陣做為有向圖的一般轉移矩陣。對矩陣每個元素而言，我們有：
注意，Google矩陣每一列求和為1，這是一個隨機矩陣，它滿足一個性質，即占優特征值為1.
我們采用矩陣表示形式，即：
其中\(\boldsymbol{E}\)為元素全為1的\(n \times n\)方陣。
然后我們定義向量\(\boldsymbol{p}\)，其元素\(\boldsymbol{p}_i\)是待在頁面\(i\)上的概率。我們由前面的冪迭代算法知道，矩陣與向量重復相乘后向量會被推到特征值為1的方向。而這里，與特征值１對應的特征向量是一組頁面的穩態概率，根據定義這就是\(i\)個頁面的等級，即PageRank算法名字中的Rank的由來。（同時，這也是\(\boldsymbol{G}\)定義的馬爾科夫過程的穩態解)。故我們定義迭代過程：
注意，每輪迭代后我們要對\(\boldsymbol{p}\)向量歸一化（為了減少時間復雜度我們除以\(p\)向量所有維度元素中的最大值即可，以近似二范數歸一化）；而且，我們在所有輪次的迭代結束后也要對\(p\)向量進行歸一化（除以所有維度元素之和以保證所有維度之和為1）。
我們對該圖的PageRank算法代碼實現如下(其中移動到一個隨機頁面的概率\(q\)按慣例取0.15)：
import numpy as np
# 歸一化同時迭代，k是迭代步數
# 欲推往A特征值的方向，A肯定是方陣
def PageRank(A, p, k, q):
assert(A.shape[0]==A.shape[1])
n = A.shape[0]
M = A.T.astype(np.float32) #注意要轉為浮點型
for i in range(n):
M[:, i] = M[:, i]/np.sum(M[:, i])
G = (q/n)*np.ones((n,n)) + (1-q)*M
#G_T = G.T
p_t = p.copy()
for i in range(k):
y = G.dot(p_t)
p_t = y/np.max(y)
return p_t/np.sum(p_t)
if __name__ == '__main__':
A = np.array(
[
[0, 1, 1],
[0, 0, 1],
[1, 0, 0]
]
)
k = 20
p = np.array([1, 1, 1])
q = 0.15 #概率為1移動到一個隨機頁面通常為0.15
# 概率為1-q移動到與本頁面鏈接的頁面
R= PageRank(A, p, k, q)
print(R)
該算法運行結果如下：
[0.38779177 0.21480614 0.39740209]
可以看到20步迭代結束后網頁的Rank向量\(\boldsymbol{R}=(0.38779177, 0.21480614, 0.39740209)^T\)，這也可以看做網頁的重要性程度。
4 知名程序庫和源碼閱讀建議
PageRank算法有很多優秀的開源實現，這里推薦幾個項目：
4.1 Spark-GraphX
GraphX是一個優秀的分布式圖計算庫，從屬於Spark分布式計算框架，采用Scala語言實現了很多分布式的圖計算算法，也包括我們這里所講的PageRank算法。
文檔地址：https://spark.apache.org/graphx
源碼地址：https://github.com/apache/spark
4.2 neo4j
neo4j是一個采用Java實現的知名的圖數據庫，該數據庫也提供了PageRank算法的實現。
文檔地址：https://neo4j.com/
源碼地址：https://github.com/neo4j/neo4j.git
4.3 elastic-search
如果你有興趣深入研究搜索引擎的實現，那么向你推薦elastic-search項目，它是基於Java實現的一個搜索引擎。
文檔地址：https://www.elastic.co/cn/
源碼地址：https://github.com/elastic/elasticsearch.git
參考
- [1] Timothy sauer. 數值分析(第2版)[M].機械工業出版社, 2018.
- [2] 李航. 統計學習方法(第2版)[M]. 清華大學出版社, 2019.