线性方程组与本征值问题(05)：原始的 QR 算法

对于矩阵的谱, 有几种求法:

1. 解特征多项式
2. 任何一本线代教材都会讲到的对角化手续
3. ...

第一种方法求方程的数值解, 原则上需要迭代; 第二种不一定适用, 因为矩阵未必能够对角化. 是否可以找到一个求一般矩阵的本征值的"直接解法"呢?

很不幸的是, 这一点原则上是不可能的. 假设存在一个有限步骤的所谓"直接算法", 即通过有限多次的初等代数运算就可以获得一般矩阵的本征值, 就意味着存在一个任意阶多项式求根的公式, 而这是与 Abel 对高次方程的研究相矛盾的.

因此, 本征值问题的一般解仍需要迭代方法, 本节将要介绍的是著名的 QR 算法.

QR 分解

在引入 QR 算法前, 先介绍实矩阵的 QR 分解:

如果矩阵\(\mathbf{A}\in\mathbb{R}^{n\times n}\)并且满足 \[ \mathbf{A}=\mathbf{Q}\mathbf{R} \] 其中, \(\mathbf{Q}\)和\(\mathbf{R}\)分别是正交矩阵和上三角矩阵, 则称之为矩阵\(\mathbf{A}\)的一个 QR 分解.

非奇异矩阵

QR 分解的存在性可以利用所谓的 Schmidt 正交归一化来说明. 考虑\(\mathbf{A}\)的列向量组\(\{\pmb{\alpha}_1,\pmb{\alpha}_2,\cdots,\pmb{\alpha}_n\}\)​线性无关, 那么有所谓的 Schmidt 正交归一化手续 得到一个正交归一基: \[ \begin{align} \pmb{\gamma}_1 &= \frac{\pmb{\beta}_1}{|\!|\pmb{\beta}_1|\!|} \qquad \pmb{\beta}_1 = \pmb{\alpha}_1 \\ \pmb{\gamma}_2 &= \frac{\pmb{\beta}_2}{|\!|\pmb{\beta}_2|\!|} \qquad \pmb{\beta}_2 = \pmb{\alpha}_2-(\pmb{\alpha}_2, \pmb{\gamma}_1)\pmb{\gamma}_1 \\ &\qquad\qquad\vdots \\ \pmb{\gamma}_k &= \frac{\pmb{\beta}_k}{|\!|\pmb{\beta}_k|\!|} \qquad \pmb{\beta}_k = \pmb{\alpha}_k-\sum_{j=1}^{k-1} (\pmb{\alpha}_k, \pmb{\gamma}_j)\pmb{\gamma}_j \\ &\qquad\qquad\vdots \end{align} \] 其中, \(k=1,2,\cdots,n\), \((\sim_1,\sim_2)\)表示向量的内积. 如果定义新的矩阵元 \[ r_{jk}=\begin{align*}\begin{cases}\displaystyle (\pmb{\alpha}_k,\pmb{\gamma}_j)\qquad &j<k \\\\\displaystyle |\!|\pmb{\beta}_k|\!| &j=k \\\\\displaystyle 0 &j>k \end{cases}\end{align*} \] 则有 \[ \pmb{\alpha}_k=\sum_{j=1}^n\pmb{\gamma}_{j}r_{jk} \qquad\Rightarrow\qquad \begin{pmatrix}\pmb{\alpha}_1&\cdots&\pmb{\alpha}_n\end{pmatrix}= \begin{pmatrix}\pmb{\gamma}_1&\cdots&\pmb{\gamma}_n\end{pmatrix} \begin{pmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\ & r_{22} & \cdots & r_{2n} \\ && \ddots & \vdots \\ &&& r_{nn} \end{pmatrix} \] 即\(\mathbf{A}=\mathbf{Q}\mathbf{R}\). 可见任何非奇异方阵都存在 QR 分解.

奇异矩阵

考虑\(\mathbf{A}\in\mathbb{R}^{n\times n}\), 且\(\text{rank}(\mathbf{A})=m\). 和上面相当类似, 对于前\(k\)列构成的向量组\(\{\pmb{\alpha}_1,\cdots,\pmb{\alpha}_k\}\), 设它的秩为\(m_k\), 显然\(m_1=1\), \(m_n=m\). 可以找到\(\{\pmb{\varepsilon}_1,\cdots,\pmb{\varepsilon}_{m_k}\}\)这\(m_k\)个正交归一向量构成的基来线性表示: \[ \pmb{\alpha}_k=\sum_{j=1}^{m_k}(\pmb{\alpha}_k,\pmb{\varepsilon}_j)\pmb{\varepsilon}_j \] 因为\(m_k\le k\), 所以\((\pmb{\alpha}_k,\pmb{\varepsilon}_j)\)仍是一个上三角矩阵的矩阵元. 与非奇异的区别在于, 这里的上三角矩阵即使在上三角部分也会出现零元, 甚至必定出现零元. 它的后\(n-m\)行都是零行. 另外还需注意, 到目前为止, \(\pmb{\varepsilon}_j\)无法构成正交矩阵. 它只有\(m\)个矢量. 我们再找\(n-m\)个列矢量, 使得所有\(n\)个构成正交归一基, 这样才能形成正交矩阵.

综上所述, 此时仍存在分解, 只是\(\mathbf{Q}\)的前\(m\)列就足以表示\(\mathbf{A}\), 因而\(\mathbf{R}\)也只有前\(m\)行才可能出现非零元, 后\(n-m\)行无论对角线上下都是零元.

标准 QR 迭代

前面已经证明了方阵普遍具备 QR 分解, 现在可以介绍 QR 迭代的基本步骤了.

要求

初始矩阵\(\mathbf{A}\), 初始正交矩阵\(\mathbf{Q}^{(0)}\). 若\(\mathbf{Q}^{(0)}=\mathbf{I}\), 则\(\mathbf{T}^{(0)}=\mathbf{A}\).

步骤

1. for \(k=1,2,\cdots\) do

   1. 对矩阵\(\mathbf{T}^{(k-1)}\)进行 QR 分解, 即

   \[ \mathbf{T}^{(k-1)}=\mathbf{Q}^{(k)}\mathbf{R}^{(k)} \]

   2. 计算\(\mathbf{T}^{(k)}=\mathbf{R}^{(k)}\mathbf{Q}^{(k)}\).
   3. 这一步进行终止条件判定.

2. 输出结果

复杂度分析

\(n\)阶方阵乘法每确定一个矩阵元为\(O(n)\), 共\(n^2\)个, 因此每一轮迭代的复杂度为\(O(n^3)\).

实 Schur 分解

作为迭代算法, QR 迭代当然是近似的. 它的每一步结果是 \[ \mathbf{T}^{(k)}=\left(\mathbf{Q}^{(0)}\cdots\mathbf{Q}^{(k)}\right)^\mathbf{T}\mathbf{A}\left(\mathbf{Q}^{(0)}\cdots\mathbf{Q}^{(k)}\right) \] 实际上, 如果\(\mathbf{A}\)存在复本征值, 那么\(\mathbf{T}^{(k)}\)永不可能是真的上三角的. 一方面, 实矩阵之积必定是实矩阵, 另一方面, \(\mathbf{A}\)相似变换而成的上三角矩阵, 对角元都是本征值. 因此, 具有复本征值的矩阵不可能在实数域上严格上三角化.

尽管如此, 我们可以设法让\(\mathbf{A}\)"尽量上三角化":

\(\forall\mathbf{A}\in\mathbb{R}^{n\times n}\), 存在\(\mathbf{Q}\)为正交矩阵且 \[ \mathbf{Q}^\mathbf{T}\mathbf{A}\mathbf{Q}= \begin{pmatrix} R_{11} & R_{12} & \cdots & R_{1m} \\ 0 & R_{22} & \cdots & R_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & R_{mm} \end{pmatrix} \] 其中\(R_{ii}\in\mathbb{R}\)表示所有实本征值, 或\(R_{ii}\in\mathbb{R}^{2\times2}\). 这样的形式称作矩阵\(\mathbf{A}\)的 实 Schur 形式.

不难看出, 实 Schur 形式很可能不是上三角的, 它的对角元处允许\(2\times2\)的块矩阵存在(其实也只允许对角线下方的副对角线存在), 这代表一对共轭的复本征值^[1].

另外, 我们希望非对角块包含且仅包含无法对角化的复本征值, 后面我们会看到, 由于算法的缺陷, 实本征值也可能没有对角化而形成对角块, 这种情况不算严格意义上的实 Schur 形式.

下面的定理给出实 Schur 形式和 QR 分解的联系:

正交矩阵\(\mathbf{Q}\)可以通过 QR 迭代获得, 即 \[ \mathbf{Q}=\lim_{k\to\infty}\left[\mathbf{Q}^{(0)}\mathbf{Q}^{(1)}\cdots\mathbf{Q}^{(k)} \right] \] 这样得到一个矩阵的迭代序列\(\{\mathbf{T}^{(i)}\}\), 本征值满足\(|\lambda_1|\ge|\lambda_2|\ge\cdots\ge|\lambda_n|\). 则它的非对角元具有如下收敛速度: \[ \left|t_{ij}^{(k)}\right|\sim\left|\frac{\lambda_i}{\lambda_j}\right|^k, i>j \]

1. 我们好像没有说明共轭的根据, 这是因为本征方程是实的, 取共轭不变, 因此根也是共轭不变的. ↩︎