线性方程组与本征值问题(08)：关于线性代数的进一步讨论
\(\newcommand{\K}{\mathbb{K}}\newcommand{\R}{\mathbb{R}}\newcommand{\C}{\mathbb{C}}\newcommand{\b}{\mathbf}\newcommand{\bi}{\boldsymbol}\newcommand{\mdl}[1]{\left|\!\left|#1\right|\!\right|}\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{\set}[1]{\left\{ #1 \right\}}\newcommand{\abs}[1]{\left| #1 \right|}\)对本征值问题的讨论即将收尾, 在此之前, 有必要进一步将前置的数学工具(尽管不一定用得上)陈列如下.
对角化和 Jordan 标准型
正规矩阵
在第4节讨论数学时, 我们提到了重数与亏损还有对角化的联系: 亏损等价于不可对角化. 可对角化的矩阵非常常见, 例如自伴矩阵和酉矩阵, 它们有一个共同的特点: 都和自己的共轭对易 \[ \b{A}\b{A}^\dagger=\b{A}^\dagger\b{A} \] 这类矩阵被称为 正规矩阵. 可以证明: 能酉对角化的矩阵必定正规: \[ \begin{cases} \begin{aligned} \b{U}^\dagger\b{A}\b{U}&=\b{D} \\ &\Downarrow\dagger \\ \b{U}^\dagger\b{A}^\dagger\b{U}&=\b{D}^\dagger \end{aligned} \end{cases} \overset{multiply}{\Longrightarrow} \begin{cases} \begin{aligned} \b{U}^\dagger\b{A}\b{A}^\dagger\b{U}&=\b{D}\b{D}^\dagger \\ \\ \b{U}^\dagger\b{A}^\dagger\b{A}\b{U}&=\b{D}^\dagger\b{D} \end{aligned} \end{cases} \] 对角矩阵天然地和共轭对易, 不同顺序乘法的右侧相等, 因此左侧也相等, 得证.
同理, 正规矩阵必定能酉对角化, 先考虑 Schur 分解定理: \[ \begin{cases} \begin{aligned} \b{U}^\dagger\b{A}\b{U}&=\b{H} \\ &\Downarrow\dagger \\ \b{U}^\dagger\b{A}^\dagger\b{U}&=\b{H}^\dagger \end{aligned} \end{cases} \overset{multiply}{\Longrightarrow} \begin{cases} \begin{aligned} \b{U}^\dagger\b{A}\b{A}^\dagger\b{U}&=\b{H}\b{H}^\dagger \\ \\ \b{U}^\dagger\b{A}^\dagger\b{A}\b{U}&=\b{H}^\dagger\b{H} \end{aligned} \end{cases}\Longrightarrow \b{H}\b{H}^\dagger=\b{H}^\dagger\b{H} \] 其中\(\b{H}\in\C^{n\times n}\)是一个上三角矩阵: \[ h_{jk}\begin{cases}\displaystyle \neq0\qquad j\le k \\\displaystyle =0\qquad j>k \end{cases} \qquad (h^\dagger)_{jk}\begin{cases}\displaystyle \neq0\qquad j\ge k \\\displaystyle =0\qquad j<k \end{cases} \] 即 \[ [\b{H}\b{H}^\dagger]_{jj}=\sum_{l=j}^n\abs{h_{jl}}^2 = [\b{H}^\dagger\b{H}]_{jj}=\sum_{l=1}^j\abs{h_{lj}}^2 \]
- 取\(j=1\), 得到\(h_{12},h_{13},\cdots,h_{1n}\)的模方和为\(0\)
- 取\(j=2\), 得到\(h_{23}, h_{24},\cdots,h_{2n}\)的模方和为\(\abs{h_{12}}^2\), 此前求出为\(0\)
- 取\(j=3\), 得到\(h_{34},\cdots,h_{3n}\)的模方和为\(\abs{h_{13}}^2+\abs{h_{23}}^2\), 此前已求出为\(0\)
- ...
像这样逐步推导, 得到所有非对角元都为\(0\), 即上三角矩阵\(\b{H}\)是对角矩阵, 得证.
Jordan 标准型
正规矩阵等价于可以酉对角化, 而非正规的矩阵就无法酉对角化乃至无法(相似)对角化了, 例如下面的矩阵 \[ \b{J}_\nu(\lambda)=\begin{pmatrix} \lambda & 1 & 0 & \cdots & 0 \\ 0 & \lambda & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ \vdots & & \ddots & \lambda & 1 \\ 0 & \cdots & \cdots & 0 & \lambda \end{pmatrix} \] 其中\(\nu\)为阶数. 不难验证, 它具有\(\nu\)重本征值\(\lambda\), 但本征子空间仅为\(1\)维, 显然是亏损矩阵, 不可对角化. 这样的矩阵叫做 Jordan 块矩阵(\(\nu\ge2\)); 对于\(\nu=1\), 约定 Jordan 块矩阵退化为\(\lambda\).
我们不加证明地给出如下定理:
对于任意\(\b{A}\in\C^{n\times n}\), 总有满秩矩阵\(\b{X}\in\C^{n\times n}\), \(\text{s.t.}\) \[ \b{X}^{-1}\b{A}\b{X}=\diag{\b{J}_{\nu_1}(\lambda_1),\cdots,\b{J}_{\nu_{k}}(\lambda_k)} \] \(\lambda_1\sim\lambda_k\)是矩阵的本征值(可能相等).
这个对角形式称作矩阵\(\b{A}\)的 Jordan 标准型. 需要说明的是
一个本征值可能对应多个 Jordan 块, \(\text{e.g.}\) \[ \begin{pmatrix} \lambda & 1 & 0 \\ 0 & \lambda & 0 \\ 0 & 0 & \lambda \end{pmatrix} \] 这其实就是两个 Jordan 块, 一个二阶和一个一阶;
一个本征值对应的 Jordan 块数目等于它的独立本征矢数目(几何重数), 整个矩阵的 Jordan 块数就是总的独立本征矢数目;
矩阵可对角化\(\Leftrightarrow\)标准型中的 Jordan 块都是一阶.
从这些结论能够立即得到: 本征值全不等(即有\(n\)个不等的本征值)的矩阵一定是可对角化的 , 这个也可由重数的结论导出, 因为几何重数至少为\(1\), 因此代数重数全为\(1\)时, 几何重数不会比代数重数更小了.
简而言之, 一个矩阵是否可对角化, 与它的 Jordan 标准型挂钩.
Jordan 标准型的不稳定性
请看如下矩阵 \[ \b{A}=\begin{pmatrix} 1 & 1 \\ \epsilon & 1 \end{pmatrix} \] 当\(\epsilon=0\)时, 它已经是 Jordan 标准型了; 但\(\epsilon\neq0\)时, 它的特征值为 \[ \lambda_\pm=1\pm\sqrt{\epsilon} \] 对应 Jordan 标准型为 \[ \b{J}=\begin{pmatrix} 1+\sqrt{\epsilon} & 0 \\ 0 & 1-\sqrt{\epsilon} \end{pmatrix} \] 换言之, 矩阵元\(j_{12}(\epsilon)\)关于\(\epsilon\)是非连续的, 这表明它在数值计算中非常不稳定, 因此实际运算中, 很少利用 Jordan 标准型完成某些任务, 但前面的大量分析说明它具有一定的指导意义.
本征值的分布
Gershgorin 对于本征值在复平面上的分布提出了一系列定理, 其中最为重要的是:
对任意\(\b{A}\in\C^{n\times n}\), 构建\(n\)个圆盘 \[ \mathcal{R}_j=\set{z\in\C\Bigg|\abs{z-a_{jj}}\le\sum_{k=1,k\neq j}^{n}\abs{a_{jk}}} \] 这种圆盘叫 Gershgorin 圆盘, 则矩阵的本征谱满足 \[ \sigma(\b{A})\subset\bigcup_{j=1}^n\mathcal{R}_j \]
这个称为 Gershgorin 圆盘定理, 它说明矩阵的本征值一定落在这些圆盘中. 特别是自伴矩阵, 它的本征值一定是实数, 圆盘退化为实轴上的区间.
Krylov 子空间
解线性方程组和 Krylov 子空间
这个方法可能最早来自求解线性方程组. 有时候线性方程组的系数矩阵是大型稀疏矩阵, 而求解又不需要严格成立, 这时可以设法在一个子空间中求它的近似解.
线代中有一条 Cayley-Hamilton 定理: \(p_\b{A}(\b{A})=\b{O}\), 其中\(p_\b{A}(x)\)是矩阵\(\b{A}\)的特征多项式. 具体地, 可以写作 \[ c_{n}\b{A}^n+c_{n-1}\b{A}^{n-1}+\cdots+1=\b{O} \] 进一步假设\(\b{A}\)是可逆的, 那么\(\b{A}^{-1}=-c_{n}\b{A}^{n-1}-c_{n-1}\b{A}^{n-2}-\cdots\), 而线性方程组也可表示为 \[ \b{A}\bi{x}=\bi{b}\Rightarrow\bi{x}=-\sum_{j=0}^{n-1}c_{j+1}\b{A}^j\bi{b} \] 这实际上给了我们一定的灵感. 我们深知方程求根和优化问题的相似性: \[ \text{solve }\b{A}\bi{x}=\bi{b}\Leftrightarrow\underset{\bi{x}\in\K^n}{\displaystyle\text{min}}\mdl{\bi{b}-\b{A}\bi{x}} \] 那么近似求根, 也就相当于在更低维(独立矢量和包含信息更少)的子空间进行优化: \[ \text{approx-solve }\b{A}\bi{x}=\bi{b}\Leftrightarrow\underset{\bi{x}\in K_m}{\displaystyle\text{min}}\mdl{\bi{b}-\b{A}\bi{x}} \] 其中, \(K_m(\b{A};\bi{b})\)是一个子空间, 它由如下矢量生成: \[ K_m(\b{A};\bi{b}):=\text{span}\{\bi{b},\b{A}\bi{b},\cdots,\b{A}^{m-1}\bi{b} \} \] 称为 Krylov 子空间. 根据 Cayley-Hamilton 定理, 当\(m=n+1\)时, 这个矢量组就必定能互相线性表出, 因此\(m\)无论取多少, 总有\(\dim{K_m}\le\min{\{m,n\}}\).
Krylov 子空间的实例: 共轭梯度法
之前的数学准备不足, 现在可以介绍一下共轭梯度法和 Krylov 子空间的联系了. 我们考虑一个二次型函数 \[ f(\bi{x})=-\bi{b}^\b{T}\bi{x}+\frac{1}{2}\bi{x}^\b{T}\b{A}\bi{x} \] 它的优化条件是\(\bi{r}=\bi{b}-\b{A}\bi{x}=0\), 或者\(\underset{\bi{x}\in\R^n}{\min}\mdl{\bi{b}-\b{A}\bi{x}}\).
共轭梯度法采用的 Krylov 子空间为\(K_m(\b{A};\bi{r}_0)\), 每一步搜索方向都从这个子空间中挑选. 单纯的\(\b{A}^j\bi{r}_0\)往往并不是正交的, 而且计算不便, 所以搜索方向[1]也不是单纯的\(\b{A}^j\bi{r}_0\): \[ \bi{p}_k=g_k(\bi{r}_0,\b{A}\bi{r}_0, \cdots,\b{A}^{k-1}\bi{r}_0) \] 每次搜索\(\bi{x}_k+\!\!=\alpha_k\bi{p}_k, \bi{r}_k-\!\!=\alpha_k\b{A}\bi{p}_k\), 都保证残差\(\bi{r}\)完全失去一个方向上的值: \(\bi{r}_{k+1}^\b{T}\bi{p}_k=0\).
这时, 残差(等价于位置)的递推式已经给出: \[ \bi{r}_{k+1}=\bi{r}_k-\alpha_k\b{A}\bi{p}_k \qquad \alpha_k\text{满足:}\quad\bi{r}_{k+1}^\b{T}\bi{p}_k=0 \] 现在有必要考虑搜索方向如何更新, 一种思路是, 保证搜索方向在 Krylov 子空间构成正交矢量组 \[ \bi{p}_{k+1}^\b{T}\bi{p}_k=0 \] 这看起来非常好, 承担着 Schmidt 正交化的任务; 但是稍加思索就能发现这是不合理的. 每一步残差扣除的部分为\(\b{A}\bi{p}_k\), 需要避免新的搜索方向"破坏"这个既有成果, 因此 \[ \bi{p}_{k+1}^\b{T}\b{A}\bi{p}_k=0 \] 另外, \(\bi{p}_{k+1}\)当然和新的残差有关[2], 你总是从残差的大小和方向确定优化方向, 所以\(\bi{p}_{k+1}\)由\(\bi{r}_{k+1}\)和\(\bi{p}_k\)决定 \[ \bi{p}_{k+1}=\bi{r}_{k+1}+\beta_k\bi{p}_k\qquad\beta_k满足:\quad\bi{p}_{k+1}^\b{T}\b{A}\bi{p}_k=0 \]
这些正是共轭梯度法的基本思路. 共轭梯度法的各级搜索方向在 Krylov 子空间中构成了共轭正交的矢量组. 实际上还可以证明, 残差也在子空间中, 并且是正交矢量组.