幂法

特征值与特征向量的计算

张瑞
中国科学技术大学数学科学学院

幂法

若有数$\lambda$及$n$维非零向量$v$满足

\[Ax=\lambda x \]

则称$\lambda$是矩阵$A$的特征值，$v$称为特征向量。

矩阵的特征值是经常需要计算的量： 如迭代矩阵是否收敛矩阵取决于它的谱半径， Google搜索中，使用特征值来决定网页的相关性。

  • 在线性代数中，计算矩阵$A$的特征值需要计算矩阵的特征多项式$\det(\lambda I-A)$的根。
  • 这是一个高次多项式，它的求根非常困难。因而，需要高效的数值方法。

幂法

  • 矩阵的按模最大特征值往往比其它的特征值具有更重要的实际意义。
  • 幂法就是用来计算矩阵的按模最大特征值和相应的特征向量。

这里，需要矩阵$A$有$n$个线性无关的特征向量$v_1$, $v_2$, $\cdots$, $v_n$， 即$A$是完备的。

以一个非0向量$x^{(0)}$开始，构造序列

\[x^{(k+1)}=A x^{(k)} \]

显然序列有通项公式$x^{(k)}=A^k x^{(0)}$。下面，看一下序列$\{x^{(k)}\}$的表现。

设初值$x^{(0)}$可以表达为

\[x^{(0)}=\alpha_1v_1+\alpha_2v_2+\cdots+\alpha_nv_n \]

则有通项

\[\begin{aligned} x^{(k)}=&A^nx^{(0)}=A^k(\alpha_1v_1+\alpha_2v_2+\cdots+\alpha_nv_n) \\ =&\alpha_1 A^kv_1+\alpha_2A^kv_2+\cdots+\alpha_nA^kv_n \\ \end{aligned} \]

由$v_i$均为特征向量，因而

\[\begin{aligned} A^k v_i=&A^{k-1}\lambda_i v_i=\lambda_i A^{k-1}v_i \\ =&\lambda_iA^{k-2}\lambda_iv_i =\lambda_i^2 A^{k-2}v_i \\ =&\cdots=\lambda_i^k v_i \end{aligned} \]

因而

\[x^{(k)}=\lambda_1^k\alpha_1v_1+\lambda_2^k\alpha_2v_2+\cdots+\lambda_n^k\alpha_nv_n \]

不防假定

\[|\lambda_1|\geq|\lambda_2|\geq\cdots\geq|\lambda_n| \]

则有

\[x^{(k)}=\lambda_1^k\left(\alpha_1v_1+\left(\frac{\lambda_2}{\lambda_1}\right)^k\alpha_2v_2 +\cdots+\left(\frac{\lambda_n}{\lambda_1}\right)^k\alpha_nv_n\right) \]

(1) 若$|\lambda_1|>|\lambda_2|\geq\cdots\geq|\lambda_n|$

此时，$|\frac{\lambda_i}{\lambda_1}|<1$, $i=2,\cdots,n$。 因而，随着$k$的不断增大，$\left(\frac{\lambda_i}{\lambda_1}\right)^k$不断接近0。

当$k$充分大后，

\[\begin{aligned} x^{(k)}=&\lambda_1^k\left(\alpha_1v_1+\left(\frac{\lambda_2}{\lambda_1}\right)^k\alpha_2v_2 +\cdots+\left(\frac{\lambda_n}{\lambda_1}\right)^k\alpha_nv_n\right) \\ \approx& \lambda_1^k \alpha_1v_1 \end{aligned} \]

因而，有

\[x^{(k+1)}=A x^{(k)} \approx\lambda_1^{k+1}\alpha_1 v_1=\lambda_1 \lambda_1^k\alpha_1v_1=\lambda_1 x^{(k)} \]

也就是说，$x^{(k)}$近似为$\lambda_1$的特征向量，且$x^{(k)}$与$x^{(k+1)}$共线。

因而，对$x^{(k)}$的任意分量

\[\lambda_1\approx \frac{x^{(k+1)}_i}{x^{(k)}_i}, i=1,2,\cdots,n \]

依特征向量的定义，此时$x^{(k)}$和$x^{(k+1)}$均为$\lambda_1$的特征向量。

(2) 若$|\lambda_1|=|\lambda_2|>|\lambda_3|\geq\cdots\geq|\lambda_n|$，且$\lambda_1=-\lambda_2$

当$k$充分大后，

\[\begin{aligned} x^{(k)}=&\lambda_1^k\left(\alpha_1v_1+\left(-1\right)^k\alpha_2v_2 +\cdots+\left(\frac{\lambda_n}{\lambda_1}\right)^k\alpha_nv_n\right) \\ \approx& \lambda_1^k (\alpha_1v_1+(-1)^k\alpha_2v_2) \end{aligned} \]

从而

\[\begin{aligned} x^{(k+1)} \approx & \lambda_1^{k+1} (\alpha_1v_1+(-1)^{k+1}\alpha_2v_2) \\ x^{(k+2)} \approx & \lambda_1^{k+2} (\alpha_1v_1+(-1)^{k+2}\alpha_2v_2) \\ \approx & \lambda_1^2 x^{(k)} \end{aligned} \]

可以看到，$x^{(k)}$与$x^{(k+2)}$共线，且对$x^{(k)}$的任意分量

\[\lambda_1=\sqrt{\frac{x^{(k+2)}_i}{x^{(k)}_i}} , i=1,2,\cdots,n \]

另外，由

\[\begin{aligned} \lambda_1 x^{(k+1)} \approx & \lambda_1^{k+2} (\alpha_1v_1+(-1)^{k+1}\alpha_2v_2) \\ x^{(k+2)} \approx & \lambda_1^{k+2} (\alpha_1v_1+(-1)^{k+2}\alpha_2v_2) \\ \end{aligned} \]

得到

\[\begin{aligned} 2 \lambda_1^{k+2} \alpha_1v_1=x^{(k+2)}+\lambda_1 x^{(k+1)} \\ 2(-1)^k \lambda_1^{k+2} \alpha_2v_2=x^{(k+2)}-\lambda_1 x^{(k+1)} \\ \end{aligned} \]

依特征向量的性质，相应的特征向量可以取为

\[\begin{aligned} v_1=x^{(k+2)}+\lambda_1 x^{(k+1)} \\ v_2=x^{(k+2)}-\lambda_1 x^{(k+1)} \\ \end{aligned} \]

(3) 特征值的其它分布情况，如： $\lambda_1=\bar\lambda_2$，请参考相关书籍， 这里不再考虑

幂法的算法为：

  1. 以非0向量$x^{(0)}$，得到序列$x^{(k)}=A x^{(k-1)}$
  2. 若$x^{(k+1)}$与$x^{(k)}$共线，则
    1. $\lambda_1=\frac{x^{(k+1)}_i}{x^{(k)}_i}$(任取分量), $v_1=x^{(k+1)}$
  3. 若$x^{(k+1)}$与$x^{(k-1)}$共线，则
    1. $\lambda_1=\sqrt{\frac{x^{(k+1)}_i}{x^{(k-1)}_i}}$(任取分量), $\lambda_2=-\lambda_1$, $v_1=x^{(k+1)}+\lambda_1 x^{(k)}$, $v_2=x^{(k+1)}-\lambda_1 x^{(k)}$
  4. 若序列是其它的情况，则退出。 
    for i in range(maxsteps):
        x2=np.dot(mat,x1)
        with np.errstate(divide='ignore'): # 去掉除0的警告
            x21=x2/x1  # x^{(k+1)}/x^{(k)}
            x20=x2/x0  # x^{(k+1)}/x^{(k-1)}
        steps.append((i+2,[x2,x21,x20]))
        y0=[]
        for y00 in x21:
            if np.isnan(y00): continue
            y0.append(y00)

        if np.array(y0).ptp()<eps: #计算最大值与最小值差
            result=[1,np.array(y0).max(),[x2]]
            break

        y0=[]
        for y00 in x20:
            if np.isnan(y00): continue
            y0.append(y00)
        if np.array(y0).ptp()<eps:
            result=[2,np.sqrt(np.array(y0).max()),[x1,x2]]
            break;
 
        x0=x1
        x1=x2

例 1. 求解矩阵

\[A=\begin{pmatrix} -2. & 1. & 0. \\ 1. & 2. & 0. \\ 1. & 0. & 1.5 \end{pmatrix}, \quad\, B=\begin{pmatrix} -9. & 1. & 3. \\ 2. & 7. & 0. \\ 3. & 0. & 7. \end{pmatrix} \]

的按模最大特征值和相应的特征向量

解. 若干步计算后，得到$A$的计算结果

After  64  steps
---- Step: 61
[4.65661287e+20 9.31322575e+20 2.53997066e+20]
跟前1步的比值      [inf 2.  1.5]
跟前2步的比值      [5. 5. 5.]
---- Step: 62
[0.00000000e+00 2.32830644e+21 8.46656886e+20]
跟前1步的比值      [0.         2.5        3.33333333]
跟前2步的比值      [nan  5.  5.]
---- Step: 63
[2.32830644e+21 4.65661287e+21 1.26998533e+21]
跟前1步的比值      [inf 2.  1.5]
跟前2步的比值      [5. 5. 5.]
---- Step: 64
[0.00000000e+00 1.16415322e+22 4.23328443e+21]
跟前1步的比值      [0.         2.5        3.33333333]
跟前2步的比值      [nan  5.  5.]
---- Step: 65
[1.16415322e+22 2.32830644e+22 6.34992665e+21]
跟前1步的比值      [inf 2.  1.5]
跟前2步的比值      [5. 5. 5.]

根据数据，可以判断特征值分布为？特征向量是？

按模最大有互为反号的2个
2.23606797749979
特征向量：
[array([1.16415322e+22, 2.32830644e+22, 6.34992665e+21]),
 array([0.00000000e+00, 5.82076609e+22, 2.11664222e+22])]

$B$的计算结果 

After  125  steps
---- Step: 123
[ 2.87805266e+120 -3.45499252e+119 -5.18248878e+119]
跟前1步的比值      [-9.66025404 -9.66025404 -9.66025404]
跟前2步的比值      [93.32050808 93.32050808 93.32050808]
---- Step: 124
[-2.78027198e+121  3.33761055e+120  5.00641582e+120]
跟前1步的比值      [-9.66025404 -9.66025404 -9.66025404]
跟前2步的比值      [93.32050808 93.32050808 93.32050808]
---- Step: 125
[ 2.68581336e+122 -3.22421658e+121 -4.83632486e+121]
跟前1步的比值      [-9.66025404 -9.66025404 -9.66025404]
跟前2步的比值      [93.32050808 93.32050808 93.32050808]
---- Step: 126
[-2.59456394e+123  3.11467512e+122  4.67201268e+122]
跟前1步的比值      [-9.66025404 -9.66025404 -9.66025404]
跟前2步的比值      [93.32050808 93.32050808 93.32050808]

根据数据，可以判断特征值分布为？特征向量是？

按模最大只有一个
-9.660254037755953
特征向量：
[array([ 2.50641467e+124, -3.00885529e+123, -4.51328293e+123])]

由

\[\begin{aligned} x^{(k)}=&\lambda_1^k\left(\alpha_1v_1+\left(\frac{\lambda_2}{\lambda_1}\right)^k\alpha_2v_2 +\cdots+\left(\frac{\lambda_n}{\lambda_1}\right)^k\alpha_nv_n\right) \\ \end{aligned} \]

可以看到，

  • 若$|\lambda_1|>1$，则$k$充分大以后，$x^{(k)}$的每个分量都会很大。
  • 若$|\lambda_1|<1$，则$k$充分大以后，$x^{(k)}$的每个分量都会接近0。

规范的幂法

幂法在计算中，很可能造成计算的溢出（上溢出到$\infty$，或下溢出归0）。

规范的幂法运算就是在每一步计算中添加一个规范化的操作

\[\begin{cases} x^{(k+1)}=Ay^{(k)} \\ y^{(k+1)}=\frac{x^{(k+1)}}{\|x^{(k+1)}\|_{\infty}} \end{cases}, \quad k=0,1,2,\cdots \]

给定初值$x^{(0)}$后，取$y^{(0)}=\frac{x^{(0)}}{\|x^{(0)}\|_{\infty}}$， 代入上式构造序列。

得到的序列$y^{(k)}$满足$\|y^{(k)}\|_{\infty}=1$，它永远不会溢出。

依算法，有

\[\begin{aligned} y^{(k)}=&\frac{x^{(k)}}{\|x^{(k)}\|_{\infty}} =\frac{Ay^{(k-1)}}{\|x^{(k)}\|_{\infty}} =\frac{1}{\|x^{(k)}\|_{\infty}}Ay^{(k-1)} \\ =&\frac1{\|x^{(k)}\|_{\infty}}A\frac{x^{(k-1)}}{\|x^{(k-1)}\|_{\infty}} =\frac1{\|x^{(k)}\|_{\infty}}A\frac{Ay^{(k-2)}}{\|x^{(k-1)}\|_{\infty}} \\ =&\frac1{\|x^{(k)}\|_{\infty}\|x^{(k-1)}\|_{\infty}}A^2y^{(k-2)} \\ =&\cdots \\ =&\frac1{\|x^{(k)}\|_{\infty}\cdots\|x^{(1)}\|_{\infty}}A^ky^{(0)} \end{aligned} \]

由$\|y^{(k)}\|_{\infty}=1$，因而$y^{(k)}$可以写成

\[y^{(k)}=\frac{A^ky^{(0)}}{\|A^ky^{(0)}\|_{\infty}} \]

不防设

\[y^{(0)}=\alpha_1 v_1+\cdots +\alpha_nv_n \]

则

\[A^ky^{(0)}=\lambda_1^k\left(\alpha_1v_1+\left(\frac{\lambda_2}{\lambda_1}\right)^k\alpha_2v_2 +\cdots+\left(\frac{\lambda_n}{\lambda_1}\right)^k\alpha_nv_n\right) \]

(1) 若$|\lambda_1|>|\lambda_2|\geq\cdots\geq|\lambda_n|$

则$k$充分大以后，有

\[A^ky^{(0)}\approx \lambda_1^k \alpha_1v_1 \]

从而

\[y^{(k)}\approx\frac{\lambda_1^k \alpha_1v_1}{\|\lambda_1^k \alpha_1v_1\|_{\infty}} =\frac{\lambda_1^k}{|\lambda_1|^k}\frac{\alpha_1v_1}{\|\alpha_1v_1\|_{\infty}} \]
  • 若$\lambda_1>0$，则有$y^{(k)}=\frac{\alpha_1v_1}{\|\alpha_1v_1\|_{\infty}}$， 即$y^{(k)}$序列收敛到一个向量。
  • 若$\lambda_1<0$，则有$y^{(k)}=(-1)^k\frac{\alpha_1v_1}{\|\alpha_1v_1\|_{\infty}}$， 即$y^{(k)}$的奇序列和偶序列收敛各自收敛，且收敛到的向量互为反方向。

此时，$y^{(k)}$就是特征向量。利用它可以得到特征值。

注意到

\[x^{(k+1)}=Ay^{(k)}=\lambda_1 y^{(k)} \]

因而

\[|\lambda_1|=\|x^{(k+1)}\|_{\infty} \]

$\lambda_1$的符号由$y^{(k)}$序列的表现决定。

(2) 若$|\lambda_1|=|\lambda_2|>|\lambda_3|\geq\cdots\geq|\lambda_n|$， 且$\lambda_1=-\lambda_2$

则$k$充分大以后，有

\[A^ky^{(0)}\approx \lambda_1^k (\alpha_1v_1+(-1)^k\alpha_2v_2) \]

从而

\[y^{(k)}\approx\frac{\lambda_1^k (\alpha_1v_1+(-1)^k\alpha_2v_2)}{\|\lambda_1^k (\alpha_1v_1+(-1)^k\alpha_2v_2)\|_{\infty}} \]

因此，

  • $y^{(k)}$的奇序列和偶序列收敛各自收敛，且收敛到的向量不共线。

此时，需要借助幂法的操作来得到特征值和特征向量。

不防记

\[y^{(k)}\approx \lambda_1^k (\tilde\alpha_1v_1+(-1)^k\tilde\alpha_2v_2) \]

则做2步幂法，

\[\begin{aligned} w^{(1)}=A y^{(k)}= & \lambda_1^{k+1} (\tilde\alpha_1v_1+(-1)^{k+1}\tilde\alpha_2v_2) \\ w^{(2)}=A w^{(1)} = & \lambda_1^{k+2} (\tilde\alpha_1v_1+(-1)^{k+2}\tilde\alpha_2v_2) \\ =&\lambda_1^2 y^{(k)} \end{aligned} \]

类似与幂法时的分析，可以得到

此时$w^{(2)}$与$y^{(k)}$是共线的，取任意分量计算可以得到

\[\lambda_1=\sqrt{\frac{w^{(2)}_i}{y^{(k)}_i}}=\sqrt{\|w^{(2)}\|_{\infty}} \]

及

\[\begin{aligned} v_1=& w^{(2)}+\lambda_1w^{(1)} \\ v_2=& w^{(2)}-\lambda_1w^{(1)} \\ \end{aligned} \]

规范幂法的算法：

  1. 以非0向量$x^{(0)}$，得到序列$y^{(k)}=\frac{A x^{(k-1)}}{\|A x^{(k-1)}\|_{\infty}}$
  2. 若$y^{(k)}$序列收敛，则
    1. 作$x^{(k+1)}=Ay^{(k)}$，有$\lambda_1=\|x^{(k+1)}\|_{\infty}$, $v_1=y^{(k)}$
  3. 若$y^{(k)}$的奇、偶序列各自收敛，则
    1. 若收敛的两个向量互为反方向，则作$x^{(k+1)}=Ay^{(k)}$，有$\lambda_1=-\|x^{(k+1)}\|_{\infty}$, $v_1=y^{(k)}$
    2. 若收敛的两个向量不共线，则作
      \[w^{(1)}=A y^{(k)}, \quad w^{(2)}=A w^{(1)} \]
      有 $\lambda_1=\sqrt{\frac{w^{(2)}_i}{y^{(k)}_i}}=\sqrt{\|w^{(2)}\|_{\infty}}$, $\lambda_2=-\lambda_1$, $v_1= w^{(2)}+\lambda_1w^{(1)}$,$v_2= w^{(2)}-\lambda_1w^{(1)}$
  4. 若序列是其它的情况，则退出。

例 2. 用规范的幂法求解矩阵

\[A=\begin{pmatrix} -2. & 1. & 0. \\ 1. & 2. & 0. \\ 1. & 0. & 1.5 \end{pmatrix}, \quad\, B=\begin{pmatrix} -9. & 1. & 3. \\ 2. & 7. & 0. \\ 3. & 0. & 7. \end{pmatrix} \]
\[C=\begin{pmatrix} 9. & 1. & 3. \\ 2. & 7. & 0. \\ 3. & 0. & 8. \end{pmatrix} \]

的按模最大特征值和相应的特征向量

解. 若干步计算后，得到$A$的计算结果

After  58  steps
---- Step: 55
[array([0.5       , 1.        , 0.27272727])]
---- Step: 56
[array([0.        , 1.        , 0.36363636])]
---- Step: 57
[array([0.5       , 1.        , 0.27272727])]
---- Step: 58
[array([0.        , 1.        , 0.36363636])]
---- Step: 59
[array([0.5       , 1.        , 0.27272727])]

根据数据，可以判断特征值分布为？特征向量是？

按模最大有互为反号的2个
2.23606797749979
特征向量：
[array([1.        , 2.        , 0.54545455]),
array([0.        , 5.        , 1.81818182])]

与幂法的结果比较 

按模最大有互为反号的2个
2.23606797749979
特征向量：
[array([1.16415322e+22, 2.32830644e+22, 6.34992665e+21]),
 array([0.00000000e+00, 5.82076609e+22, 2.11664222e+22])]

$B$的计算结果 

Regular Power method:
---- Step: 99
[array([ 1.        , -0.12004619, -0.18006928])]
---- Step: 100
[array([-1.        ,  0.12004619,  0.18006928])]
---- Step: 101
[array([ 1.        , -0.12004619, -0.18006928])]
---- Step: 102
[array([-1.        ,  0.12004619,  0.18006928])]
---- Step: 103
[array([ 1.        , -0.12004619, -0.18006928])]

根据数据，可以判断特征值分布为？特征向量是？

按模最大只有一个且<0
-9.66025403842121
特征向量：
[-1.          0.12004619  0.18006928]

与幂法的结果相比 

按模最大只有一个
-9.660254037755953
特征向量：
[array([ 2.50641467e+124, -3.00885529e+123, -4.51328293e+123])]

$C$的计算结果 

Regular Power method:
---- Step: 40
[array([1.        , 0.4174243 , 0.79128785])]
---- Step: 41
[array([1.        , 0.4174243 , 0.79128785])]
---- Step: 42
[array([1.        , 0.4174243 , 0.79128785])]
---- Step: 43
[array([1.        , 0.4174243 , 0.79128785])]
---- Step: 44
[array([1.        , 0.4174243 , 0.79128785])]

根据数据，可以判断特征值分布为？特征向量是？

按模最大只有一个且>0
11.791287847526174
特征向量：
[1.         0.4174243  0.79128785]

与幂法的结果相比 

按模最大只有一个
11.79128784753576
特征向量：
[array([2.75104171e+55, 1.14835167e+55, 2.17686587e+55])]

例 3. 求矩阵

\[D=\begin{pmatrix} 0.35 & 0.001& 0.1 \\ 0.1 & 0.2 & 0.13 \\ 0. & 0. & 0.34 \\ \end{pmatrix} \]

的按模最大特征值和相应的特征向量

解. 用Python的 numpy.linalg.eig 得到的 

特征值： [0.35066373 0.19933627 0.34      ]

用幂法的计算结果

幂法 

到达步数上限
---- Step: 497
[6.73673188e-226 4.47136945e-226 1.39643836e-233]
跟前1步的比值      [0.35066373 0.35066373 0.34      ]
跟前2步的比值      [0.12296505 0.12296505 0.1156    ]
---- Step: 498
[2.36232754e-226 1.56794710e-226 4.74789042e-234]
跟前1步的比值      [0.35066373 0.35066373 0.34      ]
跟前2步的比值      [0.12296505 0.12296505 0.1156    ]
---- Step: 499
[8.28382592e-227 5.49822179e-227 1.61428274e-234]
跟前1步的比值      [0.35066373 0.35066373 0.34      ]
跟前2步的比值      [0.12296505 0.12296505 0.1156    ]
---- Step: 500
[2.90483731e-227 1.92802697e-227 5.48856133e-235]
跟前1步的比值      [0.35066373 0.35066373 0.34      ]
跟前2步的比值      [0.12296505 0.12296505 0.1156    ]

规范的幂法 

Regular Power method:
---- Step: 553
[array([1.00000000e+00, 6.63729754e-01, 3.67709048e-09])]
---- Step: 554
[array([1.00000000e+00, 6.63729754e-01, 3.56526968e-09])]
---- Step: 555
[array([1.00000000e+00, 6.63729754e-01, 3.45684936e-09])]
---- Step: 556
[array([1.00000000e+00, 6.63729754e-01, 3.35172612e-09])]
---- Step: 557
[array([1.00000000e+00, 6.63729754e-01, 3.24979968e-09])]
按模最大只有一个且>0
0.35066373006863033
特征向量：
[1.00000000e+00 6.63729754e-01 3.15097285e-09]

由

\[A^ky^{(0)}=\lambda_1^k\left(\alpha_1v_1+\left(\frac{\lambda_2}{\lambda_1}\right)^k\alpha_2v_2 +\cdots+\left(\frac{\lambda_n}{\lambda_1}\right)^k\alpha_nv_n\right) \]

可以看到，规范幂法的收敛速度由$\frac{\lambda_2}{\lambda_1}$或$\frac{\lambda_3}{\lambda_1}$的大小决定， 它越接近0，收敛越快。

前面算例中的数值实验结果也可以看到这点。

特征值 $\frac{\lambda_2}{\lambda_1}$ 收敛步数
([ 1.5 , 2.23606798, -2.23606798]) 0.67 58
([11.79128785, 5. , 7.20871215]) 0.62 43
([-9.66025404, 7.66025404, 7. ]) 0.79 102
([0.19933627 , 0.34 , 0.35066373]) 0.97 556

可以使用原点位移法来加速幂法的收敛。

注意到矩阵$B=A-pI$的特征值与矩阵$A$的特征值相差$p$，即$\lambda_B=\lambda_A-p$。 因而，可以选取合适的$p$，使得

\[\left|\frac{\lambda_2-p}{\lambda_1-p}\right|<\left|\frac{\lambda_2}{\lambda_1}\right| \]

这样，对矩阵$B$的幂法的收敛速度会比矩阵$A$的幂法快。

  • $p$的选取需要保证$\lambda_1-p$也是$B$的按模最大特征值。
  • 只有对矩阵的特征值分布有个大致了解，才能选取有效的$p$。

例 4. 求矩阵

\[D=\begin{pmatrix} 0.35 & 0.001& 0.1 \\ 0.1 & 0.2 & 0.13 \\ 0. & 0. & 0.34 \\ \end{pmatrix}, \quad\, E=\begin{pmatrix} 0.15 & 0.001& 0.1 \\ 0.1 & 0.0 & 0.13 \\ 0. & 0.0 & 0.14 \\ \end{pmatrix} \]

的按模最大特征值和相应的特征向量

解. 可以看到$D$的特征值与$E$的相差$0.2$，它们的收敛步数为 

(array([0.19933627 ,  0.34      ,  0.35066373]), 步数: 556, 比值: 0.97)
(array([-0.00066373,  0.14      ,  0.15066373]), 步数: 245, 比值: 0.93)

问题. 若初值

\[x^{(0)}=\alpha_1v_1+\alpha_2v_2+\cdots+\alpha_nv_n \]

中$\alpha_1=0$， 会对计算结果产生怎样的影响？

反幂法

反幂法用来求解矩阵的按模最小的特征值和相应的特征向量。

若矩阵$A$可逆，$\lambda$的$v$分别是$A$的特征值和相应的特征向量，则

\[\begin{aligned} &Av = \lambda v \\ \Rightarrow\quad & \frac1\lambda A^{-1}Ax=\frac1\lambda A^{-1}\lambda v \\ \Rightarrow\quad & A^{-1}v=\frac1\lambda v \\ \end{aligned} \]

即，$A$与$A^{-1}$的特征值互为倒数。

这样，通过求$A^{-1}$的按模最大特征值就可以得到$A$的按模最小特征值。

迭代序列

\[\begin{cases} x^{(k+1)}=A^{-1}y^{(k)} \\ y^{(k+1)}=\frac{x^{(k+1)}}{\|x^{(k+1)}\|_{\infty}} \end{cases}, k=0,1,\cdots \]

就是规范的幂法运算求解$A^{-1}$的按模最大特征值的序列。

通常用解线性方程组来避免求解矩阵的逆，

\[\begin{cases} Ax^{(k+1)}=y^{(k)} \\ y^{(k+1)}=\frac{x^{(k+1)}}{\|x^{(k+1)}\|_{\infty}} \end{cases}, k=0,1,\cdots \]

这里，一般用LU分解法来解方程组。

若$A$的特征值为$|\lambda_1|>\cdots>|\lambda_{n-1}|>|\lambda_n|$， 则反幂法的收敛速度由

\[\frac{\frac1{\lambda_{n-1}}}{\frac1{\lambda_n}}=\frac{\lambda_n}{\lambda_{n-1}} \]

的大小来决定。

可以看到，$\lambda_n$越接近0，收敛速度越快。

例 5. 求矩阵

\[D=\begin{pmatrix} 0.35 & 0.001& 0.1 \\ 0.1 & 0.2 & 0.13 \\ 0. & 0. & 0.34 \\ \end{pmatrix}, \quad\, E=\begin{pmatrix} 0.15 & 0.001& 0.1 \\ 0.1 & 0.0 & 0.13 \\ 0. & 0.0 & 0.14 \\ \end{pmatrix} \]

的按模最大特征值和相应的特征向量

解. 可以看到$D$的特征值与$E$的相差$0.2$，它们的收敛步数为 

(特征值([0.19933627 ,  0.34      ,  0.35066373]), 步数: 45)
(特征值([-0.00066373,  0.14      ,  0.15066373]), 步数: 5)

矩阵$D$的运行结果 

特征值： [0.35066373 0.19933627 0.34      ]
特征向量：
 [[ 0.83317794 -0.00663715 -0.84807616]
 [ 0.55300499  0.99997797 -0.52217004]
 [ 0.          0.          0.09002932]]
[12.50144587 -0.61337544 11.10749297]
After  45  steps
按模最小只有一个且>0
0.1993362702596363
特征向量：
[ 6.63729738e-03 -1.00000000e+00  2.05746046e-11]

矩阵$E$的运行结果， 

特征值： [ 0.15066373 -0.00066373  0.14      ]
特征向量：
 [[ 0.83317794 -0.00663715 -0.84807616]
 [ 0.55300499  0.99997797 -0.52217004]
 [ 0.          0.          0.09002932]]
[12.50144587 -0.61337544 11.10749297]
After  5  steps
按模最小只有一个且<0
-0.0006637297521077966
特征向量：
[-6.63729752e-03  1.00000000e+00  8.77659247e-17]

将反幂法用到矩阵$A-pI$上，

  • 则可以得到$A-pI$的按模最小特征值$\lambda_p$，
  • 从而，$A$的特征值中最接近$p$的特征值是$p+\lambda_p$，设为$\lambda_i$，
  • 而且$\lambda_p$越接近0，收敛速度越快。 即对$A$的特征值$\lambda_i$的估计值$p$越精确，收敛速度越快。

例 6. 求矩阵

\[F=\begin{pmatrix} 0.05 & 0.001& 0.1 \\ 0.1 & -0.1 & 0.13 \\ 0. & 0. & 0.04 \\ \end{pmatrix}, \quad\, G=\begin{pmatrix} 0.02 & 0.001& 0.1 \\ 0.1 & -0.13 & 0.13 \\ 0. & 0.0 & 0.01 \\ \end{pmatrix} \]

的按模最大特征值和相应的特征向量

解. 可以看到$F$,$G$的特征值与前面的$D$的特征值相差$0.3$,$0.31$，它们的收敛步数为 

(特征值([ 0.05066373 -0.10066373  0.04      ]), 步数: 81)
(特征值([ 0.02066373 -0.13066373  0.01      ]), 步数: 27)

矩阵$F$的运行结果 

特征值： [ 0.05066373 -0.10066373  0.04      ]
特征向量：
 [[ 0.83317794 -0.00663715 -0.84807616]
 [ 0.55300499  0.99997797 -0.52217004]
 [ 0.          0.          0.09002932]]
[12.50144587 -0.61337544 11.10749297]
=================================
Inverse Power method:
After  81  steps
按模最小只有一个且>0
0.039999999964299746
特征向量：
[-1.         -0.61571125  0.10615711]

矩阵$G$的运行结果， 

特征值： [ 0.02066373 -0.13066373  0.01      ]
特征向量：
 [[ 0.83317794 -0.00663715 -0.84807616]
 [ 0.55300499  0.99997797 -0.52217004]
 [ 0.          0.          0.09002932]]
[12.50144587 -0.61337544 11.10749297]
=================================
Inverse Power method:
After  27  steps
按模最小只有一个且>0
0.009999999991478744
特征向量：
[-1.         -0.61571125  0.10615711]

目录

本节读完

例 7.

7.