幂法
含义
幂法是计算矩阵的按模最大的特征值和相应特征向量的一种向量迭代法
适用于大型稀疏矩阵
基础知识
求矩阵A特征值问题等价于求A的特征方程。
∣
λ
E
−
A
∣
=
[
λ
−
a
11
−
a
12
⋯
−
a
1
n
−
a
21
λ
−
a
22
⋯
−
a
2
n
⋯
⋯
⋯
⋯
λ
−
a
n
1
−
a
n
2
⋯
λ
−
a
n
n
]
=
0
|\lambda E-A| = \left[ \begin{matrix} \lambda-a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda-a_{22} & \cdots & -a_{2n} \\ \cdots & \cdots & \cdots & \cdots \\ \lambda-a_{n1} & -a_{n2} & \cdots & \lambda-a_{nn} \\ \end{matrix} \right] = 0
∣λE−A∣=⎣⎢⎢⎡λ−a11−a21⋯λ−an1−a12λ−a22⋯−an2⋯⋯⋯⋯−a1n−a2n⋯λ−ann⎦⎥⎥⎤=0
求A的属于特征值\lambda的特征向量等价求
(
λ
E
−
A
)
x
=
0
(\lambda E -A)x =0
(λE−A)x=0 的非零解
设
λ
为
A
∈
R
x
∗
n
的
特
征
值
,
x
称
为
A
的
与
特
征
值
λ
相
对
应
的
一
个
特
征
向
量
,
即
A
x
=
λ
x
,
(
x
≠
0
)
,
则
有
:
设\lambda 为A\in R^{x*n}的特征值,x称为A的与特征值\lambda 相对应的一个特征向量,即Ax = \lambda x,(x\neq0),则有:
设λ为A∈Rx∗n的特征值,x称为A的与特征值λ相对应的一个特征向量,即Ax=λx,(x̸=0),则有:
(1)
c
x
(
c
≠
0
)
也
是
A
的
与
特
征
值
λ
相
对
应
的
一
个
特
征
向
量
,
即
A
c
x
=
λ
c
x
cx(c\neq0)也是A的与特征值\lambda 相对应的一个特征向量,即Acx = \lambda cx
cx(c̸=0)也是A的与特征值λ相对应的一个特征向量,即Acx=λcx
(2)
λ
k
为
A
k
的
特
征
值
,
A
k
x
=
λ
k
x
\lambda ^k为A^k的特征值,A^kx=\lambda ^kx
λk为Ak的特征值,Akx=λkx
思想
设n阶矩阵A的特征值和特征向量为
λ
1
,
λ
2
,
.
.
.
,
λ
n
v
1
,
v
2
,
.
.
.
v
n
\lambda_1,\lambda_2,...,\lambda_n v_1,v_2,...v_n
λ1,λ2,...,λn v1,v2,...vn
v
1
,
v
2
,
.
.
.
v
n
线
性
无
关
,
且
v_1,v_2,...v_n线性无关,且
v1,v2,...vn线性无关,且
∣
λ
1
∣
>
∣
λ
2
∣
>
.
.
.
>
∣
λ
n
∣
|\lambda_1|>|\lambda_2|>...>|\lambda_n|
∣λ1∣>∣λ2∣>...>∣λn∣
称
λ
1
为
A
的
按
模
最
大
特
征
值
(
主
特
征
值
)
称\lambda_1为A的按模最大特征值(主特征值)
称λ1为A的按模最大特征值(主特征值)
任取非零向量
x
0
=
(
x
1
(
0
)
,
x
x
(
0
)
,
⋯
 
,
x
n
(
0
)
)
x_0 =(x_1^{(0)},x_x^{(0)},\cdots,x_n^{(0)})
x0=(x1(0),xx(0),⋯,xn(0))
特征向量构成n维空间的一组基.任取的初始向量X(0)由它们的线性组合给出
x
0
=
a
1
v
1
+
a
2
v
2
+
…
+
a
n
v
n
x_0 =a_1v_1+a_2v_2+…+a_nv_n
x0=a1v1+a2v2+…+anvn
假
设
a
1
≠
0
,
由
此
构
造
向
量
序
列
x
k
=
A
x
k
−
1
,
k
=
1
,
2
,
3...
假设a_1\neq0,由此构造向量序列x_k = Ax_{k-1},k=1,2,3...
假设a1̸=0,由此构造向量序列xk=Axk−1,k=1,2,3...
其中,
x
1
=
A
x
0
=
a
1
A
v
1
+
a
2
A
v
2
+
.
.
.
+
a
n
A
v
n
=
x_1=Ax_0=a_1Av_1+a_2Av_2+...+a_nAv_n=
x1=Ax0=a1Av1+a2Av2+...+anAvn=
=
a
1
λ
1
v
1
+
a
2
λ
2
v
2
+
.
.
.
+
a
n
λ
n
v
n
=a_1\lambda_1v_1+a_2\lambda_2v_2+...+a_n\lambda_nv_n
=a1λ1v1+a2λ2v2+...+anλnvn
x
2
=
A
x
1
=
a
1
λ
1
A
v
1
+
a
2
λ
2
A
v
2
+
.
.
.
+
a
n
λ
n
A
v
n
=
x_2=Ax_1=a_1\lambda_1Av_1+a_2\lambda_2Av_2+...+a_n\lambda_nAv_n=
x2=Ax1=a1λ1Av1+a2λ2Av2+...+anλnAvn=
=
a
1
λ
1
2
v
1
+
a
2
λ
2
2
v
2
+
.
.
.
+
a
n
λ
n
2
v
n
=a_1\lambda_1^2v_1+a_2\lambda_2^2v_2+...+a_n\lambda_n^2v_n
=a1λ12v1+a2λ22v2+...+anλn2vn
得到向量序列:
x
1
=
a
1
λ
1
v
1
+
a
2
λ
2
v
2
+
.
.
.
+
a
n
λ
n
v
n
x_1=a_1\lambda_1v_1+a_2\lambda_2v_2+...+a_n\lambda_nv_n
x1=a1λ1v1+a2λ2v2+...+anλnvn
x
2
=
a
1
λ
1
2
v
1
+
a
2
λ
2
2
v
2
+
.
.
.
+
a
n
λ
n
2
v
n
x_2= a_1\lambda_1^2v_1+a_2\lambda_2^2v_2+...+a_n\lambda_n^2v_n
x2=a1λ12v1+a2λ22v2+...+anλn2vn
…
\dots
…
x
k
=
a
1
λ
1
k
v
1
+
a
2
λ
2
k
v
2
+
.
.
.
+
a
n
λ
n
k
v
n
x_k= a_1\lambda_1^kv_1+a_2\lambda_2^kv_2+...+a_n\lambda_n^kv_n
xk=a1λ1kv1+a2λ2kv2+...+anλnkvn
下面按模最大特征值λ1是单根的情况讨论:
x
k
=
a
1
λ
1
k
v
1
+
a
2
λ
2
k
v
2
+
.
.
.
+
a
n
λ
n
k
v
n
=
λ
1
k
[
a
1
v
1
+
a
2
(
λ
2
λ
1
)
k
v
2
+
.
.
.
+
a
n
(
λ
n
λ
1
)
2
v
n
]
=
λ
1
k
(
a
1
v
1
+
ε
k
)
x_k= a_1\lambda_1^kv_1+a_2\lambda_2^kv_2+...+a_n\lambda_n^kv_n=\lambda_1^k[a_1v_1 +a_2(\frac{\lambda_2}{\lambda_1})^kv_2+...+a_n(\frac{\lambda_n}{\lambda_1})^2v_n]=\lambda_1^k(a_1v_1+\varepsilon_k)
xk=a1λ1kv1+a2λ2kv2+...+anλnkvn=λ1k[a1v1+a2(λ1λ2)kv2+...+an(λ1λn)2vn]=λ1k(a1v1+εk)
ε
k
=
a
2
(
λ
2
λ
1
)
k
v
2
+
.
.
.
+
a
n
(
λ
n
λ
1
)
2
v
n
\varepsilon_k=a_2(\frac{\lambda_2}{\lambda_1})^kv_2+...+a_n(\frac{\lambda_n}{\lambda_1})^2v_n
εk=a2(λ1λ2)kv2+...+an(λ1λn)2vn
由于
∣
λ
1
∣
>
∣
λ
2
∣
>
.
.
.
>
∣
λ
n
∣
故
∣
λ
j
λ
1
∣
<
1
(
j
=
2
,
3
,
.
.
.
n
)
,
故
:
|\lambda_1|>|\lambda_2|>...>|\lambda_n| 故|\frac{\lambda_j}{\lambda_1}|<1(j=2,3,...n),故:
∣λ1∣>∣λ2∣>...>∣λn∣ 故∣λ1λj∣<1(j=2,3,...n),故:
lim
k
→
+
∞
ε
k
=
0
则
有
lim
k
→
+
∞
1
λ
1
k
x
k
=
lim
k
→
+
∞
λ
1
k
(
a
1
v
1
+
ε
k
)
λ
1
k
=
a
1
v
1
\lim_{k\rightarrow+\infty}\varepsilon_k=0 则有\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^k}x_k=\lim_{k\rightarrow+\infty}\frac{\lambda_1^k(a_1v_1+\varepsilon_k)}{\lambda_1^k}=a_1v_1
k→+∞limεk=0 则有k→+∞limλ1k1xk=k→+∞limλ1kλ1k(a1v1+εk)=a1v1
lim
k
→
+
∞
1
λ
1
k
x
k
=
a
1
v
1
就
是
λ
1
对
应
的
特
征
向
量
\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^k}x_k=a_1v_1 就是\lambda_1对应的特征向量
k→+∞limλ1k1xk=a1v1就是λ1对应的特征向量
用
向
量
x
k
+
1
的
第
i
个
分
量
与
向
量
x
k
的
第
i
个
分
量
之
比
的
极
限
:
用向量x_{k+1}的第i个分量与向量x_k的第i个分量之比的极限:
用向量xk+1的第i个分量与向量xk的第i个分量之比的极限:
lim
k
→
+
∞
(
x
k
+
1
)
i
(
x
k
)
i
=
lim
k
→
+
∞
λ
1
k
+
1
(
a
1
v
1
+
ε
k
+
1
)
i
λ
1
k
(
a
1
v
1
+
ε
k
)
i
=
lim
k
→
+
∞
λ
1
[
a
1
(
v
1
)
i
+
(
ε
k
+
1
)
i
]
a
1
(
v
1
)
i
+
(
ε
k
)
i
=
λ
1
\lim_{k\rightarrow+\infty}\frac{(x_{k+1})_i}{(x_k)_i}=\lim_{k\rightarrow+\infty}\frac{\lambda_1^{k+1}(a_1v_1+\varepsilon_{k+1})_i}{\lambda_1^k(a_1v_1+\varepsilon_k)_i}=\lim_{k\rightarrow+\infty}\frac{\lambda_1[a_1(v_1)_i+(\varepsilon_{k+1})_i]}{a_1(v_1)_i+(\varepsilon_k)_i}=\lambda_1
k→+∞lim(xk)i(xk+1)i=k→+∞limλ1k(a1v1+εk)iλ1k+1(a1v1+εk+1)i=k→+∞lima1(v1)i+(εk)iλ1[a1(v1)i+(εk+1)i]=λ1
得到两个极限:
lim
k
→
+
∞
1
λ
1
k
x
k
=
a
1
v
1
lim
k
→
+
∞
(
x
k
+
1
)
i
(
x
k
)
i
=
λ
1
\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^{k}}x_k=a_1v_1 \lim_{k\rightarrow+\infty}\frac{(x_{k+1})_i}{(x_k)_i}=\lambda_1
k→+∞limλ1k1xk=a1v1 k→+∞lim(xk)i(xk+1)i=λ1
因此当k趋近无穷时,就能近似得到按模最大的特征值和对应的特征向量:
λ
1
≈
(
x
k
+
1
)
i
(
x
k
)
i
x
k
≈
λ
1
k
a
1
v
1
\lambda_1 \approx \frac{(x_{k+1})_i}{(x_k)_i} x_k\approx\lambda_1^k a_1v_1
λ1≈(xk)i(xk+1)i xk≈λ1ka1v1
因此得幂方法的迭代公式:
x
k
=
A
x
k
−
1
k
=
1
,
2
,
.
.
.
x_k = Ax_{k-1} k=1,2,...
xk=Axk−1 k=1,2,...
当k充分大时,有:
{
λ
1
≈
(
x
k
+
1
)
i
(
x
k
)
i
x
k
≈
λ
1
k
a
1
v
1
\begin{cases} \lambda_1 \approx \frac{(x_{k+1})_i}{(x_k)_i} \\ x_k\approx\lambda_1^k a_1v_1 \end{cases}
{λ1≈(xk)i(xk+1)ixk≈λ1ka1v1
存在的问题是:当 ∣ λ 1 ∣ > 1 时 , x k 中 不 为 0 的 分 量 会 随 着 k → ∞ 而 趋 于 ∞ , 计 算 机 计 算 时 会 造 成 溢 出 , 当 |\lambda_1|>1时,x_k中不为0的分量会随着k\rightarrow \infty而趋于 \infty,计算机计算时会造成溢出,当 ∣λ1∣>1时,xk中不为0的分量会随着k→∞而趋于∞,计算机计算时会造成溢出,当|\lambda_1|<1时,x_k中的分量会随着k\rightarrow \infty而趋于0$,
实际计算时,为了避免计算过程中出现绝对值过大或过小的数参加运算,通常在每步迭代时,将向量“归一化”即用的按模最大的分量 。
修改计算公式:
{
y
k
=
A
x
k
−
1
m
k
=
m
a
x
(
y
k
)
(
k
=
1
,
2
,
.
.
.
)
x
k
=
y
k
/
m
k
\begin{cases} y_k=Ax_{k-1}\\ m_k=max(y_k) & (k=1,2,...)\\ x_k=y_k/m_k \end{cases}
⎩⎪⎨⎪⎧yk=Axk−1mk=max(yk)xk=yk/mk(k=1,2,...)
上式中
m
k
是
向
量
y
k
中
绝
对
值
最
大
的
一
个
分
量
,
这
时
x
k
分
量
的
模
最
大
为
1
m_k是向量y_k中绝对值最大的一个分量,这时x_k分量的模最大为1
mk是向量yk中绝对值最大的一个分量,这时xk分量的模最大为1
当k充分大时,有:
{
λ
1
≈
m
k
v
1
≈
x
k
\begin{cases} \lambda_1 \approx m_k \\ v_1 \approx x_k \end{cases}
{λ1≈mkv1≈xk
实例
求矩阵
[
2
3
2
10
3
4
3
6
1
]
的
主
特
征
值
和
其
对
应
的
特
征
向
量
\left[ \begin{matrix} 2 &3 & 2 \\ 10 & 3 & 4 \\ 3 & 6 & 1 \end{matrix} \right] 的主特征值和其对应的特征向量
⎣⎡2103336241⎦⎤的主特征值和其对应的特征向量
根据迭代公式:
{
y
k
=
A
x
k
−
1
m
k
=
m
a
x
(
y
k
)
(
k
=
1
,
2
,
.
.
.
)
x
k
=
y
k
/
m
k
\begin{cases} y_k=Ax_{k-1}\\ m_k=max(y_k) & (k=1,2,...)\\ x_k=y_k/m_k \end{cases}
⎩⎪⎨⎪⎧yk=Axk−1mk=max(yk)xk=yk/mk(k=1,2,...)
取
x
0
=
(
0
,
0
,
1
)
T
,
则
有
:
x_0=(0,0,1)^T,则有:
x0=(0,0,1)T,则有:
y
1
=
(
2
,
4
,
1
)
T
,
m
1
=
4
,
x
1
=
1
m
1
y
1
=
(
0.5
,
1
,
0.25
)
T
y_1 = (2,4,1)^T,m_1 = 4,x_1=\frac{1}{m_1}y_1=(0.5,1,0.25)^T
y1=(2,4,1)T,m1=4,x1=m11y1=(0.5,1,0.25)T
直到k=8
Python实现
#-*- coding:utf-8 -*-
#@author:whs
#@time: 2019/3/2111:25
import numpy as np
def Solve(mat, max_itrs, min_delta):
"""
mat 表示矩阵
max_itrs 表示最大迭代次数
min_delta 表示停止迭代阈值
"""
itrs_num = 0
delta = float('inf')
N = np.shape(mat)[0]
# 所有分量都为1的列向量
x = np.ones(shape=(N, 1))
#x = np.array([[0],[0],[1]])
while itrs_num < max_itrs and delta > min_delta:
itrs_num += 1
y = np.dot(mat, x)
#print(y)
m = y.max()
#print("m={0}".format(m))
x = y / m
print("***********第{}次迭代*************".format(itrs_num))
print("y = ",y)
print("m={0}".format(m))
print("x^T为:",x)
print(
"===================================")
IOS = np.array([[2, 3, 2],
[10, 3, 4],
[3, 6, 1]], dtype=float)
Solve(IOS, 10, 1e-10)
输出:
***********第1次迭代*************
y = [[ 7.]
[17.]
[10.]]
m=17.0
x^T为: [[0.41176471]
[1. ]
[0.58823529]]
===================================
***********第2次迭代*************
y = [[5. ]
[9.47058824]
[7.82352941]]
m=9.470588235294118
x^T为: [[0.52795031]
[1. ]
[0.82608696]]
===================================
***********第3次迭代*************
y = [[ 5.70807453]
[11.58385093]
[ 8.40993789]]
m=11.58385093167702
x^T为: [[0.49276139]
[1. ]
[0.72600536]]
===================================
***********第4次迭代*************
y = [[ 5.43753351]
[10.83163539]
[ 8.20428954]]
m=10.831635388739945
x^T为: [[0.50200485]
[1. ]
[0.75743775]]
===================================
***********第5次迭代*************
y = [[ 5.5188852 ]
[11.04979951]
[ 8.2634523 ]]
m=11.049799514875502
x^T为: [[0.49945569]
[1. ]
[0.74783731]]
===================================
***********第6次迭代*************
y = [[ 5.49458599]
[10.98590609]
[ 8.24620437]]
m=10.985906091381928
x^T为: [[0.50014864]
[1. ]
[0.75061668]]
===================================
***********第7次迭代*************
y = [[ 5.50153064]
[11.00395312]
[ 8.2510626 ]]
m=11.003953118800313
x^T为: [[0.49995948]
[1. ]
[0.74982713]]
===================================
***********第8次迭代*************
y = [[ 5.49957322]
[10.99890329]
[ 8.24970556]]
m=10.998903290037243
x^T为: [[0.50001105]
[1. ]
[0.75004801]]
===================================
***********第9次迭代*************
y = [[ 5.50011813]
[11.00030258]
[ 8.25008117]]
m=11.000302582696586
x^T为: [[0.49999699]
[1. ]
[0.74998675]]
===================================
***********第10次迭代*************
y = [[ 5.49996747]
[10.99991685]
[ 8.24997771]]
m=10.999916852432536
x^T为: [[0.50000082]
[1. ]
[0.75000364]]
===================================