幂法

含义

幂法是计算矩阵的按模最大的特征值和相应特征向量的一种向量迭代法

适用于大型稀疏矩阵

基础知识

求矩阵A特征值问题等价于求A的特征方程。

$|\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的特征向量等价求
$(\lambda E -A)x =0$(λE−A)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）,则有：
（１）$cx（c\neq0）也是A的与特征值\lambda 相对应的一个特征向量,即Acx = \lambda cx$cx（c̸​=0）也是A的与特征值λ相对应的一个特征向量,即Acx=λcx
（２）$\lambda ^k为A^k的特征值，A^kx=\lambda ^kx$λk为Ak的特征值，Akx=λkx

思想

设n阶矩阵A的特征值和特征向量为
$\lambda_1,\lambda_2,...,\lambda_n　　　 v_1,v_2,...v_n$λ1​,λ2​,...,λn​　　　v1​,v2​,...vn​
$v_1,v_2,...v_n线性无关，且$v1​,v2​,...vn​线性无关，且
$|\lambda_1|>|\lambda_2|>...>|\lambda_n|$∣λ1​∣>∣λ2​∣>...>∣λn​∣
$称\lambda_1为A的按模最大特征值（主特征值）$称λ1​为A的按模最大特征值（主特征值）

任取非零向量$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_1v_1+a_2v_2+…+a_nv_n$x0​=a1​v1​+a2​v2​+…+an​vn​
$假设a_1\neq0,由此构造向量序列x_k = Ax_{k-1}，k=1,2,3...$假设a1​̸​=0,由此构造向量序列xk​=Axk−1​，k=1,2,3...
其中，$x_1=Ax_0=a_1Av_1+a_2Av_2+...+a_nAv_n=$x1​=Ax0​=a1​Av1​+a2​Av2​+...+an​Avn​=
$＝a_1\lambda_1v_1+a_2\lambda_2v_2+...+a_n\lambda_nv_n$　　　　　　　　　　　　　　　　　　　　　＝a1​λ1​v1​+a2​λ2​v2​+...+an​λn​vn​

$x_2=Ax_1=a_1\lambda_1Av_1+a_2\lambda_2Av_2+...+a_n\lambda_nAv_n=$x2​=Ax1​=a1​λ1​Av1​+a2​λ2​Av2​+...+an​λn​Avn​=
$＝a_1\lambda_1^2v_1+a_2\lambda_2^2v_2+...+a_n\lambda_n^2v_n$　　　　　　　　　　　　　　　　　　　　　＝a1​λ12​v1​+a2​λ22​v2​+...+an​λn2​vn​

得到向量序列：
$x_1＝a_1\lambda_1v_1+a_2\lambda_2v_2+...+a_n\lambda_nv_n$x1​＝a1​λ1​v1​+a2​λ2​v2​+...+an​λn​vn​
$x_2= a_1\lambda_1^2v_1+a_2\lambda_2^2v_2+...+a_n\lambda_n^2v_n$x2​=a1​λ12​v1​+a2​λ22​v2​+...+an​λn2​vn​
$\dots$　　　　　　…
$x_k= a_1\lambda_1^kv_1+a_2\lambda_2^kv_2+...+a_n\lambda_n^kv_n$xk​=a1​λ1k​v1​+a2​λ2k​v2​+...+an​λnk​vn​

下面按模最大特征值λ1是单根的情况讨论:
$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​λ1k​v1​+a2​λ2k​v2​+...+an​λnk​vn​=λ1k​[a1​v1​+a2​(λ1​λ2​​)kv2​+...+an​(λ1​λn​​)2vn​]=λ1k​(a1​v1​+εk​)

$\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​
由于$|\lambda_1|&gt;|\lambda_2|&gt;...&gt;|\lambda_n|　　故|\frac{\lambda_j}{\lambda_1}|&lt;1(j=2,3,...n),故：$∣λ1​∣>∣λ2​∣>...>∣λn​∣　　故∣λ1​λj​​∣<1(j=2,3,...n),故：
$\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​λ1k​1​xk​=k→+∞lim​λ1k​λ1k​(a1​v1​+εk​)​=a1​v1​

$\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^k}x_k=a_1v_1 就是\lambda_1对应的特征向量$k→+∞lim​λ1k​1​xk​=a1​v1​就是λ1​对应的特征向量
$用向量x_{k+1}的第i个分量与向量x_k的第i个分量之比的极限：$用向量xk+1​的第i个分量与向量xk​的第i个分量之比的极限：
$\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​(a1​v1​+εk​)i​λ1k+1​(a1​v1​+εk+1​)i​​=k→+∞lim​a1​(v1​)i​+(εk​)i​λ1​[a1​(v1​)i​+(εk+1​)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​λ1k​1​xk​=a1​v1​　　　k→+∞lim​(xk​)i​(xk+1​)i​​=λ1​

因此当k趋近无穷时，就能近似得到按模最大的特征值和对应的特征向量：
$\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​≈λ1k​a1​v1​

因此得幂方法的迭代公式：$x_k = Ax_{k-1}　　k=1,2,...$xk​=Axk−1​　　k=1,2,...
当k充分大时，有：
$\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​)i​​xk​≈λ1k​a1​v1​​

存在的问题是：当$|\lambda_1|&gt;1时，x_k中不为0的分量会随着k\rightarrow \infty而趋于 \infty，计算机计算时会造成溢出，当$∣λ1​∣>1时，xk​中不为0的分量会随着k→∞而趋于∞，计算机计算时会造成溢出，当|\lambda_1|<1时,x_k中的分量会随着k\rightarrow \infty而趋于0$，

实际计算时，为了避免计算过程中出现绝对值过大或过小的数参加运算，通常在每步迭代时，将向量“归一化”即用的按模最大的分量 。
修改计算公式：
$\begin{cases} y_k=Ax_{k-1}\\ m_k=max(y_k) &amp; (k=1,2,...)\\ x_k=y_k/m_k \end{cases}$⎩⎪⎨⎪⎧​yk​=Axk−1​mk​=max(yk​)xk​=yk​/mk​​(k=1,2,...)​
上式中$m_k是向量y_k中绝对值最大的一个分量，这时x_k分量的模最大为1$mk​是向量yk​中绝对值最大的一个分量，这时xk​分量的模最大为1
当k充分大时，有：
$\begin{cases} \lambda_1 \approx m_k \\ v_1 \approx x_k \end{cases}${λ1​≈mk​v1​≈xk​​

实例

求矩阵
$\left[ \begin{matrix} 2 &amp;3 &amp; 2 \\ 10 &amp; 3 &amp; 4 \\ 3 &amp; 6 &amp; 1 \end{matrix} \right] 的主特征值和其对应的特征向量$⎣⎡​2103​336​241​⎦⎤​的主特征值和其对应的特征向量
根据迭代公式：
$\begin{cases} y_k=Ax_{k-1}\\ m_k=max(y_k) &amp; (k=1,2,...)\\ x_k=y_k/m_k \end{cases}$⎩⎪⎨⎪⎧​yk​=Axk−1​mk​=max(yk​)xk​=yk​/mk​​(k=1,2,...)​
取$x_0=(0,0,1)^T,则有：$x0​=(0,0,1)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​=m1​1​y1​=(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]]
===================================