数值分析:幂迭代和PageRank算法(Numpy实现)

算法的完整实现代码我已经上传到了GitHub仓库:NumericalAnalysis-Python(包括其它数值分析算法),感兴趣的童鞋可以前往查看。

1 幂迭代算法(简称幂法)

1.1 占优特征值和占优特征向量

已知方阵ARn×n, A的占优特征值是比A的其他特征值(的绝对值)都大的特征值λ,若这样的特征值存在,则与λ相关的特征向量我们称为占优特征向量。

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

为什么会出现这种情况呢?因为对xRn都可以表示为A所有特征向量的线性组合(假设所有特征向量张成Rn空间)。我们设x(0)=(5,5)T,则幂迭代的过程可以如下表示:

(1)x(1)=Ax(0)=4(1,1)T2(3,2)Tx(2)=A2x(0)=42(1,1)T+2(3,2)T...x(4)=A4x(0)=44(1,1)T+2(3,2)T=256(1,1)T+2(3,2)T

可以看出是和占优特征值对应的特征向量在多次计算后会占优。在这里4是最大的特征值,而计算就越来越接近占优特征向量(1,1)T的方向。
不过这样重复与矩阵连乘和容易导致数值上溢,我们必须要在每步中对向量进行归一化。就这样,归一化和与矩阵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轮迭代的第一行,得到的是归一化后的第t1轮迭代的特征向量近似值u(t1)(想一想,为什么),然后按照x(t)=Au(t1)计算出第t轮迭代未归一化的特征向量近似值x(t),需要计算出第t轮迭代对应的特征值。这里我们我们直接直接运用了结论λ(t)=(u(t1))Tx(t)。该结论的推导如下:

证明


对于第t轮迭代,其中特征值的近似未λ(t),我们想解特征方程

(2)x(t1)λ(t)=Ax(t1)

以得到第t轮迭代时该特征向量对应的特征值λ(t)。我们采用最小二乘法求方程(2)的近似解。我们可以写出该最小二乘问题的正规方程为

(3)(x(t1))Tx(t1)λ(t1)=(x(t1))T(Ax(t1))

然而我们可以写出该最小二乘问题的解为

(4)λ(t)=(x(t1))TAx(t1)(x(t1))Tx(t1)

式子(4)就是瑞利(Rayleigh)商。给定特征向量的近似,瑞利商式特征值的最优近似。又由归一化的定义有

(5)u(t1)=x(t1)||x(t1)||=x(t1)[(x(t1))Tx(t1)]12

则我们可以将式(4)写作:

(6)λ(t)=(u(t1))TAu(t1)

又因为前面已经计算出x(t)=Au(t1),为了避免重复计算,我们将x(t)代入式(5)得到:

(7)λ(t)=(u(t1))Tx(t)

证毕。


可以看出,幂迭代本质上每步进行归一化的不动点迭代。

2 逆向幂迭代

上面我们的幂迭代算法用于求解(绝对值)最大的特征值。那么如何求最小的特征值呢?我们只需要将幂迭代用于矩阵的逆即可。

我们有结论,矩阵A1的最大特征值就是矩阵A的最小特征值的倒数。事实上,对矩阵ARn×n ,令其特征值表示为λ1,λ2,...,λn,如果其逆矩阵存在,则逆矩阵A的特征值为λ11,λ21,...,λn1, 特征向量和矩阵A相同。该定理证明如下:

证明


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

(8)Av=λv

这蕴含着

(9)v=λA1v

因而

(10)A1v=λ1v

得证。


对逆矩阵A1使用幂迭代,并对所得到的的A1的特征值求倒数,就能得到矩阵A的最小特征值。幂迭代式子如下:

(11)x(t+1)=A1x(t)

但这要求我们对矩阵A求逆,当矩阵A过大时计算复杂度过高。于是我们需要稍作修改,对式(11)的计算等价于

(12)Ax(t+1)=x(t)

这样,我们就可以采用高斯消元对x(t+1)进行求解,
不过,我们现在这个算法用于找出矩阵最大和最小的特征值,如何找出其他特征值呢?
如果我们要找出矩阵A在实数s附近的特征值,可以对矩阵做出接近特征值的移动。我们有定理:对于矩阵ARn×n,设其特征值为λ1,λ2,...,λn,则其转移矩阵AsI的特征值为λ1s,λ2s,...,λns,而特征向量和矩阵A相同。该定理证明如下:

证明


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

(13)Av=λv

从两侧减去sIv,得到

(14)(AsI)v=(λs)v

因而矩阵AsI的特征值为λs,特征向量仍然为v,得证。


这样,我们想求矩阵A在实数s附近的特征值,可以先对矩阵(AsI)1使用幂迭代求出(AsI)1的最大特征值b(因为我们知道转移后的特征值为(λs)1,要使λ尽可能接近s,就得取最大的特征值),其中每一步的x(t)可以对(AsI)x(t)=u(t1)进行高斯消元得到。最后,我们计算出λ=b1+s即为矩阵As附近的特征值。该算法的实现如下:

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算法作用在有向图上的迭代算法,收敛后可以给每个节点赋一个表示重要性程度的值,该值越大表示节点在图中显得越重要。
比如,给定以下有向图:

其邻接矩阵为:

(15)(011001100)

我们将邻接矩阵转置后,再沿着列归一化,就得到了马尔可夫概率转移矩阵M

(16)(00112001210)

我们定义上网者从一个页面转移到另一个随机页面的概率是q,点击本页面上链接的概率是1q。设图的节点数为n,然后我们可以计算Google矩阵做为有向图的一般转移矩阵。对矩阵每个元素而言,我们有:

(17)Gij=qn+(1q)Mij

注意,Google矩阵每一列求和为1,这是一个随机矩阵,它满足一个性质,即占优特征值为1.
我们采用矩阵表示形式,即:

(18)G=qnE+(1q)M

其中E为元素全为1的n×n方阵。
然后我们定义向量p,其元素pi是待在页面i上的概率。我们由前面的幂迭代算法知道,矩阵与向量重复相乘后向量会被推到特征值为1的方向。而这里,与特征值1对应的特征向量是一组页面的稳态概率,根据定义这就是i个页面的等级,即PageRank算法名字中的Rank的由来。(同时,这也是G定义的马尔科夫过程的稳态解)。故我们定义迭代过程:

(19)pt+1=Gpt

注意,每轮迭代后我们要对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向量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

如果你有兴趣深入研究搜索引擎的实现,那么向你推荐elastic-search项目,它是基于Java实现的一个搜索引擎。
文档地址https://www.elastic.co/cn/
源码地址https://github.com/elastic/elasticsearch.git

参考

  • [1] Timothy sauer. 数值分析(第2版)[M].机械工业出版社, 2018.
  • [2] 李航. 统计学习方法(第2版)[M]. 清华大学出版社, 2019.
posted @   orion-orion  阅读(1662)  评论(2编辑  收藏  举报
编辑推荐:
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
点击右上角即可分享
微信分享提示