8.1 前言#
在许多科学和工程问题中, 例如振动问题和稳定性问题等, 都会遇到特征值和特征向量的计算问题. 所谓矩阵特征值问题是指: 给定 $n$ 阶方阵 $\mathbf{A} \in \mathbb{R}^{n \times n}$ , 寻找常数 $\lambda \in \mathbb{C}$ 和非零向量 $\mathbf{x} \in \mathbb{C}^n$ , 使得
$$ \boldsymbol{A} \boldsymbol{x} = \lambda \boldsymbol{x},\tag{8.1} $$其中 $\mathbb{C}$ 和 $\mathbb{C}^n$ 分别表示复数域和 $n$ 维复向量空间.
对于 n 阶方阵 A, 称
$$ p (\lambda) = \det (\mathbf{A} - \lambda \mathbf{I})\tag{8.2} $$为矩阵 A 的特征多项式, 其中 I 为 n 阶单位矩阵, $\det(A - \lambda I)$ 表示方阵 A - $\lambda I$ 的行列式. 由线性代数知识知, $p(\lambda)$ 是关于 $\lambda$ 的 n 次多项式, 且 p 的零点即为矩阵 A 的全部特征值. 设 $\lambda_{i}$ 为矩阵 A 的特征值, 由于 $\det(A - \lambda_{i}I) = 0$ , 故方程 $(A - \lambda I)x = 0$ 必有非零解 $x^{(i)}$ , 称 $x^{(i)}$ 为矩阵 A 对应于特征值 $\lambda_{i}$ 的特征向量. 上述内容实际上就是线性代数中介绍的 n 阶矩阵 A 的特征值及其对应特征向量的计算方法.
例8.1.1 计算矩阵 $\left( \begin{array}{ccc}1 & 0 & 2\\ 0 & 1 & -1\-1 & 1 & 1 \end{array} \right)$ 的全部特征值.
解: 矩阵 A 的特征多项式为
$$ p (\lambda) = \det (\boldsymbol{A} - \lambda \boldsymbol{I}) = \det \left(\left[ \begin{array}{c c c} 1 - \lambda & 0 & 2 \\ 0 & 1 - \lambda & - 1 \-1 & 1 & 1 - \lambda \end{array} \right]\right) = (1 - \lambda) (\lambda^{2} - 2 \lambda + 4).\tag{8.3} $$矩阵 A 的特征值为 $p(\lambda)=0$ 的解, 即 $\lambda_{1}=1, \lambda_{2,3}=1\pm\sqrt{3}i.$
当矩阵阶数较小时, 可以通过上述方法计算矩阵 A 的全部特征值. 然而, 当矩阵阶数较大时, 特征多项式 (8.2) 的零点没有简单的解析表达式, 只能通过近似计算得到. 而高次方程近似求根的稳定性差, 求得的近似解会有较大误差. 因此对于高阶矩阵的特征值问题, 用上述方法求解特征值问题就不太可取了, 必须构造稳定而有效的数值方法.
下面叙述一些与特征值有关的结论.
若 $\lambda_{i}(i=1,2,\cdots,n)$ 是矩阵 A 的特征值, 则有
(1) $\sum_{i=1}^{n} \lambda_i = \operatorname{tr}(A)$ ,
(2) $\prod_{i=1}^{n}\lambda_{i}=\det(\boldsymbol{A}),$
其中 $\operatorname{tr}(\mathbf{A}) = \sum_{i=1}^{n} a_{ii}$ 表示矩阵 $\mathbf{A}$ 的迹.
设 A, B 都是 n 阶矩阵, 若存在可逆矩阵 P, 使得
$$ \boldsymbol{P}^{- 1} \boldsymbol{A} \boldsymbol{P} = \boldsymbol{B},\tag{8.4} $$则称矩阵 A 与矩阵 B 相似.
定理8.1.5 (Gerschgorin圆盘定理）设 $A = (a_{ij})_{n\times n}$ ，则 $\pmb{A}$ 的每一个特征值必属于下述某个圆盘之中：
$$ | \lambda - a_{ii} | \leqslant \sum_{j = 1, j \neq i}^{n} | a_{ij} |, \quad i = 1, 2, \dots , n,\tag{8.5} $$其中式 (8.5) 表示复平面上以 $a_{ii}$ 为中心, 以 $\sum_{j=1, j \neq i}^{n} |a_{ij}|$ 为半径的圆盘. 证: 设 $\lambda$ 是 $\mathbf{A}$ 的任一特征值, $\mathbf{x}$ 为其对应的特征向量, 即
$$ (\boldsymbol{A} - \lambda \boldsymbol{I}) \boldsymbol{x} = \boldsymbol{0}.\tag{8.6} $$因为 x 是特征向量, 故 $\boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n})^{\mathrm{T}}\neq\mathbf{0}$ .
记 $|x_{i}| = \max_{1\leqslant j\leqslant n}|x_{j}|\neq 0,$ 则上面方程组中的第 $i$ 个方程为
$$ (a_{i 1} x_{1} + a_{i 2} x_{2} + \dots + a_{i, i - 1} x_{i - 1}) - (\lambda - a_{ii}) x_{i} + (a_{i, i + 1} x_{i + 1} + \dots + a_{i, n} x_{n}) = 0,\tag{8.7} $$即
$$ (\lambda - a_{ii}) x_{i} = a_{i 1} x_{1} + a_{i 2} x_{2} + \dots + a_{i, i - 1} x_{i - 1} + a_{i, i + 1} x_{i + 1} + \dots + a_{i, n} x_{n}.\tag{8.8} $$利用三角不等式, 可得
$$ \left| \lambda - a_{ii} \right| \left| x_{i} \right| \leqslant \left| a_{i 1} \right| \left| x_{1} \right| + \left| a_{i 2} \right| \left| x_{2} \right| + \dots + \left| a_{i, i - 1} \right| \left| x_{i - 1} \right| + \left| a_{i, i + 1} \right| \left| x_{i + 1} \right| + \dots + \left| a_{i, n} \right| \left| x_{n} \right|.\tag{8.9} $$因为 $|x_{j}|\leqslant |x_{i}|$ ，所以
$$ | \lambda - a_{ii} | \leqslant \sum_{j = 1, j \neq i}^{n} | a_{ij} | \frac{| x_{j} |}{| x_{i} |} \leqslant \sum_{j = 1, j \neq i}^{n} | a_{ij} |.\tag{8.10} $$设 A 为 n 阶实对称矩阵, $\lambda_{1} \geqslant \lambda_{2} \geqslant \cdots \geqslant \lambda_{n}$ 为其全部特征值, 则
$$ \begin{array}{l} \text{(1)} \lambda_{n} \leqslant \frac{(\boldsymbol{Ax} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})} \leqslant \lambda_{1}, \quad \forall \boldsymbol{x} \in \mathbb{R}^{n}, \boldsymbol{x} \neq \boldsymbol{0}; \\ \text{(2)} \lambda_{n} = \min_{\boldsymbol{x} \neq 0} \frac{(\boldsymbol{Ax} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})}; \\ \text{(3)} \lambda_{1} = \max_{\boldsymbol{x} \neq 0} \frac{(\boldsymbol{Ax} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})}. \end{array} $$证：因为 A 是实对称矩阵，故有完全的特征向量系。设 $\boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \cdots, \boldsymbol{x}^{(n)}$ 是 A 的规范化正交特征向量组，即 $(\boldsymbol{x}^{(i)}, \boldsymbol{x}^{(j)}) = \delta_{ij}$ 。任一非零向量 x 都能表示为 $x = \sum_{i=1}^{n} a_i x^{(i)}$ ，且
$$ \begin{array}{c} (\boldsymbol{x}, \boldsymbol{x}) = \left(\sum_{1 \leqslant i \leqslant n} \alpha_{i} \boldsymbol{x}^{(i)}, \sum_{1 \leqslant i \leqslant n} \alpha_{i} \boldsymbol{x}^{(i)}\right) = \sum_{1 \leqslant i \leqslant n} \alpha_{i}^{2}, \\ (\boldsymbol{Ax}, \boldsymbol{x}) = \left(\sum_{1 \leqslant i \leqslant n} \boldsymbol{A} \alpha_{i} \boldsymbol{x}^{(i)}, \sum_{1 \leqslant i \leqslant n} \alpha_{i} \boldsymbol{x}^{(i)}\right) = \sum_{1 \leqslant i \leqslant n} \lambda_{i} \alpha_{i}^{2}. \end{array}\tag{8.11} $$故
$$ \frac{(\boldsymbol{A} \boldsymbol{x} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})} = \frac{\sum_{1 \leqslant i \leqslant n} \lambda_{i} \alpha_{i}^{2}}{\sum_{1 \leqslant i \leqslant n} \alpha_{i}^{2}} \leqslant \frac{\lambda_{1} \sum_{1 \leqslant i \leqslant n} \alpha_{i}^{2}}{\sum_{1 \leqslant i \leqslant n} \alpha_{i}^{2}} = \lambda_{1}.\tag{8.12} $$同理可证
$$ \frac{(\boldsymbol{A} \boldsymbol{x} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})} \geqslant \lambda_{n}.\tag{8.13} $$在瑞利商中, 分别取 $\boldsymbol{x} = \boldsymbol{x}^{(1)}$ 和 $\boldsymbol{x} = \boldsymbol{x}^{(n)}$ , 可达到瑞利商的最大值 $\lambda_{1}$ 和最小值 $\lambda_{n}$ , 故
$$ \lambda_{1} = \max_{\boldsymbol{x} \neq 0} \frac{(\boldsymbol{Ax} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})}, \quad \lambda_{n} = \min_{\boldsymbol{x} \neq 0} \frac{(\boldsymbol{Ax} , \boldsymbol{x})}{(\boldsymbol{x} , \boldsymbol{x})}.\tag{8.14} $$8.2 幂 方法#
8.2.1 乘幂法#
乘幂法是一种计算实矩阵按模最大的特征值及其对应的特征向量的迭代方法. 对该方法进行简单修改, 就可计算其他特征值. 幂方法的优点之一就是它不但能给出特征值, 还能计算对应的特征向量. 实用中, 幂方法常用来计算已知特征值所对应的特征向量.
为应用乘幂法, 我们需要对实矩阵 A 作如下假设:
(1) 具有 n 个线性无关的特征向量 $\boldsymbol{x}^{(i)}, i = 1, 2, \cdots, n;$
(2) $A\boldsymbol{x}^{(i)} = \lambda_{i}\boldsymbol{x}^{(i)}$ ，其中 $\lambda_{i}$ 是实特征值且满足 $|\lambda_{1}| > |\lambda_{2}| \geqslant \cdots \geqslant |\lambda_{n}|$ .
选定初始向量 $\boldsymbol{v}^{(0)}\in\mathbb{R}^{n}$ . 由于 $\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)},\cdots,\boldsymbol{x}^{(n)}$ 线性无关, 故
$$ \boldsymbol{v}^{(0)} = \sum_{i = 1}^{n} \alpha_{i} \boldsymbol{x}^{(i)} = \alpha_{1} \boldsymbol{x}^{(1)} + \alpha_{2} \boldsymbol{x}^{(2)} + \dots \alpha_{n} \boldsymbol{x}^{(n)}.\tag{8.15} $$利用迭代关系式 $\boldsymbol{v}^{(k+1)} = \boldsymbol{A}\boldsymbol{v}^{(k)}$ 反复计算, 可得向量序列 $\{\boldsymbol{v}^{(k)}\}_{k=0}^{\infty}$ .
下面对向量序列 $\{\pmb{v}^{(k)}\}_{k = 0}^{\infty}$ 进行分析. 因为 $A\pmb{x}^{(i)} = \lambda_i\pmb{x}^{(i)}$ ，所以
$$ \begin{array}{r} \pmb{v}^{(1)} = \pmb{A} \pmb{v}^{(0)} = \lambda_{1} \alpha_{1} \pmb{x}^{(1)} + \lambda_{2} \alpha_{2} \pmb{x}^{(2)} + \dots + \lambda_{n} \alpha_{n} \pmb{x}^{(n)}, \\ \pmb{v}^{(2)} = \pmb{A} \pmb{v}^{(1)} = \lambda_{1}^{2} \alpha_{1} \pmb{x}^{(1)} + \lambda_{2}^{2} \alpha_{2} \pmb{x}^{(2)} + \dots + \lambda_{n}^{2} \alpha_{n} \pmb{x}^{(n)}. \end{array}\tag{8.16} $$一般地, 由归纳法可得
$$ \begin{array}{c} \boldsymbol{v}^{(k)} = \boldsymbol{A} \boldsymbol{v}^{(k - 1)} = \lambda_{1}^{k} \alpha_{1} \boldsymbol{x}^{(1)} + \lambda_{2}^{k} \alpha_{2} \boldsymbol{x}^{(2)} + \dots + \lambda_{n}^{k} \alpha_{n} \boldsymbol{x}^{(n)} \=lambda_{1}^{k} \left(\alpha_{1} \boldsymbol{x}^{(1)} + \left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k} \alpha_{2} \boldsymbol{x}^{(2)} + \dots + \left(\frac{\lambda_{n}}{\lambda_{1}}\right)^{k} \alpha_{n} \boldsymbol{x}^{(n)}\right). \end{array}\tag{8.17} $$当 k 充分大时, 因为 $\left|\frac{\lambda_{i}}{\lambda_{1}}\right| < 1 (i = 2, 3, \cdots, n)$ , 所以
$$ \pmb{v}^{(k)} \approx \lambda_{1}^{k} \alpha_{1} \pmb{x}^{(1)}.\tag{8.18} $$于是
$$ \boldsymbol{v}^{(k + 1)} \approx \lambda_{1}^{k + 1} \alpha_{1} \boldsymbol{x}^{(1)} \approx \lambda_{1} \boldsymbol{v}^{(k)}.\tag{8.19} $$注意到
$$ \boldsymbol{v}^{(k + 1)} = \boldsymbol{A} \boldsymbol{v}^{(k)},\tag{8.20} $$则有
$$ \pmb{A} \pmb{v}^{(k)} \approx \lambda_{1} \pmb{v}^{(k)}.\tag{8.21} $$这表明 $\boldsymbol{v}^{(k)}$ 可以近似地作为矩阵 A 的按模最大特征值 $\lambda_{1}$ 所对应的特征向量. 若将向量 $\boldsymbol{v}^{(k)}$ 的第 j 个分量记为 $\boldsymbol{v}_{j}^{(k)}$ ，则由 (8.20) 及 (8.21) 可知
$$ \lambda_{1} \approx{\frac{{\pmb v}_{j}^{(k + 1)}}{{\pmb v}_{j}^{(k)}}} \quad{\text{或}} \quad \lambda_{1} \approx{\frac{1}{n}} \sum_{j = 1}^{n} {\frac{{\pmb v}_{j}^{(k + 1)}}{{\pmb v}_{j}^{(k)}}}, \qquad{\pmb v}_{j}^{(k)} \neq 0.\tag{8.22} $$根据上述分析, 可构造如下计算矩阵按模最大特征值的乘幂法.
(1) 给定初始非零向量 $v^{(0)}$ 及控制参数 $\varepsilon;$
(2) $\pmb{v}^{(1)} = \pmb{A}\pmb{v}^{(0)},\pmb{v}^{(2)} = \pmb{A}\pmb{v}^{(1)};$
(3) $M_{2}$ 为向量 $\boldsymbol{v}^{(1)}$ 中非零分量所对应的指标集, $N_{2}$ 为向量 $\boldsymbol{v}^{(1)}$ 中非零分量的个数;
(4) $\lambda_{2} = \frac{1}{N_{2}}\sum_{j\in M_{2}}\frac{\boldsymbol{v}_{j}^{(2)}}{\boldsymbol{v}_{i}^{(1)}};$
(5) 对 $k = 2, 3, \cdots$ , 做:
$$ \boldsymbol{v}^{(k + 1)} = \boldsymbol{A} \boldsymbol{v}^{(k)}; $$$M_{k+1}$ 为向量 $\boldsymbol{v}^{(k)}$ 中非零分量所对应的指标集, $N_{k+1}$ 为向量 $\boldsymbol{v}^{(k)}$ 中非零分量的个数;
$$ \lambda_{k + 1} = \frac{1}{N_{k + 1}} \sum_{j \in M_{k + 1}} \frac{\boldsymbol{v}_{j}^{(k + 1)}}{\boldsymbol{v}_{i}^{(k)}}; $$计算 $err_{k + 1} = |\lambda_{k + 1} - \lambda_k|$ ;
当 $err_{k + 1} < \varepsilon$ 时, 停止迭代;
(6) 返回近似按模最大特征值 $\lambda_{k+1}$ 和近似特征向量 $v^{(k+1)}$ .
例8.2.2 用乘幂法计算矩阵 $\mathbf{A} = \left( \begin{array}{ccc}12 & 6 & -6\\ 6 & 16 & 2\-6 & 2 & 16 \end{array} \right)$ 的按模最大特征值及其对应的特征向量.
解: 取初始向量 $(1.0,0.5,-0.5)^{\mathrm{T}}$ 和控制参数 $\varepsilon=1.0\times10^{-4}$ . 利用迭代关系式 $\boldsymbol{v}^{(k+1)}=\boldsymbol{A}\boldsymbol{v}^{(k)}$ 进行计算, 其结果见表 8-1. 故可取按模最大特征值 $\lambda_{1}\approx21.5440$ , 对应的特征向量为
$$ x^{(1)} \approx (36784710912, 2925607232, - 2925607232)^{\mathrm{T}}. $$从上例可以看出, 随着迭代步数 $k$ 的增加, 向量 $\pmb{v}^{(k)}$ 的分量逐渐增大. 事实上, 当 $|\lambda_1| > 1$ 或 $|\lambda_1| < 1$ 时, 向量序列 $\{\pmb{v}^{(k)}\}$ 中, 各个不为零的分量将会随着 $k$ 的增大, 而或无限增大或趋于零. 在计算机上计算时, 就有可能产生所谓“上溢”或“下溢”现象, 从而得到错误的结果.
为了避免这种情况, 常采取如下措施, 即在每次迭代之前, 先将 $\boldsymbol{v}^{(k)}$ 规范化. 规范化的方法为用 $\boldsymbol{v}^{(k)}$ 的按模最大分量, 记为 $\max[\boldsymbol{v}^{(k)}]$ , 去除所有分量, 得到向量 $\boldsymbol{u}^{(k)} = \frac{\boldsymbol{v}^{(k)}}{\max[\boldsymbol{v}^{(k)}]}$ .
例如, 若 $\boldsymbol{v}^{(k)}=(-2,3,-5,1)^{\mathrm{T}}$ , 则 $\max[\boldsymbol{v}^{(k)}]=-5,\boldsymbol{u}^{(k)}=(0.4,-0.6,1,-0.2)^{\mathrm{T}}$ . 利用规范化的方法可将向量 $\boldsymbol{u}^{(k)}$ 的各个分量控制在 $[-1,1]$ 中.
表 8-1 乘幂法的计算结果
| k | $v^{(k)}$ | $\lambda_k$ | $err_k$ |
| 0 | $(1.0, 0.5, -0.5)^{\text{T}}$ | -- | -- |
| 1 | $(18, 13, -13)^{\text{T}}$ | 23.333 3 | -- |
| 2 | $(372, 290, -290)^{\text{T}}$ | 21.760 7 | 1.572 6 |
| 3 | $(7\ 944, 6\ 292, -6\ 292)^{\text{T}}$ | 21.582 6 | 0.178 0 |
| 4 | $(170\ 832, 135\ 752, -135\ 752)^{\text{T}}$ | 21.551 7 | 0.030 9 |
| 5 | $(3\ 679\ 008, 2\ 925\ 520, -2\ 925\ 520)^{\text{T}}$ | 21.545 6 | 0.006 1 |
| 6 | $(79\ 254\ 336, 63\ 031\ 328, -63\ 031\ 328)^{\text{T}}$ | 21.544 3 | 0.001 3 |
| 7 | $(1\ 707\ 427\ 968, 135\ 796\ 408, -1\ 357\ 964\ 608)^{\text{T}}$ | 21.544 1 | 0.000 3 |
| 8 | $(36\ 784\ 710\ 912, 2\ 925\ 607\ 232, -2\ 925\ 607\ 232)^{\text{T}}$ | 21.544 0 | <0.000 1 |
下面再次利用迭代方法并结合规范化方法给出一个向量序列．选定初始向量 $\boldsymbol{v}^{(0)}$ ，利用关系式 $\boldsymbol{u}^{(k)} = \frac{\boldsymbol{v}^{(k)}}{\max[\boldsymbol{v}^{(k)}]}$ 和 $\boldsymbol{v}^{(k+1)} = \boldsymbol{A}\boldsymbol{u}^{(k)}$ 反复迭代计算，可得向量序列 $\{\boldsymbol{u}^{(k)}\}_{k=0}^{\infty}$ 和 $\{\boldsymbol{v}^{(k)}\}_{k=0}^{\infty}$ . 由于
$$ \pmb{v}^{(0)} = \alpha_{1} \pmb{x}^{(1)} + \alpha_{2} \pmb{x}^{(2)} + \dots + \alpha_{n} \pmb{x}^{(n)},\tag{8.23} $$故有
$$ \begin{array}{l} \boldsymbol{u}^{(0)} = \frac{\boldsymbol{v}^{(0)}}{\max [ \boldsymbol{v}^{(0)} ]} = \frac{\alpha_{1} \boldsymbol{x}^{(1)} + \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \alpha_{n} \boldsymbol{x}^{(n)}}{\max [ \alpha_{1} \boldsymbol{x}^{(1)} + \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \alpha_{n} \boldsymbol{x}^{(n)} ]}, \\ \boldsymbol{v}^{(1)} = A u^{(0)} = \frac{\lambda_{1} \alpha_{1} \boldsymbol{x}^{(1)} + \lambda_{2} \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \lambda_{n} \alpha_{n} \boldsymbol{x}^{(n)}}{\max [ \alpha_{1} \boldsymbol{x}^{(1)} + \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \alpha_{n} \boldsymbol{x}^{(n)} ]}. \end{array}\tag{8.24} $$一般地，
$$ \begin{array}{r} \pmb{u}^{(k)} = \frac{\lambda_{1}^{k} \alpha_{1} \pmb{x}^{(1)} + \lambda_{2}^{k} \alpha_{2} \pmb{x}^{(2)} + \cdots + \lambda_{n}^{k} \alpha_{n} \pmb{x}^{(n)}}{\max [ \lambda_{1}^{k} \alpha_{1} \pmb{x}^{(1)} + \lambda_{2}^{k} \alpha_{2} \pmb{x}^{(2)} + \cdots + \lambda_{n}^{k} \alpha_{n} \pmb{x}^{(n)} ]}, \\ \pmb{v}^{(k + 1)} = \pmb{Au}^{(k)} = \frac{\lambda_{1}^{k + 1} \alpha_{1} \pmb{x}^{(1)} + \lambda_{2}^{k + 1} \alpha_{2} \pmb{x}^{(2)} + \cdots + \lambda_{n}^{k + 1} \alpha_{n} \pmb{x}^{(n)}}{\max [ \lambda_{1}^{k} \alpha_{1} \pmb{x}^{(1)} + \lambda_{2}^{k} \alpha_{2} \pmb{x}^{(2)} + \cdots + \lambda_{n}^{k} \alpha_{n} \pmb{x^{(n)}} ]}. \end{array}\tag{8.25} $$简单计算可得,
$$ \begin{array}{l} \boldsymbol{u}^{(k)} = \frac{\alpha_{1} \boldsymbol{x}^{(1)} + \left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k} \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \left(\frac{\lambda_{n}}{\lambda_{1}}\right)^{k} \alpha_{n} \boldsymbol{x}^{(n)}}{\max [ \alpha_{1} \boldsymbol{x}^{(1)} + \left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k} \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \left(\frac{\lambda_{n}}{\lambda_{1}}\right)^{k} \alpha_{n} \boldsymbol{x}^{(n)} ]}, \\ \max [ \boldsymbol{v}^{(k)} ] = \lambda_{1} \frac{\max [ \alpha_{1} \boldsymbol{x}^{(1)} + \left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k + 1} \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \left(\frac{\lambda_{n}}{\lambda_{1}}\right)^{k + 1} \alpha_{n} \boldsymbol{x}^{(n)} ]}{\max [ \alpha_{1} \boldsymbol{x}^{(1)} + \left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k} \alpha_{2} \boldsymbol{x}^{(2)} + \cdots + \left(\frac{\lambda_{n}}{\lambda_{1}}\right)^{k} \alpha_{n} \boldsymbol{x}^{\left(n\right)} ]}. \end{array}\tag{8.26} $$当 k 充分大时, 因为 $\left|\frac{\lambda_{i}}{\lambda_{1}}\right| < 1 (i = 2, 3, \cdots, n)$ , 所以
$$ \pmb{u}^{(k + 1)} \approx \pmb{u}^{(k)} \approx \frac{\pmb{x}^{(1)}}{\max [ \pmb{x}^{(1)} ]}, \qquad \max [ \pmb{v}^{(k)} ] \approx \lambda_{1}.\tag{8.27} $$(1) 给定初始非零向量 $\boldsymbol{v}^{(0)}$ 及控制参数 $\varepsilon;$
(2) $\lambda_{0}$ 为向量 $\boldsymbol{v}^{(0)}$ 的按模最大分量；
(3) $\boldsymbol{u}^{(0)} = \boldsymbol{v}^{(0)} / \lambda_{0};$
(4) 对 $k = 1, 2, \cdots$ ，做：
$\boldsymbol{v}^{(k+1)} = \boldsymbol{A}\boldsymbol{u}^{(k)};$ $\lambda_{k+1}$ 为向量 $\boldsymbol{v}^{(k+1)}$ 的按模最大分量；
$\boldsymbol{u}^{(k+1)}=\boldsymbol{v}^{(k+1)}/\lambda_{k+1};$
计算 $err_{k+1}=|\lambda_{k+1}-\lambda_{k}|;$
当 $err_{k+1}<\varepsilon$ 时，停止迭代；从而
$$ \boldsymbol{u}^{(k + 1)} = \frac{\boldsymbol{v}^{(k + 1)}}{\max [ \boldsymbol{v}^{(k)} ]} \approx \frac{\boldsymbol{Au}^{(k)}}{\lambda_{1}} \approx \frac{\boldsymbol{Au}^{(k + 1)}}{\lambda_{1}}.\tag{8.28} $$于是, $\max[\boldsymbol{v}^{(k)}]$ 可作为按模最大特征值 $\lambda_{1}$ 的近似, $\boldsymbol{u}^{(k+1)}$ 可作为对应的近似特征向量.
改进乘幂法的算法描述如下.
算法 8.2.3(改进乘幂法)
(5) 返回近似按模最大特征值 $\lambda_{k+1}$ 和近似特征向量 $\boldsymbol{u}^{(k+1)}$ .
以下是改进乘幂法的 Matlab 程序.
function [t,y] = eigIPower(a,xinit,ep)
v0 = xinit;
[tv,ti] = max(abs(v0));
lam0 = v0(ti);
u0 = v0/lam0;
flag = 0;
while (flag==0)
v1 = a*u0;
[tv,ti] = max(abs(v1));
lam1 = v1(ti);
u0 = v1/lam1;
err = abs(lam0-lam1);
if (err<=ep)
flag = 1;
end
lam0 = lam1;
end
t = lam1;
y = u0;例 8.2.4 用改进乘幂法计算例 8.2.2 中矩阵 A 的按模最大特征值及其对应的特征向量. 解: 取初始向量 $\boldsymbol{v}^{(0)} = (1.0, 0.5, -0.5)^{\mathrm{T}}$ 和控制参数 $\varepsilon = 1.0 \times 10^{-4}$ . 则 $\lambda_{0} = 1$ 和 $\boldsymbol{u}^{(0)} = \boldsymbol{v}^{(0)}$ . 其他计算结果见表 8-2. 当 k = 9 时 $err_{9} < \varepsilon$ ，迭代停止. 矩阵 A 的按模最大特征值 $\lambda_{1} \approx$
21.544 0, 对应的特征向量 $\boldsymbol{x}^{(1)} \approx (1.000\ 0, 0.795\ 3, -0.795\ 3)^{\mathrm{T}}$ .
表 8-2 改进乘幂法的计算结果
| k | $v^{(k)}$ | $\lambda_k$ | $u^{(k)}$ | $err_k$ |
| 0 | $(1.0,0.5,-0.5)^{\text{T}}$ | 1.000 0 | $(1.0000,0.5000,-0.5000)^{\text{T}}$ | - |
| 1 | $(18,13,-13)^{\text{T}}$ | 18.000 0 | $(1.0000,0.7222,-0.7222)^{\text{T}}$ | 17 |
| 2 | $(20.6667,16.1111,-16.1111)^{\text{T}}$ | 20.666 7 | $(1.0000,0.7796,-0.7796)^{\text{T}}$ | 2.666 7 |
| 3 | $(21.3548,16.9140,-16.9140)^{\text{T}}$ | 21.354 8 | $(1.0000,0.7920,-0.7920)^{\text{T}}$ | 0.688 2 |
| 4 | $(21.5045,17.0886,-17.0886)^{\text{T}}$ | 21.504 5 | $(1.0000,0.7947,-0.7947)^{\text{T}}$ | 0.149 7 |
| 5 | $(21.5358,17.1251,-17.1251)^{\text{T}}$ | 21.535 8 | $(1.0000,0.7952,-0.7952)^{\text{T}}$ | 0.031 3 |
| 6 | $(21.5423,17.1327,-17.1327)^{\text{T}}$ | 21.542 3 | $(1.0000,0.7953,-0.7953)^{\text{T}}$ | 0.006 5 |
| 7 | $(21.5437,17.1343,-17.1343)^{\text{T}}$ | 21.543 7 | $(1.0000,0.7953,-0.7953)^{\text{T}}$ | 0.001 3 |
| 8 | $(21.5439,17.1346,-17.1346)^{\text{T}}$ | 21.543 9 | $(1.0000,0.7953,-0.7953)^{\text{T}}$ | 0.000 3 |
| 9 | $(21.5440,17.1347,-17.1347)^{\text{T}}$ | 21.544 0 | $(1.0000,0.7953,-0.7953)^{\text{T}}$ | <0.000 1 |
从上面的例子可以看出, 乘幂法和改进乘幂法是一类通过矩阵的特征向量来求特征值的迭代法. 它的优点是算法比较简单, 容易在计算机上实现, 计算量也较小, 对于高阶稀疏矩阵较为合适. 其缺点在于收敛速度依赖于比值 $\left|\frac{\lambda_{2}}{\lambda_{1}}\right|$ , 比值越小收敛越快, 当比值接近于 1 时, 收敛会变得很慢. 此时可以使用加速收敛的方法, 如瑞利商加速法等.
8.2.2 反幂法#
反幂法是一种计算矩阵按模最小特征值及其对应特征向量的迭代方法.
为应用反幂法, 我们需要对实矩阵 A 作如下假设:
(1) 具有 n 个线性无关的特征向量 $\boldsymbol{x}^{(i)}$ ;
(2) $\boldsymbol{A}\boldsymbol{x}^{(i)} = \lambda_{i}\boldsymbol{x}^{(i)}$ ，其中 $\lambda_{i}$ 是实特征值且满足 $|\lambda_{1}| \geqslant |\lambda_{2}| \geqslant \cdots > |\lambda_{n}| > 0$ .
由于 $\lambda_{i} \neq 0$ , 故 $\mathbf{A}$ 是非奇异矩阵. 设 $\lambda$ 是非奇异矩阵 $\mathbf{A}$ 的特征值, $\mathbf{x}$ 是 $\lambda$ 对应的特征向量. 容易证明: $\lambda^{-1}$ 是矩阵 $\mathbf{A}^{-1}$ 的特征值, $\mathbf{x}$ 是 $\lambda^{-1}$ 所对应的特征向量. 因此, 由算法假设可知, $\mathbf{A}^{-1}$ 的特征值为 $\{\lambda_{i}^{-1}\}_{i=1}^{n}$ , 且满足
$$ | \lambda_{1}^{- 1} | \leqslant | \lambda_{2}^{- 1} | \leqslant \dots < | \lambda_{n}^{- 1} |.\tag{8.29} $$其对应的特征向量仍为 $\{\pmb{x}^{(i)}\}_{i = 1}^n$ 。这表明非奇异矩阵的按模最小特征值问题等价于其逆阵的按模最大特征值问题。于是，可利用乘幂法求得矩阵 $A^{-1}$ 的按模最大特征值 $\mu_{n}$ 及其对应的特征向量 $\pmb{x}^{(n)}$ 。而矩阵 $A$ 的按模最小特征值 $\lambda_{n} = \mu_{n}^{-1}$ ，对应的特征向量为 $\pmb{x}^{(n)}$ 。
反幂法的算法描述如下:
算法 8.2.5 (反幂法)
(1) 给定初始非零向量 $v^{(0)}$ 及控制参数 $\varepsilon;$
(2) $\lambda_{0}$ 为向量 $v^{(0)}$ 的按模最大分量;
(3) $\pmb{u}^{(0)} = \pmb{v}^{(0)} / \lambda_0$ ;
(4) 对 $k = 1, 2, \cdots$ , 做:
$\boldsymbol{A}\boldsymbol{v}^{(k + 1)} = \boldsymbol{u}^{(k)};$
$\lambda_{k+1}$ 为向量 $\boldsymbol{v}^{(k+1)}$ 的按模最大分量;
$$ \boldsymbol{u}^{(k + 1)} = \boldsymbol{v}^{(k + 1)} / \lambda_{k + 1}; $$计算 $err_{k + 1} = |\lambda_{k + 1}^{-1} - \lambda_k^{-1}|$ ;
当 $err_{k + 1} < \varepsilon$ 时，停止迭代；
(5) 返回近似按模最小特征值 $\lambda_{k+1}^{-1}$ 和近似特征向量 $\boldsymbol{u}^{(k+1)}$ .
在上述算法中, 为避免矩阵求逆, $\boldsymbol{v}^{(k+1)} = \boldsymbol{A}^{-1}\boldsymbol{u}^{(k)}$ 改为等价形式 $\boldsymbol{A}\boldsymbol{v}^{(k+1)} = \boldsymbol{u}^{(k)}$ . 对于该线性方程组, 如果系数矩阵 $\boldsymbol{A}$ 满足三角分解的条件, 则可以采用三角分解法求解. 由于矩阵三角分解仅需计算一次, 所以计算量可大为减少. 类似于乘幂法, 反幂法的收敛速度也依赖于比值 $\left|\frac{\lambda_n}{\lambda_{n-1}}\right|$ 的大小, 比值越小, 收敛越快.
例 8.2.6 用反幂法计算例 8.2.2 中矩阵 A 的按模最小特征值及其对应的特征向量.
解：取初始向量 $\boldsymbol{v}^{(0)}=(1.0,-0.5,0.5)^{\mathrm{T}}$ 和控制参数 $\varepsilon=1.0\times10^{-4}$ . 则 $\lambda_{0}=1$ 和 $\boldsymbol{u}^{(0)}=\boldsymbol{v}^{(0)}$ . 其他计算结果如表 8-3.
表 8-3 反幂法的计算结果
| k | $v^{(k)}$ | $\lambda_{k}^{-1}$ | $u^{(k)}$ | $err_{k}$ |
| 0 | $(1.0, -0.5, 0.5)^{\text{T}}$ | 1 | $(1, -0.5, 0.5)^{\text{T}}$ | -- |
| 1 | $(0.2083, -0.1250, 0.1250)^{\text{T}}$ | 4.8000 | $(1.0000, -0.6000, 0.6000)^{\text{T}}$ | 3.8000 |
| 2 | $(0.2208, -0.1375, 0.1375)^{\text{T}}$ | 4.5283 | $(1.0000, -0.6226, 0.6226)^{\text{T}}$ | 0.2717 |
| 3 | $(0.2237, -0.1403, 0.1403)^{\text{T}}$ | 4.4710 | $(1.0000, -0.6274, 0.6274)^{\text{T}}$ | 0.0573 |
| 4 | $(0.2243, -0.1409, 0.1409)^{\text{T}}$ | 4.4591 | $(1.0000, -0.6284, 0.6284)^{\text{T}}$ | 0.0119 |
| 5 | $(0.2244, -0.1411, 0.1411)^{\text{T}}$ | 4.4566 | $(1.0000, -0.6286, 0.6286)^{\text{T}}$ | 0.0025 |
| 6 | $(0.2244, -0.1411, 0.1411)^{\text{T}}$ | 4.4561 | $(1.0000, -0.6287, 0.6287)^{\text{T}}$ | 0.0005 |
| 7 | $(0.2244, -0.1411, 0.1411)^{\text{T}}$ | 4.4560 | $(1.0000, -0.6287, 0.6287)^{\text{T}}$ | 0.0001 |
| 8 | $(0.2244, -0.1411, 0.1411)^{\text{T}}$ | 4.4560 | $(1.0000, -0.6287, 0.6287)^{\text{T}}$ | <0.0001 |
因此, 矩阵 A 的按模最小特征值约为 4.4561, 对应的特征向量约为 $(1.000\ 0,-0.628\ 7,0.628\ 7)^{\mathrm{T}}$ .
8.2.3 结合原点平移的反幂法#
若已知矩阵 A 的某特征值的粗略近似值 p，则可利用结合原点平移技巧的反幂法来计算该矩阵与数 p 最接近的特征值及其对应的特征向量.
为应用结合原点平移的反幂法, 我们假设
数 p 是矩阵 A 的第 i 个特征值 $\lambda_{i}$ 的近似, 且满足 $0 < |\lambda_{i} - p| < |\lambda_{j} - p|$ , $j \neq i$ .
由算法假设可知, 矩阵 B = A - pI 满足反幂法所要求的假设, 且 $\mu_{i} = \lambda_{i} - p$ 是矩阵 B 的按模最小特征值. 因此, 可以对 B 应用反幂法求出 $\mu_{i}$ 及其对应特征向量 $\boldsymbol{x}^{(i)}$ 的近似值, 而矩阵 A 与数 p 最接近的特征值为 $\lambda_{i} = \mu_{i} + p$ , 其对应的特征向量为 $\boldsymbol{x}^{(i)}$ .
结合原点平移的反幂法的算法描述如下:
算法 8.2.7(带原点平移的反幂法)
(1) 给定初始非零向量 $\boldsymbol{v}^{(0)}$ ，控制参数 $\varepsilon$ 以及数 p;
(2) $\lambda_{0}$ 为 $v^{(0)}$ 中按模最大的那个分量;
(3) $\pmb{u}^{(0)} = \pmb{v}^{(0)} / \lambda_0$ ;
(4) 对 $k = 1, 2, \cdots$ ，做:
$$ (\boldsymbol{A} - p \boldsymbol{I}) \boldsymbol{v}^{(k + 1)} = \boldsymbol{u}^{(k)}; $$$\lambda_{k+1}$ 为向量 $\boldsymbol{v}^{(k+1)}$ 的按模最大分量;
$$ \boldsymbol{u}^{(k + 1)} = \boldsymbol{v}^{(k + 1)} / \lambda_{k + 1}; $$计算 $err_{k + 1} = |\lambda_{k + 1}^{-1} - \lambda_k^{-1}|$ ;
当 $err_{k+1} < \varepsilon$ 时, 停止迭代;
(5) 返回近似特征值 $\lambda_{k+1}^{-1} + p$ 和近似特征向量 $\boldsymbol{u}^{(k+1)}$ .
例 8.2.8 求矩阵 $A = \begin{pmatrix} -12 & 3 & 3 \\ 3 & 1 & -2 \\ 3 & -2 & 7 \end{pmatrix}$ 与数 p = -13 最接近的特征值及其对应的特征向量. 解：对矩阵 B = A - pI 进行 LU 分解, 有
$$ \boldsymbol{B} = \left( \begin{array}{c c c} 1 & 3 & 3 \\ 3 & 14 & - 2 \\ 3 & - 2 & 20 \end{array} \right) = \left( \begin{array}{c c c} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 3 & - \frac{11}{5} & 1 \end{array} \right) \left( \begin{array}{c c c} 1 & 3 & 3 \\ 0 & 5 & - 11 \\ 0 & 0 & - \frac{66}{5} \end{array} \right).\tag{8.30} $$取初始向量 $\boldsymbol{v}^{(0)} = \boldsymbol{u}^{(0)} = (1,1,1)^{\mathrm{T}}$ ，利用结合原点平移的反幂法进行计算。计算结果如表 8-4。因此矩阵 A 与 p = -13 最接近的特征值约为 -13.2202，对应的特征向量约为 $(1.000\ 0, -0.235\ 1, -0.171\ 6)^{\mathrm{T}}$ 。
表 8-4 带原点平移的反幂法的计算结果
| k | $\boldsymbol{v}^{(k)}$ | $\boldsymbol{u}^{(k)}$ | $\lambda_{k+1}^{-1} + p$ | $err_k$ |
| 0 | $(1, 1, 1)^{\text{T}}$ | $(1, 1, 1)^{\text{T}}$ | -12 | - |
| 1 | $(-2.454\ 5, 0.666\ 7, 0.484\ 8)^{\text{T}}$ | $(1.000\ 0, -0.271\ 6, -0.197\ 5)^{\text{T}}$ | -13.407 4 | 1.407 4 |
| 2 | $(-4.597\ 1, 1.078\ 2, 0.787\ 5)^{\text{T}}$ | $(1.000\ 0, -0.234\ 5, -0.171\ 3)^{\text{T}}$ | -13.217 5 | 0.189 9 |
| 3 | $(-4.540\ 9, 1.067\ 6, 0.779\ 3)^{\text{T}}$ | $(1.000\ 0, -0.235\ 1, -0.171\ 6)^{\text{T}}$ | -13.220 2 | 0.002 7 |
| 4 | $(-4.541\ 8, 1.067\ 8, 0.779\ 5)^{\text{T}}$ | $(1.000\ 0, -0.235\ 1, -0.171\ 6)^{\text{T}}$ | -13.220 2 | <0.000 1 |
上例表明, 若已知矩阵 A 的某近似特征值, 则可以利用结合原点平移技巧的反幂法计算其对应的特征向量.
8.3 QR 方法#
本节主要介绍计算矩阵全部特征值的 QR 方法. 该方法以矩阵的 QR 分解为基础, 是计算中小型矩阵全部特征值的最有效方法之一, 且具有收敛速度快、算法稳定等特点.
首先给出 QR 方法的算法描述.
算法 8.3.1(QR 方法)
(1) 令 $A_{1}=A$ ，并给定控制参数 $\varepsilon;$
(2) 对 $k = 1, 2, \cdots$ , 做:
计算矩阵 $A_{k}$ 的 QR 分解 $A_{k}=Q_{k}R_{k};$
计算 $A_{k+1}=R_{k}Q_{k};$
计算 $err_{k + 1} = \| \mathrm{diag}(\pmb{A}_{k + 1} - \pmb{A}_k)\|$
当 $err_{k+1} < \varepsilon$ 时, 停止迭代;
(3) 返回 $A_{k + 1}$
这里 $\mathrm{diag}(A)$ 表示由矩阵 A 对角线元素构成的列向量.
给定矩阵 $A_{1}=A$ ，由递推关系式 $A_{k}=Q_{k}R_{k}$ 和 $A_{k+1}=R_{k}Q_{k}$ 可得矩阵序列 $\{A_{k}\}_{k=1}^{\infty}$ . 下面对该矩阵序列进行分析. 由于
$$ \boldsymbol{A}_{k + 1} = \boldsymbol{Q}_{k}^{- 1} \boldsymbol{A}_{k} \boldsymbol{Q}_{k} = \boldsymbol{Q}_{k}^{- 1} \boldsymbol{Q}_{k - 1}^{- 1} \boldsymbol{A}_{k - 1} \boldsymbol{Q}_{k - 1} \boldsymbol{Q}_{k} = \dots = \boldsymbol{Q}_{k}^{- 1} \boldsymbol{Q}_{k - 1}^{- 1} \dots \boldsymbol{Q}_{1}^{- 1} \boldsymbol{A} \boldsymbol{Q}_{1} \dots \boldsymbol{Q}_{k - 1} \boldsymbol{Q}_{k}.\tag{8.31} $$所以矩阵序列中的每个矩阵 $A_{k}$ 都与矩阵 A 相似. 由于相似矩阵具有相同的特征值, 所以 $A_{k}$ 和 A 具有相同的特征值.
另一方面, 可以证明矩阵序列 $\{\pmb{A}_k\}$ 本质收敛于上三角矩阵或块上三角矩阵, 且对角块为 $1 \times 1$ 矩阵或 $2 \times 2$ 矩阵, 其 $1 \times 1$ 矩阵就是 $\pmb{A}$ 的特征值, 每个 $2 \times 2$ 矩阵含有 $\pmb{A}$ 的一对复特征值. 这种对角块为 $1 \times 1$ 矩阵或 $2 \times 2$ 矩阵的分块上三角矩阵称为拟上三角矩阵. 所谓本质收敛, 即对角块和三角部分是收敛的, 其他元素不一定收敛.
因此, 当 k 充分大时, 由近似上三角阵或近似块上三角阵 $A_{k}$ 可得到矩阵 A 的全部特征值. 若还需计算对应的特征向量, 则可以通过结合原点平移技巧的反幂法来求得.
例 8.3.2 用 QR 方法计算下面矩阵的全部特征值:
$$ \boldsymbol{A} = \left( \begin{array}{c c c c} 5 & - 2 & - 5 & - 1 \\ 1 & 0 & - 3 & 2 \\ 0 & 2 & 2 & - 3 \\ 0 & 0 & 1 & - 2 \end{array} \right). $$解: 令 $A_{1}=A$ , 利用 QR 分解算法对矩阵 $A_{1}$ 进行 QR 分解.
$$ \begin{array}{l} \boldsymbol{A}_{1} = \boldsymbol{Q}_{1} \boldsymbol{R}_{1} \=left( \begin{array}{c c c c} 0.9806 & - 0.0377 & - 0.1923 & - 0.1038 \\ 0.1961 & 0.1887 & - 0.8804 & - 0.4192 \\ 0.0000 & 0.9813 & 0.1761 & 0.0740 \\ 0.0000 & 0.0000 & 0.3962 & - 0.8989 \end{array} \right) \\ \times \left( \begin{array}{c c c c} 5.0992 & - 1.9612 & - 5.4912 & - 0.3922 \\ 0.0000 & 2.0381 & 1.5852 & - 2.5288 \\ 0.0000 & 0.0000 & 2.5242 & - 3.2736 \\ 0.0000 & 0.0000 & 0.0000 & \text{} \quad \text{} \quad \text{} \quad \text{} \quad \text{} \end{array} \right). \end{array} $$然后, 求得 $A_{2} = R_{1}Q_{1}$ ，再对 $A_{2}$ 进行 QR 分解. 一直进行下去, 得到
$$ \boldsymbol{A}_{12} = \left( \begin{array}{c c c c} 4.0000 & * & * & * \\ & 1.8789 & - 3.5910 & * \\ & 1.3290 & 0.1211 & * \\ & & & - 1.0000 \end{array} \right). $$所以, A 的两个特征值为 4.0000 和 -1.0000. 其他两个特征值是方程
$$ \left| \begin{array}{c c} 1.8789 - \lambda & - 3.5910 \\ 1.3290 & 0.1211 - \lambda \end{array} \right| = 0 $$的根, 即
$$ \lambda^{2} - 2 \lambda - 1.8789 \times 0.1211 - 1.3290 \times (- 3.5910) = 0. $$计算可得
$$ \lambda \approx 1 \pm 2 \mathrm{i}. $$事实上, 矩阵 A 的四个特征值分别为 $-1, 4, 1 \pm 2i$ .
评 注#
本章主要介绍几种常用的矩阵特征值问题的数值解法。乘幂法是一种计算矩阵按模最大特征值的方法，该方法简单有效并能给出对应的特征向量。反幂法可用来计算非奇异矩阵的按模最小特征值及其对应的特征向量。此类方法统称为幂方法。幂方法的收敛速度依赖于最大与次大（或最小与次小）特征值的比值，因此有时收敛速度比较慢。为了提高收敛速度，可采用加速方法，如原点平移法和瑞利商加速法等。QR方法是计算矩阵全部特征值的数值方法，该方法利用QR分解将矩阵变为与其相似的上三角阵或拟上三角阵，然后利用相似矩阵具有相同特征值这一结论得到矩阵全部特征值。
对于特征值问题, Matlab 软件提供了 eig 函数. 其格式如下:
$$
[ \mathrm{V}, \mathrm{D} ] = \operatorname{eig} (\mathrm{A}) $$
其中对角阵 D 表示矩阵 A 的全部特征值, 矩阵 V 的第 i 列为 D 中第 i 列对角线元素所对应的特征向量.
习题八#
- 用乘幂法求下列矩阵的按模最大特征值及其对应的特征向量
(1) $\left( \begin{array}{ccc}2 & 3 & 2\\ 10 & 3 & 4\\ 3 & 6 & 1 \end{array} \right),$ (2) $\left( \begin{array}{ccc}3 & -4 & 3\-4 & 6 & 3\\ 3 & 3 & 1 \end{array} \right),$
当特征值有三位小数稳定时迭代终止.
用反幂法求矩阵 $A = \begin{pmatrix} 2 & 0 & 0 \\ 2 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}$ 的按模最小特征值及其对应的特征向量, 当特征值有三位小数稳定时迭代终止.
用幂方法求矩阵 $\mathbf{A} = \left( \begin{array}{ccc}6 & 2 & 1\\ 2 & 3 & 1\\ 1 & 1 & 1 \end{array} \right)$ 的最接近于6的特征值及其对应的特征向量.
已知矩阵 $A = \begin{pmatrix} 2 & 1 & 0 \\ 1 & 3 & 1 \\ 0 & 1 & 4 \end{pmatrix}$ 的近似特征值 $\overline{\lambda}_{3} = 1.2679$ (准确特征值为 $\lambda_{3} = 3 - \sqrt[3]{3}$ ), 试求该特征值对应的特征向量.
用 QR 方法计算下列矩阵的全部特征值
(1)
$$ \left( \begin{array}{c c c} 3 & 1 & 0 \\ 1 & 4 & 2 \\ 0 & 2 & 3 \end{array} \right), $$(2)
$$ \left( \begin{array}{c c c} 1 & 4 & 5 \\ 2 & 5 & 6 \\ 2 & 2 & 0 \end{array} \right). $$数值实验八#
- 考虑 $n$ 阶的三对角矩阵
分别计算当 $n=2,4,6,8,\cdots$ 时, 矩阵 A 的条件数 $\mathrm{cond}_{2}(A)=\frac{\lambda_{\max}(A)}{\lambda_{\min}(A)}$ . 2. 试求 n = 10 阶矩阵
$$ \boldsymbol{A} = \left(a_{ij}\right) = \left(\frac{1}{i + j - 1}\right), \quad i, j = 1, 2, \dots , n. $$的全部特征值及其对应的特征向量.
- 考虑例 8.2.2, 以向量序列 $\{\pmb{u}^{(k)}\}$ 的瑞利商作为近似特征值. 将计算结果与改进乘幂法得到的近似解进行比较, 并解释原因.