数值分析：幂迭代和PageRank算法

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\)，则幂迭代的过程可以如下表示：

\[\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}\)的连乘构成了如下所示的幂迭代算法：

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\)轮迭代的特征