6 矩阵特征值的数值计算

6 矩阵特征值的数值计算

6.1 特征值与特征向量

设\(A\)是\(n\)阶矩阵，\(x\)是非零列向量，如果存在数\(\lambda\)满足

\[Ax=\lambda x \]

那么称\(\lambda\)是矩阵\(A\)的一个特征值，\(x\)则是属于\(\lambda\)的一个特征向量。

理论上要求解矩阵\(A\)的所有特征值需要求解特征多项式\(f(\lambda)=\det(\lambda I-A)\)的所有根，而\(f(\lambda)\)是一个\(n\)阶多项式，因为5次以上代数多项式是没有求根公式的，因此直接求解特征多项式的难度过大。于是，就发展出了一系列的数值算法来计算矩阵的特征值。

格施哥林(Gerschgorin)圆盘定理 设\(A=[a_{ij}]\)是一个\(n\)阶矩阵，则\(A\)的每一个特征值必定属于下述某个圆盘之中：

\[|\lambda-a_{ii}|\le\sum_{k=1\\k\ne i}^n|a_{ik}|\quad(i=1,2,\cdots,n) \]

通过上述圆盘定理，我们可以估计矩阵特征值的大致范围。

6.2 幂法

幂法是通过矩阵特征向量来求出特征值的一种迭代法，主要用来求解矩阵按模最大的特征值（称为主特征值）及其对应的特征向量。

假设n阶矩阵A的特征值为\(\lambda_1,\lambda_2,\cdots,\lambda_n\)满足

\[|\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_n| \]

并且假设A具有n个线性无关的特征向量\(x_1,x_2,\cdots,x_n\)，需要注意的是这里只考虑了主特征值唯一的情况。幂法的目标就是求解\(\lambda_1,x_1\)。

幂法的基本思想是从任意选取的初始向量\(V^{(0)}\)出发，构造如下向量序列

\[V^{(0)}\ne0,V^{(1)}=AV^{(0)},\cdots,V^{(k+1)}=AV^{(k)},\cdots \]

因为A的特征向量构成了\(\mathbb{R}^n\)空间中的一组基，所以\(V^{(0)}\)可以表示为

\[V^{(0)}=a_1x_1+a_2x_2+\cdots+a_nx_n \]

这里假设\(a_1\ne0\)，于是

\[\begin{aligned} V^{(k)}&=AV^{(k-1)}=\cdots=A^{k}V^{(0)}\\ &=a_1A^kx_1+a_2A^kx_2+\cdots+a_nA^kx_n\\ &=a_1\lambda_1^kx_1+a_2\lambda_2^kx_2+\cdots+a_n\lambda_n^kx_n\\ &=\lambda_1^k(a_1x_1+a_2\left(\frac{\lambda_2}{\lambda_1}\right)^kx_2+\cdots+ a_n\left(\frac{\lambda_n}{\lambda_1}\right)^kx_n) \end{aligned} \]

当\(k\)充分大的时候，我们有

\[\left(\frac{\lambda_i}{\lambda_1}\right)^k\approx0,i=2,3,\cdots,n \]

于是，我们得到\(V^{(k+1)}\approx a_1\lambda_1^{k+1}x_1,V^{(k)}\approx a_1\lambda_1^{k}x_1\)。此时这两个向量的第\(j\)个分量，有如下关系：

\[V^{(k+1)}_j\approx{a_1\lambda_1^{k+1}[x_1]_j}=\lambda_1(a_1\lambda^k_1[x_1]_j) \approx \lambda_1V_j^{(k)} \]

所以

\[\lambda_1\approx\frac{V^{(k+1)}_j}{V^{(k)}_j} \]

并且此时\(V^{(k)}\)就可以作为属于特征值\(\lambda_1\)的一个近似特征向量，这是因为特征向量乘以任意常数仍然都是当前所属特征值的一个特征向量。

上述分析已经说明了幂法的本质，但如果直接使用上述方法则会出现严重的计算溢出问题。因为，如果\(|\lambda_1|>1\)，就会有\(\lim_{k\rightarrow\infty}\lambda_1^k=\infty\)。如果\(|\lambda_1|<1\)，则会有\(\lim_{k\rightarrow\infty}\lambda_1^k=0\)。于是在计算过程中会发生上溢或下溢。因此实际应用时，应该采用改进过后的归一化幂法，算法如下：

\[\left\{\begin{aligned} &j=\mathop{\arg\max}_i|V_i^{(k)}|\\ &m_k=V_j^{(k)}\\ &U^{(k)}=V^{(k)}/m_k\\ &V^{(k+1)}=AU^{(k)} \end{aligned}\right.\quad k=0,1,2,\cdots \]

其中\(m_k\)表示第\(k\)次迭代过程中\(V^{(k)}\)的按模最大的分量。当\(|m_k-m_{k-1}|<\varepsilon\)或达到最大迭代上限时，算法停止并返回\(m_k\)作为主特征值，\(U^{(k)}\)作为特征向量。

算法分析：使用归一化幂法，当\(k\)充分大的时候有

\[U^{(k-1)}\approx\frac{a_1\lambda_1^{k-1}}{M_{k-1}}x_1 \]

其中\(M_{k-1}=m_0m_1\cdots m_{k-1}\)。由于\(U^{(k-1)}\)的每一个分量的绝对值都是小于等于1，并且绝对值最大的分量一定等于1。因为

\[V^{(k)}=AU^{(k-1)}\approx\frac{a_1\lambda_1^{k}}{M_{k-1}}x_1\approx\lambda_1U^{(k-1)} \]

所以\(V^{(k)}\)的按模最大的分量就约等于\(\lambda_1\)，这意味着\(m_k\approx\lambda_1\)。

综上可知，当\(k\)充分大的时候，可以将\(m_k\)作为主特征值的近似值，\(U^{(k)}\)则是属于主特征值的近似特征向量。

上述分析对于主特征值是重根的情况也是成立的，例如

\[\lambda_1=\lambda_2=\cdots=\lambda_r \]

并且满足

\[|\lambda_1|=\cdots=|\lambda_r|>|\lambda_{r+1}|\ge\cdots\ge|\lambda_n| \]

那么此时\(x_1,x_2,\cdots,x_r\)均是是属于\(\lambda_1\)的特征向量，只要此时\(x_1,x_2,\cdots,x_n\)仍然构成\(\mathbb{R}^n\)空间中的一组基，那么幂法仍然是可用的。

幂法的几个问题：

1. 在幂法中如果初始向量选取不当，可能会出现\(a_1=0\)的情况，但是考虑到计算过程中的舍入误差的影响，迭代到某一步骤的时候，在\(x_1\)方向的分量就不为零了。因此，实际初始向量的选择可以是任意的。
2. 影响幂法的收敛速度主要是\(\lambda_2/\lambda_1\)的值，如果这个值越小，收敛速度越快。

6.3 反幂法

设矩阵\(A\in\mathbb{R}^{n\times n}\)是非奇异的，\(\lambda\)表示\(A\)的特征值，\(x\)则是相应的特征向量，如果有\(\lambda\ne0\)，那么

\[Ax=\lambda x\Leftrightarrow A^{-1}x=\lambda^{-1}x \]

这说明\(\lambda^{-1}\)是\(A^{-1}\)的特征值，\(x\)则是相应的特征向量。

基于上述这个性质，反幂法常用来求矩阵按模最小的特征值及其对应的特征向量。

假设\(A\)的特征值满足

\[|\lambda_1|\ge|\lambda_2|\ge\cdots\ge|\lambda_{n-1}|>|\lambda_n|>0 \]

那么\(A^{-1}\)的特征值（\(\lambda_1^{-1},\cdots,\lambda_n^{-1}\)）满足

\[\frac{1}{|\lambda_n|}>\frac{1}{|\lambda_{n-1}|}\ge\cdots\ge \frac{1}{|\lambda_1|} \]

通过幂法求出\(A^{-1}\)按模最大的特征值就是\(A\)按模最小的特征值。算法如下：

\[\left\{\begin{aligned} &j=\mathop{\arg\max}_i|V_i^{(k)}|\\ &m_k=V_j^{(k)}\\ &U^{(k)}=V^{(k)}/m_k\\ &V^{(k+1)}=A^{-1}U^{(k)} \end{aligned}\right.\quad k=0,1,2,\cdots \]

当\(|1/m_k-1/m_{k-1}|<\varepsilon\)或达到最大迭代上限时，算法停止并返回\(1/m_k\)作为按模最小特征值，\(U^{(k)}\)作为特征向量。

直接计算\(A\)的逆矩阵非常耗时，通常通过从方程\(AV^{(k+1)}=U^{(k)}\)中求解出\(V^{(k+1)}\)，来代替算法的最后一步。此外，为了进一步节省工作量可以事先将\(A\)进行LU分解。那么之后迭代过程中，就只需要求解两个三角形方程组。当矩阵阶数很大时，将有效减少计算量。

6.4 原点平移法

设矩阵\(A\in\mathbb{R}^{n\times n}\)的特征值为\(\lambda_1,\cdots,\lambda_n\)，\(x_1,\cdots,x_n\)则是相应的特征向量，那么

\[Ax_i=\lambda_ix_i\Leftrightarrow (A-qI)x_i=(\lambda_i-q)x_i \]

其中\(q\)是一个实数。

通过上述式子可知，矩阵\(A-qI\)的特征值为\(\lambda_1-q,\cdots,\lambda_n-q\)，这相当于是\(A\)特征值进行了步长为\(q\)的平移，并且保持特征向量不变。

利用这一点性质，如果我们想求出矩阵\(A\)距离\(q\)最近的特征值，那么就可以通过原点平移法将这个最接近\(q\)的特征值平移到0附近。然后使用反幂法求解即可，具体做法就是：

1. 令\(B=A-qI\)。
2. 使用反幂法求解\(B\)的按模最小特征值\(1/m_k\)，特征向量\(U^{(k)}\)。
3. 返回\(1/m_k+q\)作为特征值，\(U^{(k)}\)作为特征向量。

6.5 矩阵的正交分解

6.5.1 Householder变换

定义 设\(u\in\mathbb{R}^n\)，且\(||u||_2=1\)，则矩阵

\[H=I-2uu^T \]

称为Householder矩阵或者Householder变换，也可以简称为H矩阵或H变换。

定理1 H矩阵具有如下性质：

1. \(H\)是对称矩阵：\(H=H^T\)。
2. \(H\)是正交矩阵：\(H^TH=I,H^T=H^{-1}\)。
3. \(H\)变换保持向量长度不变：\(\forall v\in\mathbb{R}^n,||Hv||_2=||v||_2\)。
4. 设\(S\)是经过原点的超平面，且法向量为\(u\)，那么对于任意非零向量\(v\in\mathbb{R}^n\)，都有\(Hv\)与\(v\)关于超平面\(S\)对称。

证明：

性质1，2可以很容易被验证。

因为

\[||Hv||_2^2=(Hv)^T(Hv)=v^T(H^TH)v=v^Tv=||v||_2^2 \]

即\(||Hv||_2=||v||_2\)，性质3成立。

为了证明性质4，假设\(v=v_{\parallel}+v_{\perp}\)，其中\(v_{\parallel}\)是与超平面平行的分量，那么\(u^Tv_{\parallel}=0\)；\(v_{\perp}\)是与超平面垂直的分量，因此必然存在实数\(\lambda\)，使得\(v_{\perp}=\lambda u\)。

\[\begin{aligned} Hv&=Hv_{\parallel}+Hv_{\perp}\\ &=(v_{\parallel}-2uu^Tv_{\parallel})+(v_{\perp}-2uu^Tv_{\perp})\\ &=v_{\parallel}+(\lambda{u}-2\lambda{uu^Tu})\\ &=v_{\parallel}-\lambda{u}\\ &=v_{\parallel}-v_{\perp} \end{aligned} \]

\(Hv=v_{\parallel}-v_{\perp}\)表明\(Hv\)和\(v\)是关于超平面\(S\)对称的。

定理2 任意非零向量\(v\in\mathbb{R}^n\)，都可以构造出一个合适的H矩阵使得\(Hv=ce_1\)。

其中\(c=-sign(v_1)||v||_2,e_1=[1,0,\cdots,0]^T\)。\(v_1\)表示向量\(v\)的第一个分量。

构造方法：

1. 计算\(w=v-ce_1=v+sign(v_1)||v||_2e_1\)，此时\(w\)就是超平面\(S\)的一个法向量。
2. 令\(u=w/||w||_2\)，便可以计算出\(H\)。

\[H=I-2\frac{ww^T}{||w||_2^2} \]

令\(\beta=2/||w||_2^2\)，则

\[\begin{aligned} \beta=\frac{2}{||w||_2^2}&=\frac{2}{(v-ce_1)^T(v-ce_1)}\\ &=\frac{2}{v^Tv-cv^Te_1-ce_1^Tv+c^2e_1^Te_1}\\ &=\frac{1}{||v||_2(||v||_2+|v_1|)} \end{aligned} \]

于是\(H=I-\beta{ww^T}\)。

定理2的推广 可以找到一个H矩阵，将一个非零向量的某几个连续分量变为零。

假设\(v\in\mathbb{R}^n,c=-sign(v_i)\sqrt{v_i^2+v_{i+1}^2+\cdots+v_j^2},2\le i,j\le n\)。我们的目标就是通过H变换将\(v\)变换为\(Hv=[v_1,\cdots,c,0,\cdots,0,v_{j+1},\cdots,v_n]\)，即将第\(i+1\)到\(j\)个分量变成0。

构造方法：

1. 计算\(w=[0,\cdots,0,v_i-c,v_{i+1},\cdots,v_j,0,\cdots,0]\)。
2. 令 \(u=w/||w||_2\)即可算出\(H\)。

令\(\bar{v}=[v_i,\cdots,v_j]^T\)，与定理2中的分析相同，我们可以得到

\[\beta=\frac{1}{||\bar{v}||_2(||\bar{v}||_2+|v_i|)} \]

于是\(H=I-\beta{ww^T}\)。

6.5.2 Hessenberg矩阵

定义 设矩阵\(A=[a_{ij}]_{n\times{n}}\)，当\(i>j+1\)时\(a_{ij}=0\)，则称矩阵\(A\)为上hessenberg矩阵（或拟上三角形矩阵）；当\(j>i+1\)时\(a_{ij}=0\)，则称矩阵\(A\)为下hessenberg矩阵（或称拟下三角形矩阵）。

以4阶矩阵为例：

\[\left[\begin{matrix} a_{11}&a_{12}&a_{13}&a_{14}\\ a_{21}&a_{22}&a_{23}&a_{24}\\ 0&a_{32}&a_{33}&a_{34}\\ 0&0&a_{43}&a_{44}\\ \end{matrix}\right]\quad \left[\begin{matrix} a_{11}&a_{12}&0&0\\ a_{21}&a_{22}&a_{23}&0\\ a_{31}&a_{32}&a_{33}&a_{34}\\ a_{41}&a_{42}&a_{43}&a_{44}\\ \end{matrix}\right] \]

左侧是上hessenberg矩阵，右侧是下hessenberg矩阵。

定理3 对于任意n阶矩阵\(A\)，总存在正交矩阵\(Q\)，使得\(B=QAQ^{-1}\)为上hessenberg矩阵。即任意n阶矩阵总相似于一个上hessenberg矩阵。

证明：证明的思路就是构造householder矩阵，将任意的n阶矩阵逐步变成上hessenberg矩阵，如果

\[A=\left[\begin{matrix} a_{11}&\cdots&\cdots&\cdots&a_{1n}\\ a_{21}&\cdots&\cdots&\cdots&a_{2n}\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ a_{n1}&\cdots&\cdots&\cdots&a_{nn} \end{matrix}\right] \]

我们的目标就是找到一个\(\bar{H_1}\in\mathbb{R}^{(n-1)\times(n-1)}\)使得

\[\bar{H_1}\left[\begin{matrix} a_{21}\\a_{31}\\\vdots\\a_{n1} \end{matrix}\right]= \left[\begin{matrix} c_1\\0\\\vdots\\0 \end{matrix}\right] \]

显然根据定理2我们可以找到这样的\(\bar{H_1}\)，然后构造一个n阶矩阵

\[H_1=\left[\begin{matrix} 1&0\\ 0&\bar{H_1} \end{matrix}\right] \]

可以验证\(H_1\)是正交且对称的矩阵（\(H_1=H_1^T=H_1^{-1}\)），于是利用矩阵分块乘积可以得到

\[\begin{aligned} H_1AH_1&=\left[\begin{matrix} 1&0\\ 0&\bar{H_1} \end{matrix}\right] \left[\begin{matrix} a_{11}&\cdots&\cdots&\cdots&a_{1n}\\ a_{21}&\cdots&\cdots&\cdots&a_{2n}\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ a_{n1}&\cdots&\cdots&\cdots&a_{nn} \end{matrix}\right] \left[\begin{matrix} 1&0\\ 0&\bar{H_1} \end{matrix}\right]\\ &=\left[\begin{matrix} *&\cdots&\cdots&\cdots&*\\ c_1&\cdots&\cdots&\cdots&*\\ 0&\cdots&\cdots&\cdots&\vdots\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ 0&\cdots&\cdots&\cdots&* \end{matrix}\right]=A_1 \end{aligned} \]

其中\(*\)表示一个数。

因此\(A_1=H_1AH_1\)，即矩阵\(A_1\)相似于\(A\)。

基于这个思想，第\(k\)步变换的矩阵为

\[H_k=\left[\begin{matrix} I_{k\times k}&0\\ 0&\bar{H_k} \end{matrix}\right] \]

一共需要\(n-2\)步就能够得到上hessenberg矩阵，即

\[H_{n-2}\cdots H_2H_1AH_1H_2\cdots H_{n-2}=\left[\begin{matrix} *&\cdots&\cdots&\cdots&*\\ c_1&\cdots&\cdots&\cdots&*\\ 0&c_2&\cdots&\cdots&\vdots\\ \vdots&\ddots&\ddots&\cdots&\vdots\\ 0&\cdots&0&c_{n-1}&* \end{matrix}\right]=A_{n-2} \]

于是\(A_{n-2}=HAH^{-1}\)，其中\(H=H_{n-2}\cdots H_2H_1\)。这说明经过\(n-2\)步变换之后，一定可以得到上hessenberg矩阵，并且\(A_{n-2}\)相似于\(A\)。

因此所求的正交矩阵\(Q=H\)。

6.5.3 QR分解

定理4 对于任意n阶矩阵\(A\)，总存在正交矩阵\(Q\)和上三角阵\(R\)，使得\(A=QR\)。

证明：证明方法与定理3中的类似，根据householder变换，构造出一系列的\(H_1,H_2,\cdots,H_{n-1}\)，使得

\[H_{n-1}\cdots H_2H_1A=R\Leftrightarrow A=(H_1H_2\cdots H_{n-1})R \]

其中\(Q=H_1H_2\cdots H_{n-1}\)，\(R\)就是上三角阵。

需要注意的是，任意的n阶矩阵\(A\)总可以进行QR分解，但是其分解一般是不唯一的。但是如果\(A\)是可逆矩阵，且上三角阵\(R\)的对角线元素均取正值，则\(A\)的QR分解是唯一的。

QR方法 是幂法的一种推广和变形，可以用来求解任意矩阵的特征值，能够有效地解决阶数不高的任意实方阵的全部特征值问题。

因为实矩阵\(A\)总可以进行QR分解，于是\(A=QR\)。并且令\(B=RQ\)，那么

\[B=RQ=(Q^{-1}A)Q=Q^{-1}AQ \]

即\(B\)与\(A\)相似，因此它们有相同的特征值。基于这个想法，可以得到基本的QR算法：

1. 令\(A_1=A\)。
2. 对于\(k=1,2,\cdots\)，作\(A_k\)的QR分解：\(A_k=Q_kR_k\)。
3. \(A_{k+1}=R_kQ_k\)。

可以验证：\(A_{k+1}=Q_k^{-1}A_kQ_k\)，于是\(A_{k+1}\)总是和\(A_k\)相似。因此，在整个迭代过程中总是保持矩阵\(A\)的特征值不变。

实Schur分解定理 对于任意方阵\(A\in\mathbb{R}^{n\times n}\)，QR方法产生的\(\{A_k\}\)序列收敛于形如下式的分块上三角阵：

\[\left[\begin{matrix} R_{11}&R_{12}&\cdots&R_{1n}\\ 0&R_{22}&\cdots&R_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&R_{nn}\\ \end{matrix}\right] \]

其中\(R_{ii}\)为实数或者\(2\times2\)的矩阵块。如果\(R_{ii}\)是一个实数，那么\(R_{ii}\)就是矩阵\(A\)的一个特征值；如果\(R_{ii}\)是一个\(2\times2\)的矩阵块，那么\(R_{ii}\)的特征值就是\(A\)的一对共轭复数特征值。

QR方法需要注意的问题

1. 在实际中直接使用QR方法时，为了加快收敛，可以先将\(A\)转化为上hessenberg矩阵，然后在使用QR分解求特征值。
2. 迭代停止的条件：（1）设置最大迭代上限。（2）设定一个误差上限\(\varepsilon\)，只要某一个范数满足\(||diag(A_{k+1}-A_k)||<\varepsilon\)成立，停止迭代。