數值分析：冪迭代和PageRank算法

• 2021 年 10 月 14 日
• 筆記
• 數值分析, 矩陣計算, 算法

1. 冪迭代算法（簡稱冪法）

(1) 佔優特徵值和佔優特徵向量

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

(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 \bm{x} \in \R^{n}\)都可以表示為\(A\)所有特徵向量的線性組合(假設所有特徵向量張成\(\R^n\)空間)。我們設\(\bm{x}^{(0)} = (-5, 5)^T\)，則冪迭代的過程可以如下表示：

可以看出是和佔優特徵值對應的特徵向量在多次計算後會佔優。在這裡4是最大的特徵值，而計算就越來越接近特徵向量\((1, 1)^T\)的方向。
不過這樣重複與矩陣連乘和容易導致數值上溢，我們必須要在每步中對向量進行歸一化。就這樣，歸一化和與矩陣\(\bm{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\)輪迭代的特徵向量近似值\(\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)}}\)。該結論的推導如下：

證明

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

2. 逆向冪迭代

上面我們的冪迭代算法用於求解（絕對值）最大的特徵值。那麼如何求最小的特徵值呢？我們只需要將冪迭代用於矩陣的逆即可。

我們有結論，矩陣\(\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}^{-1}\)使用冪迭代，並對所得到的的\(\bm{A}^{-1}\)的特徵值求倒數，就能得到矩陣\(\bm{A}\)的最小特徵值。冪迭代式子如下：

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

這樣，我們就可以採用高斯消元對\(\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\)相同。該定理證明如下：

證明

這樣，我們想求矩陣\(\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\)附近的特徵值。該算法的實現如下：

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算法作用在有向圖上的迭代算法，收斂後可以給每個節點賦一個表示重要性程度的值，該值越大表示節點在圖中顯得越重要。
比如，給定以下有向圖：

其鄰接矩陣為：

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

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

注意，Google矩陣每一列求和為1，這是一個隨機矩陣，它滿足一個性質，即最大特徵為1.
我們採用矩陣表示形式，即：

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

注意，每輪迭代後我們要對\(\bm{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.copy().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_T.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 #移動到一個隨機頁面概率通常為0.15
    # 概率為1-q移動到與本頁面鏈接的頁面
    R= PageRank(A, p, k, q)
    print(R)
# [0.38779177 0.21480614 0.39740209]

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

知名程序庫和源碼閱讀建議

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

(1) Spark-GraphX

GraphX是一個優秀的分佈式圖計算庫，從屬於Spark分佈式計算框架，採用Scala語言實現了很多分佈式的圖計算算法，也包括我們這裡所講的PageRank算法。
文檔地址：//spark.apache.org/graphx
源碼地址：//github.com/apache/spark

(2) neo4j

neo4j是一個採用Java實現的知名的圖數據庫，該數據庫也提供了PageRank算法的實現。
文檔地址：//neo4j.com/
源碼地址：//github.com/neo4j/neo4j.git

(3) elastic-search

如果你有興趣深入研究搜索引擎的實現，那麼向你推薦elastic-search項目，它是基於Java實現的一個搜索引擎。
文檔地址：//www.elastic.co/cn/
源碼地址：//github.com/elastic/elasticsearch.git

參考文獻

• [1] Timothy sauer. 數值分析(第2版)[M].機械工業出版社, 2018.
• [2] 李航. 統計學習方法(第2版)[M]. 清華大學出版社, 2019.