四、矩阵特征值与特征向量的计算

1 幂法和反幂法

1.1 幂法

幂法主要用于求矩阵的按模[1]最大的特征值与相应的特征向量。它是通过迭代产生向量序列，由此计算特征值和特征向量。

设n×nn\times nn×n阶实矩阵AAA的特征值λi(i=1,2,…,n)\lambda_i\ (i=1,2,\dots,n)λi (i=1,2,…,n)满足：

且与λi(i=1,2,…,n)\lambda_i(i=1,2,\dots,n)λi(i=1,2,…,n)相应的特征向量u1,u2,…,unu_1,u_2,\dots,u_nu1,u2,…,un线性无关。（？为什么线性无关）

给定初始向量x(0)≠0x^{(0)}\ne0x(0)=0，由迭代公式

x(k+1)=Ax(k)(k=0,1,2,…)(4−1)x^{(k+1)}=Ax^{(k)}\quad(k=0,1,2,\dots)\quad\quad\quad\quad\quad(4-1)x(k+1)=Ax(k)(k=0,1,2,…)(4−1)

产生向量序列{x(k)}\{x^{(k)}\}{x(k)}，可以证明，当kkk充分大时，有λi≈xi(k+1)/xi(k)\lambda_i\approx x_i^{(k+1)}/x_i^{(k)}λi≈xi(k+1)/xi(k)，相应的特征向量为x(k+1)x^{(k+1)}x(k+1)。

为简便起见，不妨设∥ui∥=1(i=1,2,…,n)\Vert u_i\Vert=1\ (i=1,2,\dots,n)∥ui∥=1 (i=1,2,…,n)。因为uiu_iui线性无关（线性无关的概念），故必存在n个不全为零的数αi(i=1,2,…,n)\alpha_i\ (i=1,2,\dots,n)αi (i=1,2,…,n)，使得x(0)=∑i=1nαiuix^{(0)}=\sum\limits_{i=1}^n\alpha_iu_ix(0)=i=1∑nαiui。由式（4-1）得：

x(k+1)=Ax(k)=AK+1x(0)=∑i=1nAk+1(αiui)x^{(k+1)}=Ax^{(k)}=A^{K+1}x^{(0)}=\sum\limits_{i=1}^nA^{k+1}(\alpha_iu_i)x(k+1)=Ax(k)=AK+1x(0)=i=1∑nAk+1(αiui)

即：

x(k+1)=λ1k+1[α1u1+(λ2λ1)α2u2+⋯+(λnλ1)αnun](4−2)x^{(k+1)}=\lambda_1^{k+1}[\alpha_1u_1+(\frac{\lambda_2}{\lambda_1})\alpha_2u_2+\cdots+(\frac{\lambda_n}{\lambda_1})\alpha_nu_n]\quad\quad\quad\quad\quad(4-2)x(k+1)=λ1k+1[α1u1+(λ1λ2)α2u2+⋯+(λ1λn)αnun](4−2)

设αi≠0\alpha_i\ne0αi=0，由∣λ1∣>∣λi∣(i=2,3,…,n)|\lambda_1|>|\lambda_i|\ (i=2,3,\dots,n)∣λ1∣>∣λi∣ (i=2,3,…,n)得：

lim⁡k→∞(λiλ1)k+1αiui=0\lim\limits_{k\to\infty}(\frac{\lambda_i}{\lambda_1})^{k+1}\alpha_iu_i=0k→∞lim(λ1λi)k+1αiui=0

于是：

lim⁡k→∞∑i=2n(λiλ1)αiui=0\lim\limits_{k\to\infty}\sum\limits_{i=2}^n(\frac{\lambda_i}{\lambda_1})\alpha_iu_i=0k→∞limi=2∑n(λ1λi)αiui=0

故只要k充分大，就有：

x(k+1)=λ1k+1[α1u1+∑i=2n(λiλ1)αiui]≈λ1k+1α1u1(4−3)x^{(k+1)}=\lambda_1^{k+1}[\alpha_1u_1+\sum\limits_{i=2}^n(\frac{\lambda_i}{\lambda_1})\alpha_iu_i]\approx\lambda_1^{k+1}\alpha_1u_1\quad\quad\quad\quad(4-3)x(k+1)=λ1k+1[α1u1+i=2∑n(λ1λi)αiui]≈λ1k+1α1u1(4−3)

因此，可把x(k+1)x^{(k+1)}x(k+1)作为与λ1\lambda_1λ1相应的特征向量的近似。由式：

x(k+1)≈λ1k+1α1u1,x(k)≈λ1kα1u1x^{(k+1)}\approx\lambda_1^{k+1}\alpha_1u_1,\quad x^{(k)}\approx\lambda_1^k\alpha_1u_1x(k+1)≈λ1k+1α1u1,x(k)≈λ1kα1u1

此时我们可以得出：

λ1≈xi(k+1)xi(k)(i=1,2,…,n)(4−4)\lambda_1\approx\frac{x_i^{(k+1)}}{x_i^{(k)}}\quad(i=1,2,\dots,n)\quad\quad\quad\quad\quad\quad(4-4)λ1≈xi(k)xi(k+1)(i=1,2,…,n)(4−4)

其中xi(k)x^{(k)}_ixi(k)为x(k)x^{(k)}x(k)的第i个分量。按式（4-1）和式（4-4）计算矩阵AAA按模最大的特征值与相应的特征向量的方法称为幂法。

需要说明的是，如果x(0)x^{(0)}x(0)的选取恰恰使得α1=0\alpha_1=0α1=0，幂法计算仍能进行。因为计算过程中舍入误差的影响，迭代若干次后，必然会产生一个向量x(k)x^{(k)}x(k)，它在u1u_1u1方向上的分量不为零。这样，以后的计算就满足所设条件。另外需要注意的一点是，因为x(k)≈λ1kα1u1x^{(k)}\approx\lambda_1^k\alpha_1u_1x(k)≈λ1kα1u1，计算过程中可能会出现溢出(∣λ1∣>1|\lambda_1|>1∣λ1∣>1)或成为0(∣λ1∣<1|\lambda_1|<1∣λ1∣<1)的情形。为避免这一情形的出现，实际计算时每次迭代所求的的向量都要归一化。因此，幂法实际使用的计算公式是：

{y(k)=x(k)xr(k),∣xr(k)∣=max⁡1≤i≤n∣xi(k)∣(k=0,1,2,…)x(k+1)=Ay(k)λ1≈xr(k+1)(4−5)\begin{cases}y^{(k)}=\frac{x^{(k)}}{x_r^{(k)}},\ |x_r^{(k)}|=\max\limits_{1\le i\le n}|x_i^{(k)}|\quad(k=0,1,2,\dots)\\x^{(k+1)}=Ay^{(k)}\\\lambda_1\approx x_r^{(k+1)}\\\end{cases}\quad\quad\quad\quad\quad(4-5)⎩⎪⎪⎨⎪⎪⎧y(k)=xr(k)x(k), ∣xr(k)∣=1≤i≤nmax∣xi(k)∣(k=0,1,2,…)x(k+1)=Ay(k)λ1≈xr(k+1)(4−5)

待补充 92

1.2 幂法的加速

因为幂法的收敛条件是线性的，而且依赖于收敛因子∣λ2/λ1∣|\lambda_2/\lambda_1|∣λ2/λ1∣，当收敛因子接近于1时，幂法收敛很慢。幂法的加速有许多方法，下面介绍其中两种。

原点移位法

容易看出，矩阵AAA与A−λ0IA-\lambda_0IA−λ0I的特征值有以下关系：若λi\lambda_iλi是AAA的特征值，则λi−λ0\lambda_i-\lambda_0λi−λ0就是A−λ0IA-\lambda_0IA−λ0I的特征值，而且相应的特征向量不变。如果对矩阵A−λ0IA-\lambda_0IA−λ0I按式（4-1）计算，则有：

xk+1=(A−λ0I)x(k)=(λ1−λ0)k+1[α1u1+((λ2−λ0)(λ1−λ0))α2u2+⋯+((λn−λ0)(λ1−λ0))αnun]x^{k+1}=(A-\lambda_0I)x^{(k)}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\\=(\lambda_1-\lambda_0)^{k+1}[\alpha_1u_1+(\frac{(\lambda_2-\lambda_0)}{(\lambda_1-\lambda_0)})\alpha_2u_2+\cdots+(\frac{(\lambda_n-\lambda_0)}{(\lambda_1-\lambda_0)})\alpha_nu_n]xk+1=(A−λ0I)x(k)=(λ1−λ0)k+1[α1u1+((λ1−λ0)(λ2−λ0))α2u2+⋯+((λ1−λ0)(λn−λ0))αnun]

适当选取λ0\lambda_0λ0，使得∣λ1−λ0∣>∣λi−λ0∣|\lambda_1-\lambda_0|>|\lambda_i-\lambda_0|∣λ1−λ0∣>∣λi−λ0∣，且：

∣λi−λ0λ1−λ0∣<∣λ2λ1∣(i=2,3,…,n)|\frac{\lambda_i-\lambda_0}{\lambda_1-\lambda_0}|<|\frac{\lambda_2}{\lambda_1}|\quad(i=2,3,\dots,n)∣λ1−λ0λi−λ0∣<∣λ1λ2∣(i=2,3,…,n)

这样，用幂法计算A−λ0IA-\lambda_0IA−λ0I的最大模特征值λ1−λ0\lambda_1-\lambda_0λ1−λ0及相应特征向量的收敛速度比对AAA用幂法计算要快。这种加速收敛的方法称为原点移位法。

待补充

幂法的艾特肯（Aitken）加速

如果序列{ak}\{a_k\}{ak}线性收敛到aaa，即：

lim⁡k→∞ak+1−aak−a=c≠0\lim\limits_{k\to\infty}\frac{a_{k+1}-a}{a_k-a}=c\ne0k→∞limak−aak+1−a=c=0

则当kkk充分大时，有

ak+2−aak+1−a≈ak+1−aak+1−a≈ak+1−aak−a\frac{a_{k+2}-a}{a_{k+1}-a}\approx\frac{a_{k+1}-a}{a_{k+1}-a}\approx\frac{a_{k+1}-a}{a_k-a}ak+1−aak+2−a≈ak+1−aak+1−a≈ak−aak+1−a

由此可解出：

a≈ak+2ak−ak+12ak+2−2ak+1+ak=ak−(ak+1−ak)2ak+2−2ak+1+ak=ak^(4−10)a\approx\frac{a_{k+2}a_k-a_{k+1}^2}{a_{k+2}-2a_{k+1}+a_k}=a_k-\frac{(a_{k+1}-a_k)^2}{a_{k+2}-2a_{k+1}+a_k}=\hat{a_k}\quad\quad\quad\quad\quad(4-10)a≈ak+2−2ak+1+akak+2ak−ak+12=ak−ak+2−2ak+1+ak(ak+1−ak)2=ak^(4−10)

序列{ak^}\{\hat{a_k}\}{ak^}比{ak}\{a_k\}{ak}更快地收敛到aaa，这就是Aitken加速法。将这一方法用于幂法所产生的序列{ak}\{a_k\}{ak}，可加快幂法的收敛速度。

待补充

1.3 反幂法

反幂法是计算矩阵按模最小的特征值及特征向量的方法，也是修正特征值、求相应特征向量的最有效的方法。

设 AAA 为 n×nn\times nn×n 阶非奇异矩阵，λ,u\lambda,uλ,u为AAA的特征值与相应的特征向量，则A−1A^{-1}A−1的特征值是AAA的特征值的倒数，而相应的特征向量不变，即：

A−1u=1λuA^{-1}u=\frac{1}{\lambda}uA−1u=λ1u

因此，若对矩阵A−1A^{-1}A−1用幂法，即可计算出A−1A^{-1}A−1的按模最大的特征值，其倒数恰为AAA的按模最小的特征值。这就是反幂法的基本思想。

因为A−1A^{-1}A−1的计算比较麻烦，而且往往不能保持矩阵AAA的一些好性质（如稀疏性），因此，反幂法在实际运算时以求解方程组：

Ax(k+1)=x(k)(4−12)Ax^{(k+1)}=x^{(k)}\quad\quad\quad\quad\quad\quad(4-12)Ax(k+1)=x(k)(4−12)

代替幂法迭代：

x(k+1)=A−1x(k)x^{(k+1)}=A^{-1}x^{(k)}x(k+1)=A−1x(k)

求得x(k+1)x^{(k+1)}x(k+1)，每迭代一次要解一个线性方程组（4-12）。由于矩阵在迭代过程中不变，故可对AAA先进行三角分解。这样，每次迭代只要解两三个三角方程组。

待补充

补充内容