QR算法的收敛原理
QR算法求特征值的过程就是，不断地对目标矩阵进行 QR 分解，同时交换 Q 矩阵和 R 矩阵相乘得到新的目标矩阵。 持续迭代之后，目标矩阵就会收敛到一个上三角矩阵上。然后我们就可以直接从该上三角矩阵的对角线上读出特征值，并且还是按照绝对值大小排序的，最大的在左上角，最小的在右下角。 这看起来真的十分优雅。初次见到这个算法的时候，我就有几个疑问：
- 其一，这个迭代的目标矩阵与原矩阵的特征值一样吗？关于这个问题， 我们在介绍特征值与特征向量算法的时候讨论过，因为它们是相似的，所以特征值保持不变。
- 其二，上三角矩阵的对角线上元素就是它的特征值吗？这点通过特征多项式和行列式很容易证明，本文将作简单介绍。
- 其三，这个算法的适用范围有多大？假如说矩阵是个方阵并且满秩，那么它就一定会收敛到上三角矩阵吗？我查阅了很多资料之后， 发现比较确定的是，如果矩阵可以对角化，那么就会收敛到上三角矩阵。
1. 上三角矩阵的特征值
对于上三角矩阵 \(\boldsymbol{U} = \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} \end{bmatrix}\) 其特征多项式是矩阵 \(\boldsymbol{U} - \lambda\boldsymbol{I}\) 的行列式 \(\begin{vmatrix} u_{11} - \lambda & u_{12} & \cdots & u_{1n} \\ & u_{22} - \lambda & \cdots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} - \lambda \end{vmatrix}\)。我们知道行列式可以按照某一行或者某一列展开成余子式的形式，由于 \(\boldsymbol{U} - \lambda\boldsymbol{I}\) 也是一个上三角矩阵， 它的第一列除了第一行其余都是 0，所以按照第一列展开之后就只有一项: $$ \det(\boldsymbol{U} - \lambda\boldsymbol{I}) = (u_{11} - \lambda) \begin{vmatrix} u_{22} - \lambda & \cdots & u_{2n} \\ & \ddots & \vdots \\ & & u_{nn} - \lambda \end{vmatrix} = (u_{11} - \lambda) \det(\boldsymbol{U} - \lambda\boldsymbol{I})_{11} $$ 接着对余子式 \(\det(\boldsymbol{U} - \lambda\boldsymbol{I})_{11}\) 展开，最终可以得到: $$ \det(\boldsymbol{U} - \lambda\boldsymbol{I}) = \prod_1^n (u_{ii} - \lambda) $$ 显然，对角线上的元素 \(u_{ii}\) 就是特征多项式 \(\prod_1^n (u_{ii} - \lambda) = 0\) 的根。所以上三角矩阵对角线上的元素就是其特征值。
2. 可逆矩阵 QR 分解的唯一性
起初我以为只要一个矩阵是可逆的那么它的 QR 分解就是唯一的，后来发现还少了一个条件，需要上三角矩阵\(\boldsymbol{R}\) 的对角线元素都是正数。 大多数 QR 分解实现也会添加这么个约束，来保证结果的唯一性。基于这样的实现，QR 算法的收敛过程才能更稳定。 所以，我们给 QR_Householder 也加上了这个约束。 下面我们来证明这个分解的唯一性。
假设矩阵 \(\boldsymbol{A}\) 可逆并且有两个 QR 分解，\(\boldsymbol{A} = \boldsymbol{Q}_1\boldsymbol{R}_1 = \boldsymbol{Q}_2 \boldsymbol{R}_2\)， 其中 \(\boldsymbol{Q}_1, \boldsymbol{Q}_2\) 正交，上三角矩阵\(\boldsymbol{R}_1,\boldsymbol{R}_2\)的对角线元素均为正数。由于矩阵 \(\boldsymbol{A}\) 可逆， 所以 \(\boldsymbol{R}_1, \boldsymbol{R}_2\) 也可逆，有: $$ \boldsymbol{Q}_1\boldsymbol{R}_1 = \boldsymbol{Q}_2 \boldsymbol{R}_2 \Longrightarrow \boldsymbol{Q}_1 = \boldsymbol{Q}_2\boldsymbol{R}_2\boldsymbol{R}_1^{-1} $$ 由于 \(\boldsymbol{Q}_1, \boldsymbol{Q}_2\) 正交 $$ \boldsymbol{I} = \boldsymbol{Q}_1^T\boldsymbol{Q}_1 = (\boldsymbol{R}_2\boldsymbol{R}_1^{-1})^T\boldsymbol{Q}_2^T \boldsymbol{Q}_2^T (\boldsymbol{R}_2\boldsymbol{R}_1^{-1}) = (\boldsymbol{R}_2\boldsymbol{R}_1^{-1})^T (\boldsymbol{R}_2\boldsymbol{R}_1^{-1}) $$ 所以矩阵 \(\boldsymbol{R}_2\boldsymbol{R}_1^{-1}\) 正交。由于上三角矩阵的逆仍然是上三角矩阵，两个上三角矩阵相乘也还是上三角的，所以 \(\boldsymbol{R}_2\boldsymbol{R}_1^{-1}\) 也是上三角的。 所以矩阵 \(\boldsymbol{R}_2\boldsymbol{R}_1^{-1}\) 只能是对角矩阵，并且对角线元素为 \(±1\)。
我们刚刚证明过，上三角矩阵对角线的元素就是它的特征值，如果 \(r_1(i, i)\) 是 \(\boldsymbol{R}_1\) 对角线的元素，那么 \(1 / r_1(i, i)\) 就是 \(\boldsymbol{R}_1^{-1}\) 对角线的元素。 矩阵 \(\boldsymbol{R}_2\boldsymbol{R}_1^{-1}\) 对角线上的元素，实际就是矩阵 \(\boldsymbol{R}_2\) 和 \(\boldsymbol{R}_1^{-1}\) 对角线上元素的乘积。 而矩阵\(\boldsymbol{R}_1, \boldsymbol{R}_2\)对角线的元素又都是正数，所以 \(\boldsymbol{R}_2\boldsymbol{R}_1^{-1}\) 对角线元素也一定是正数。
那么\(\boldsymbol{R}_2\boldsymbol{R}_1^{-1}\)只能是单位矩阵 \(\boldsymbol{I}\)，所以 \(\boldsymbol{R}_2 = \boldsymbol{R}_1\)，因此矩阵 \(\boldsymbol{Q}_1, \boldsymbol{Q}_2\) 也相等。 所以约束上三角矩阵 \(\boldsymbol{R}\) 对角线元素都是正数的情况下，可逆矩阵的 QR 分解是唯一的。
3. QR 迭代与矩阵的幂
假设我们对原始矩阵 \(\boldsymbol{A}_1\) 第一次进行 QR 分解，得到 \(\boldsymbol{A}_1 = \boldsymbol{Q}_1\boldsymbol{R}_1\)， 交换相乘得到 \(\boldsymbol{A}_2 = \boldsymbol{R}_1\boldsymbol{Q}_1\)。这个过程相当于对矩阵 \(\boldsymbol{A}_1\) 进行了一次关于 \(\boldsymbol{Q}_1\) 的相似变换， 即，\(\boldsymbol{A}_2 = \boldsymbol{Q}_1^T\boldsymbol{A}_1\boldsymbol{Q}_1\)。不难写出 QR 算法的迭代的递推关系: $$ \begin{array}{rl} \boldsymbol{A}_{k+1} & = \boldsymbol{Q}_k^T \boldsymbol{A}_k \boldsymbol{Q}_k \\ & = (\boldsymbol{Q}_k^T\boldsymbol{Q}_{k-1}^T\cdots\boldsymbol{Q}_1^T) \boldsymbol{A}_1 (\boldsymbol{Q}_1\cdots\boldsymbol{Q}_{k-1}\boldsymbol{Q}_k) \\ & = \tilde{\boldsymbol{Q}_k}^T \boldsymbol{A}_1 \tilde{\boldsymbol{Q}_k} \\ \tilde{\boldsymbol{Q}_k} & = \boldsymbol{Q}_1\cdots\boldsymbol{Q}_{k-1}\boldsymbol{Q}_k \end{array} $$ 令 \(\tilde{\boldsymbol{R}_k} = \boldsymbol{R}_k\boldsymbol{R}_{k-1}\cdots\boldsymbol{R}_1\)，则有: $$ \begin{array}{rl} \tilde{\boldsymbol{Q}_k}\tilde{\boldsymbol{R}_k} & = \tilde{\boldsymbol{Q}_{k-1}}(\boldsymbol{Q}_k\boldsymbol{R}_k)\tilde{\boldsymbol{R}_{k-1}} \\ & = \tilde{\boldsymbol{Q}_{k-1}}\boldsymbol{A}_k\tilde{\boldsymbol{R}_{k-1}} \\ & = \tilde{\boldsymbol{Q}_{k-1}}(\tilde{\boldsymbol{Q}_{k-1}}^T\boldsymbol{A}_1\tilde{\boldsymbol{Q}_{k-1}})\tilde{\boldsymbol{R}_{k-1}} \\ & = \boldsymbol{A}_1\tilde{\boldsymbol{Q}_{k-1}}\tilde{\boldsymbol{R}_{k-1}} \end{array} $$ 如此迭代下去，我们就会得到一个有意思的结论 \(\tilde{\boldsymbol{Q}_k}\tilde{\boldsymbol{R}_k} = \boldsymbol{A}_1^k\)。
4. 可对角化矩阵的 QR 迭代过程的收敛性
假设矩阵 \(\boldsymbol{A}\) 可以对角化，存在一个可逆矩阵 \(\boldsymbol{V}\) 使得 \(\boldsymbol{A} = \boldsymbol{V}^{-1}\boldsymbol{\Lambda}\boldsymbol{V}\) 成立。 其中 \(\boldsymbol{\Lambda} = \text{diag}(\lambda_1, \cdots, \lambda_n)\)，并且 \(|\lambda_1| > |\lambda_2| > \cdots > |\lambda_n| > 0\)。 如果 \(\boldsymbol{V}\) 存在 LU 分解，即 \(\boldsymbol{V} = \boldsymbol{L}_v \boldsymbol{U}_v\)，\(\boldsymbol{L}_v\) 是个单位下三角阵对角线的元素都是1，\(\boldsymbol{V}^{-1}\)的QR分解为 \(\boldsymbol{V}^{-1} = \boldsymbol{Q}_v\boldsymbol{R}_v\)。那么有: $$ \begin{equation}\label{func_1} \begin{array}{rl} \boldsymbol{A}^k & = \boldsymbol{V}^{-1} \boldsymbol{\Lambda}^k \boldsymbol{V} \\ & = \boldsymbol{Q}_v\boldsymbol{R}_v\boldsymbol{\Lambda}^k\boldsymbol{L}_v \boldsymbol{U}_v \\ & = \boldsymbol{Q}_v\boldsymbol{R}_v(\boldsymbol{\Lambda}^k\boldsymbol{L}_v \boldsymbol{\Lambda}^{-k})\boldsymbol{\Lambda}^{k} \boldsymbol{U}_v \\ \end{array} \end{equation} $$ 由于 \(\boldsymbol{L}_v\) 是个单位下三角阵，\(\boldsymbol{\Lambda}^k\) 是个对角阵，所以\(\boldsymbol{\Lambda}^k\boldsymbol{L}_v \boldsymbol{\Lambda}^{-k}\) 是一个下三角矩阵， 其 \((i,j)\) 上的元素为 $$ \left(\boldsymbol{\Lambda}^k\boldsymbol{L}_v \boldsymbol{\Lambda}^{-k}\right)(i, j) = \begin{cases} 0, & i < j \\ 1, & i = j \\ \boldsymbol{L}_v(i,j) \left(\frac{\lambda_i}{\lambda_j}\right)^k, & i > j \end{cases} $$ 对于 \(i > j\)，有 \(|\lambda_i| < |\lambda_j|\)，所以 \(k \to \infty\) 时, \(\left(\frac{\lambda_i}{\lambda_j}\right)^k \to 0\)。 因此，我们可以将\(\boldsymbol{\Lambda}^k\boldsymbol{L}_v \boldsymbol{\Lambda}^{-k}\) 写成如下的形式 $$ \boldsymbol{\Lambda}^k\boldsymbol{L}_v \boldsymbol{\Lambda}^{-k} = \boldsymbol{I} + \boldsymbol{E}_k $$ $$ \lim_{k \to \infty} \boldsymbol{E}_k = \boldsymbol{0} $$ 带入式(\(\ref{func_1}\))，有 \(\boldsymbol{A}^k = \boldsymbol{Q}_v\boldsymbol{R}_v(\boldsymbol{I} + \boldsymbol{E}_k)\boldsymbol{\Lambda}^{k} \boldsymbol{U}_v = \boldsymbol{Q}_v(\boldsymbol{I} + \boldsymbol{R}_v\boldsymbol{E}_k\boldsymbol{R}_v^{-1})\boldsymbol{R}_v \boldsymbol{\Lambda}^{k} \boldsymbol{U}_v\)。 不妨设\(\boldsymbol{I} + \boldsymbol{R}_v\boldsymbol{E}_k\boldsymbol{R}_v^{-1}\) 的 QR 分解为 \(\boldsymbol{Q}_{E_k}\boldsymbol{R}_{E_k}\)， 因为 \(\boldsymbol{E}_k \to \boldsymbol{0}\)，所以 \(\boldsymbol{Q}_{E_k} \to \boldsymbol{I}, \boldsymbol{R}_{E_k} \to \boldsymbol{I}\)。 有
$$ \begin{equation}\label{func_2} \begin{array}{rl} \boldsymbol{A}^k & = \boldsymbol{Q}_v(\boldsymbol{Q}_{E_k}\boldsymbol{R}_{E_k})\boldsymbol{R}_v \boldsymbol{\Lambda}^{k} \boldsymbol{U}_v \\ & = \left(\boldsymbol{Q}_v\boldsymbol{Q}_{E_k} \boldsymbol{D}^{-1}\right) \left(\boldsymbol{D}\boldsymbol{R}_{E_k}\boldsymbol{R}_v \boldsymbol{\Lambda}^{k} \boldsymbol{U}_v\right) \end{array} \end{equation} $$上式中矩阵 \(\boldsymbol{D}\) 是对角线元素为 \(±1\) 的对角矩阵，使得\(\boldsymbol{D}\boldsymbol{R}_{E_k}\boldsymbol{R}_v \boldsymbol{\Lambda}^{k} \boldsymbol{U}_v\) 的对角线元素均为正数。显然，\(\boldsymbol{D}\boldsymbol{R}_{E_k}\boldsymbol{R}_v \boldsymbol{\Lambda}^{k} \boldsymbol{U}_v\) 是一个上三角矩阵，所以上式构成了\(\boldsymbol{A}^k\)的一个QR分解。 根据刚刚证明的 QR 分解的唯一性，我们有 $$ \begin{align} \tilde{\boldsymbol{Q}_k} & = \boldsymbol{Q}_v\boldsymbol{Q}_{E_k} \boldsymbol{D}^{-1} \\ \tilde{\boldsymbol{R}_k} & = \boldsymbol{D}\boldsymbol{R}_{E_k}\boldsymbol{R}_v \boldsymbol{\Lambda}^{k} \boldsymbol{U}_v \end{align} $$ 根据 QR 迭代关系式，有: $$ \begin{array}{rl} \boldsymbol{A}_{k+1} & = \tilde{\boldsymbol{Q}_k}^T\boldsymbol{A}\tilde{\boldsymbol{Q}_k} \\ & = \left(\boldsymbol{D}^{-T}\boldsymbol{Q}_{E_k}^T\boldsymbol{Q}_v^T\right) \left(\boldsymbol{V}^{-1} \boldsymbol{\Lambda}^k \boldsymbol{V}\right) \left(\boldsymbol{Q}_v\boldsymbol{Q}_{E_k} \boldsymbol{D}^{-1}\right) \\ & = \left(\boldsymbol{D}^{-T}\boldsymbol{Q}_{E_k}^T\boldsymbol{Q}_v^T\right) \left(\boldsymbol{Q}_v\boldsymbol{R}_v \boldsymbol{\Lambda}^k \boldsymbol{R}_v^{-1}\boldsymbol{Q}_v^{-1}\right) \left(\boldsymbol{Q}_v\boldsymbol{Q}_{E_k} \boldsymbol{D}^{-1}\right) \\ & = \boldsymbol{D}^{-T}\boldsymbol{Q}_{E_k}^T \boldsymbol{R}_v \boldsymbol{\Lambda}^k \boldsymbol{R}_v^{-1}\boldsymbol{Q}_{E_k} \boldsymbol{D}^{-1} \end{array} $$ 由于 \(k \to \infty\) 时，\(\boldsymbol{Q}_{E_k} \to \boldsymbol{I}\)，所以 $$ \lim_{k \to \infty} \boldsymbol{A}_{k+1} = \boldsymbol{D}^{-T}\boldsymbol{R}_v \boldsymbol{\Lambda}^k \boldsymbol{R}_v^{-1}\boldsymbol{D}^{-1} $$ 显然，这是一个上三角矩阵。□