# 1 幂迭代算法（简称幂法）

## 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


\begin{aligned} & \boldsymbol{x}^{(1)} = \boldsymbol{A}\boldsymbol{x}^{(0)} = 4(1,1)^T - 2(-3, 2)^T\\ & \boldsymbol{x}^{(2)} = \boldsymbol{A}^2\boldsymbol{x}^{(0)} = 4^2(1, 1)^T + 2(-3, 2)^T\\ & ...\\ & \boldsymbol{x}^{(4)} = \boldsymbol{A}^4\boldsymbol{x}^{(0)} = 4^4(1, 1)^T + 2(-3, 2)^T = 256(1, 1)^T + 2(-3, 2)^T\\ \end{aligned} \tag{1}

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


$\boldsymbol{x^{(t-1)}} \cdot \lambda^{(t)} = \boldsymbol{A}\boldsymbol{x}^{(t-1)} \tag{2}$

$(\boldsymbol{x}^{(t-1)})^T\boldsymbol{x}^{(t-1)} \cdot \lambda^{(t-1)} = (\boldsymbol{x}^{(t-1)})^T (\boldsymbol{A}\boldsymbol{x}^{(t-1)}) \tag{3}$

$\lambda^{(t)} = \frac{(\boldsymbol{x}^{(t-1)})^T\boldsymbol{A}\boldsymbol{x}^{(t-1)}}{(\boldsymbol{x}^{(t-1)})^T\boldsymbol{x}^{(t-1)}} \tag{4}$

$\boldsymbol{u}^{(t-1)} = \frac{\boldsymbol{x}^{(t-1)}}{||\boldsymbol{x}^{(t-1)}||} = \frac{\boldsymbol{x}^{(t-1)}}{{[(\boldsymbol{x}^{(t-1)})^T\boldsymbol{x}^{(t-1)}]}^{\frac{1}{2}}} \tag{5}$

$\lambda^{(t)} = (\boldsymbol{u}^{(t-1)})^T\boldsymbol{A}\boldsymbol{u}^{(t-1)} \tag{6}$

$\lambda^{(t)} = (\boldsymbol{u}^{(t-1)})^T\boldsymbol{x}^{(t)} \tag{7}$

# 2 逆向幂迭代

$\boldsymbol{A}\boldsymbol{v} = \lambda \boldsymbol{v} \tag{8}$

$\boldsymbol{v} = \lambda \boldsymbol{A}^{-1}\boldsymbol{v} \tag{9}$

$\boldsymbol{A}^{-1}\boldsymbol{v} = \lambda^{-1}\boldsymbol{v} \tag{10}$

$\boldsymbol{x}^{(t+1)} = \boldsymbol{A}^{-1}\boldsymbol{x}^{(t)} \tag{11}$

$\boldsymbol{A}\boldsymbol{x}^{(t+1)} = \boldsymbol{x}^{(t)} \tag{12}$

$\boldsymbol{A}\boldsymbol{v} = \lambda \boldsymbol{v} \tag{13}$

$(\boldsymbol{A} - sI)\boldsymbol{v} = (\lambda - s)\boldsymbol{v} \tag{14}$

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算法

$\left( \begin{matrix} 0 & 1 & 1 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ \end{matrix} \right) \tag{15}$

$\left( \begin{matrix} 0 & 0 & 1 \\ \frac{1}{2} & 0 & 0 \\ \frac{1}{2} & 1 & 0 \\ \end{matrix} \tag{16} \right)$

$\boldsymbol{G}_{ij} = \frac{q}{n} + (1-q)\boldsymbol{M}_{ij} \tag{17}$

$\boldsymbol{G} = \frac{q}{n}\boldsymbol{E} + (1-q)\boldsymbol{M} \tag{18}$

$\boldsymbol{p}_{t+1} = \boldsymbol{G}\boldsymbol{p}_{t} \tag{19}$

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]


# 4 知名程序库和源码阅读建议

PageRank算法有很多优秀的开源实现，这里推荐几个项目：

## 4.1 Spark-GraphX

GraphX是一个优秀的分布式图计算库，从属于Spark分布式计算框架，采用Scala语言实现了很多分布式的图计算算法，也包括我们这里所讲的PageRank算法。

## 4.2 neo4j

neo4j是一个采用Java实现的知名的图数据库，该数据库也提供了PageRank算法的实现。

# 参考

• [1] Timothy sauer. 数值分析(第2版)[M].机械工业出版社, 2018.
• [2] 李航. 统计学习方法(第2版)[M]. 清华大学出版社, 2019.
posted @ 2021-10-14 11:28  orion-orion  阅读(1072)  评论(2编辑  收藏  举报