线性方程组与本征值问题(09)：实对称矩阵的本征值求法

\(\newcommand{\K}{\mathbb{K}}\newcommand{\R}{\mathbb{R}}\newcommand{\C}{\mathbb{C}}\newcommand{\b}{\mathbf}\newcommand{\bi}{\boldsymbol}\newcommand{\rank}[1]{\text{rank}\left(#1 \right)}\newcommand{\dim}[1]{\text{dim}\left(#1 \right)}\newcommand{\diag}[1]{\text{diag}\left(#1 \right)}\newcommand{\det}[1]{\text{det}\left(#1 \right)}\newcommand{\Det}[1]{\left|\begin{matrix} #1 \end{matrix}\right|}\newcommand{\set}[1]{\left\{ #1 \right\}}\newcommand{\abs}[1]{\left| #1 \right|}\newcommand{\norm}[1]{\left|\!\left| #1 \right|\!\right|}\)本节针对实对称矩阵这一特殊形式, 讨论其本征值如何求解. 任何一本线代教材都会提到实对称矩阵的对角化方案, 本文的求解办法就从这里开始.

一般实对称矩阵的Jacobi 算法

由于任何实对称矩阵原则上都能对角化, 所以直接利用 Givens 变换逐渐地消去非对角元. 假设\(\b{A}\in\R^{n\times n}\)是实对称矩阵, 对于一对指标\(p,q\)(不妨设\(p<q\)), 用 Givens 变换设法消去交叉项: \[ \begin{align}&\begin{pmatrix} c & -s \\ s & c \end{pmatrix} \begin{pmatrix} a_{pp} & a_{pq} \\ a_{pq} & a_{qq} \end{pmatrix} \begin{pmatrix} c & s \\ -s & c \end{pmatrix} \\=&\begin{pmatrix} c^2a_{pp}+s^2a_{qq}-2csa_{pq} & cs(a_{pp}-a_{qq})+(c^2-s^2)a_{pq} \\ cs(a_{pp}-a_{qq})+(c^2-s^2)a_{pq} & c^2a_{pp}+s^2a_{qq}+2csa_{pq} \end{pmatrix} \end{align} \] 其中\(c,s\)是余弦和正弦值的简写. 令对角项为\(0\), 得到 \[ \frac{\tan\theta}{\tan^2\theta-1}=\frac{a_{pq}}{a_{pp}-a_{qq}} \] 本征值为 \[ \lambda_\pm=\frac{a_{pp}+a_{qq}\pm\sqrt{(a_{pp}-a_{qq})^2+4a_{pq}^2}}{2} \] Jacobi 算法计算量还是\(O(n^3)\), 但是它比较严格, 所以十分稳定.

Sturm 序列

所谓的 Sturm 序列适用于实对称的三对角矩阵. 假设三对角矩阵\(\b{T}\)的对角元为\(d_k\), 副对角线的元为\(b_k\), 它的特征多项式可以递推求解: \[ p_k(x)=(d_k-x)p_{k-1}(x)-b_{k-1}^2p_{k-2}(x) \\ p_0(x)\equiv1,p_{1}(x)=d_1-x \] 其中\(p_k(x)\)可以理解为\(k\)阶主子矩阵的特征多项式. 这些特征多项式构成的序列, 称为 Sturm 序列, 它有一些重要的性质:

1. Sturm 序列中, 某一项\(p_k\)的根和它前一项\(p_{k-1}\)的根排序, 两者的根一定交错出现
2. 任取\(\mu\in\R\), 序列\(\set{p_0(\mu),p_1(\mu),\cdots,p_n(\mu)}\), 它从前往后变号的次数^[1]\(s(\mu)\)就是\(p_n\)的根中严格小于\(\mu\)的数目.

根据 Gershgorin 圆盘定理, \(p_n\)的根一定在如下区间 \[ \sigma(\b{T})\subset\bigcup_{k=1}^n[d_k-\abs{b_k}-\abs{b_{k-1}}, d_k+\abs{b_k}+\abs{b_{k-1}}] \] \(b\)的下标只有\(1\sim n-1\), 可以约定\(b_0=b_n=0\).

于是问题简化为了一个查找问题, 只需从这些闭区间中找到最小左值\(\alpha_0\)和最大右值\(\beta_0\), 只需要从\([\alpha_0,\beta_0]\)中搜索\(p_n(x)=0\)的根就可以了. 比较稳定的方案当然是对分法.

注意到: 前面提到的\(s(\mu)\)是一个非严格递增的函数, 可以将它作为判据, \(\text{e.g.}\)计算第\(k\)小的本征值:

1. \([\alpha,\beta]=[\alpha_0,\beta_0]\), 根据搜索区间\([\alpha,\beta]\)确定中值\(\gamma\);
2. 计算\(s(\gamma)\)并与\(k\)比较:
   1. \(s(\gamma)\ge k\), 说明\([\alpha_0,\gamma)\)内已经有不少于\(k\)个根了, 此时赋值\(\beta=\gamma\);
   2. \(s(\gamma)< k\), 说明\([\alpha_0,\gamma)\)内不足\(k\)个根, 此时赋值\(\alpha=\gamma\);
3. 判断终止条件, 不终止则按照新的\(\alpha,\beta\)继续搜索.

稀疏矩阵的 Lanczos 方法

如果待计算的矩阵是一个庞大的稀疏矩阵, 那么前面所列举的方法并不能有效地进行操作. 这时, 如果仅需要求解矩阵的部分而不是全部本征对. 那么可以使用 Krylov 子空间迭代的方法.

迭代使用的 Krylov 子空间为\(K_m(\b{A};\bi{v})\), 其中\(\b{A}\)就是待计算的稀疏矩阵, \(\bi{v}\)是可由我们任选的非零矢量. 我们希望从子空间中逐步构建正交完备基矢, 令 \[ \bi{v}_0=\b{0}\qquad \bi{v}_1=\frac{\bi{v}}{\norm{\bi{v}}} \] 从这两个矢量开始构建起一套子空间的正交基(一个长方正交矩阵) \[ \b{V}=(\bi{v}_1,\bi{v}_2,\cdots,\bi{v}_m)\in\R^{n\times m} \] 原始矩阵在这套正交基下的新表示为: \(t_{jk}=\bi{v}_j^\b{T}\b{A}\bi{v}_k\).

新矩阵当然是对称的, 更进一步地, 我们希望它是三对角的: \[ \b{T}=\begin{pmatrix} \alpha_1 & \beta_1 & \\ \beta_1 & \alpha_2 & \ddots \\ & \ddots& \ddots & \ddots & \\ & & \ddots & \alpha_{m-1} & \beta_{m-1} \\ & & &\beta_{m-1} & \alpha_m \end{pmatrix} \] 由于\(\b{A}\b{V}=\b{V}\b{T}\), 得 \[ \b{A}\bi{v}_k=\beta_{k-1}\bi{v}_{k-1}+\alpha_k\bi{v}_k+\beta_k\bi{v}_{k+1} \] 它实际上规定了迭代方程, 两边左乘\(\bi{v}_k\), 得到\(\alpha_k\): \[ \begin{align} \alpha_k&=\bi{v}_k^\b{T}\b{A}\bi{v}_k \\ \beta_k\bi{v}_{k+1}&=\b{A}\bi{v}_k-\alpha_k\bi{v}_k-\beta_{k-1}\bi{v}_{k-1} \end{align} \] \(\beta_k\)又该如何确定? 假设已知\(\beta_{k-1}\), 对右边取模即可. 为了保证递推式的成立, 规定\(\beta_0=0\).

完整递推过程如下(\(k=1,\cdots,m-1\)):

1. \(\alpha_k=\bi{v}_k^\b{T}\b{A}\bi{v}_k\);
2. \(\bi{w}_k=\b{A}\bi{v}_k-\alpha_k\bi{v}_k-\beta_{k-1}\bi{v}_{k-1}\);
3. \(\beta_k=\norm{\bi{w}_k},\bi{v}_{k+1}=\displaystyle\frac{\bi{w}_k}{\beta_{k}}\);
4. \(k+\!\!=1\);
5. 如果\(k\)循环结束, 计算\(\alpha_m=\bi{v}_m^\b{T}\b{A}\bi{v}_m\).

1. 可以约定: 零值相对前一项永远算作变号, 它的后一项相对零值永不算变号 ↩︎