QR 分解

我们在介绍方程组的秩时， 描述了基和坐标的定义。上一节中，我们在此基础上介绍了线性变换， 最后介绍了可以将空间 \(\mathbb{V}\) 中的任意一组基转换成正交基的 Gram-Schmidt 算法。实际上利用该算法可以将一个矩阵 \(\boldsymbol{A}\) 分解成一个正交矩阵 \(\boldsymbol{Q}\) 和一个上三角矩阵\(\boldsymbol{R}\)，即所谓的 QR 分解 \(\boldsymbol{A} = \boldsymbol{Q}\boldsymbol{R}\)。如果矩阵 \(\boldsymbol{A}\) 可逆，那么该分解就是惟一的。

上一节最后提到的 Gram-Schmidt 正交化方法实际上就是一个 QR 分解，只是没有把上三角矩阵保存下来。 该方法直接从向量的正交基本概念出发，依次构造出标准正交基，简单直观，但是数值稳定性较差。工程上常用 Householder 变换， 它本身就是一种正交变换，并且与向量相乘时，有消元的作用。所以可以通过一系列的 Householder 变换将 A 逐步转化为上三角阵。该方法数值稳定。

1. 基于 Gram-Schmidt 的分解

对于一个 \(n \times n\) 的矩阵 \(\boldsymbol{Q}\)，如果它的转置就是它的逆，既 \(\boldsymbol{Q}^T = \boldsymbol{Q}^{-1}\)，那么该矩阵就是一个正交矩阵。 如果我们将正交矩阵看作是一组列向量，那么每个列向量的长度都是 1，并且两两正交。此外，两个正交矩阵的乘积仍然是正交矩阵，正交矩阵左乘到一个列向量上并不会改变向量的长度。 Gram-Schmidt 算法正是基于这些性质逐步构建出正交矩阵 \(\boldsymbol{Q}\) 和上三角矩阵 \(\boldsymbol{R}\) 的。

假设我们有一个 \(n \times n\) 的矩阵 \(\boldsymbol{A}\)，写成列向量的形式为 \(\boldsymbol{A} = \begin{bmatrix} \boldsymbol{\alpha}_1 & \boldsymbol{\alpha}_2 & \boldsymbol{\alpha}_3 & \cdots & \boldsymbol{\alpha}_n \end{bmatrix}\)。对其进行 QR 分解之后得到的正交矩阵 \(\boldsymbol{Q}\)， 写成列向量的形式为 \(\boldsymbol{Q} = \begin{bmatrix} \boldsymbol{\gamma}_1 & \boldsymbol{\gamma}_2 & \boldsymbol{\gamma}_3 & \cdots & \boldsymbol{\gamma}_n \end{bmatrix}\)。有

$$ \begin{bmatrix} \boldsymbol{\alpha}_1 & \boldsymbol{\alpha}_2 & \boldsymbol{\alpha}_3 & \cdots & \boldsymbol{\alpha}_n \end{bmatrix} = \begin{bmatrix} \boldsymbol{\gamma}_1 & \boldsymbol{\gamma}_2 & \boldsymbol{\gamma}_3 & \cdots & \boldsymbol{\gamma}_n \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & r_{13} & \cdots & r_{1n} \\ & r_{22} & r_{23} & \cdots & r_{2n} \\ & & r_{33} & \cdots & r_{3n} \\ & & & \ddots & \vdots \\ & & & & r_{nn} \end{bmatrix} $$

令 \(i = 1\)，有 \(\boldsymbol{\alpha}_1 = r_{11}\boldsymbol{\gamma}_1\)。显然 \(r_{11} = \| \boldsymbol{\alpha}_1 \|\)，\(\boldsymbol{\gamma}_1 = \boldsymbol{\alpha}_1 / \| \boldsymbol{\alpha}_1 \|\)。

令 \(i = 2\)，有 \(\boldsymbol{\alpha}_2 = r_{12}\boldsymbol{\gamma}_1 + r_{22}\boldsymbol{\gamma}_2\)。因为 \(\boldsymbol{Q}\) 是个正交矩阵， 所以有\(\|\boldsymbol{\gamma}_1\| = 1, \|\boldsymbol{\gamma}_2\| = 1, \boldsymbol{\gamma}_1 \cdot \boldsymbol{\gamma}_2 = 0\)。因此有:

$$ \boldsymbol{\gamma}_1 \cdot \boldsymbol{\alpha}_2 = r_{12} \boldsymbol{\gamma}_1 \cdot \boldsymbol{\gamma}_1 + r_{22} \boldsymbol{\gamma}_1 \cdot \boldsymbol{\gamma}_2 = r_{12} $$

根据上式可以得出 \(r_{12}\)，所以有 \(r_{22}\boldsymbol{\gamma}_2 = \boldsymbol{\alpha}_2 - r_{12}\boldsymbol{\gamma}_1\)，即， \(r_{22} = \| \boldsymbol{\alpha}_2 - r_{12}\boldsymbol{\gamma}_1\|, \boldsymbol{\gamma}_2 = (\boldsymbol{\alpha}_2 - r_{12}\boldsymbol{\gamma}_1) / r_{22}\)。

令 \(i = 3\)，有 \(\boldsymbol{\alpha}_3 = r_{13}\boldsymbol{\gamma}_1 + r_{23}\boldsymbol{\gamma}_2 + r_{33}\boldsymbol{\gamma}_3\)。利用正交矩阵的性质，可以计算得到 \(r_{13}, r_{23}\)

$$ \boldsymbol{\gamma}_1 \cdot \boldsymbol{\alpha}_3 = r_{13} \boldsymbol{\gamma}_1 \cdot \boldsymbol{\gamma}_1 + r_{23} \boldsymbol{\gamma}_1 \cdot \boldsymbol{\gamma}_2 + r_{33} \boldsymbol{\gamma}_1 \cdot \boldsymbol{\gamma}_3 = r_{13} $$ $$ \boldsymbol{\gamma}_2 \cdot \boldsymbol{\alpha}_3 = r_{13} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_1 + r_{23} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_2 + r_{33} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_3 = r_{23} $$

根据上式，有 \(r_{33}\boldsymbol{\gamma}_3 = \boldsymbol{\alpha}_3 - r_{13} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_1 - r_{23} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_2\)。 那么 \(r_{33} = \|\boldsymbol{\alpha}_3 - r_{13} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_1 - r_{23} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_2\|\), \(\boldsymbol{\gamma}_3 = (\boldsymbol{\alpha}_3 - r_{13} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_1 - r_{23} \boldsymbol{\gamma}_2 \cdot \boldsymbol{\gamma}_2) / r_{33}\)。

实际上上述过程正是我们在上一节最后提到的 Gram-Schmidt 正交化方法。 对于 \(2 < k ≤ n\)，我们都可以通过 \(\boldsymbol{\alpha}_k - \sum_{i=2}^{k} \frac{\boldsymbol{\alpha}_k \cdot \boldsymbol{\beta}_i}{\boldsymbol{\beta}_i\boldsymbol{\beta}_i}\boldsymbol{\beta}_i\) 构造出 \(\boldsymbol{\beta}_k\) 与 \(\boldsymbol{\beta}_1 \cdots \boldsymbol{\beta}_{k-1}\) 都正交。 如果我们再对 \(M_2\) 中的每个向量单位化，即，令 \(\boldsymbol{\gamma}_i = \boldsymbol{\beta}_i / |\boldsymbol{\beta}_i|\)， 那么就可以得到一组单位正交基 \(M_3 = \left\{ \boldsymbol{\gamma}_1, \cdots, \boldsymbol{\gamma}_n\right\}\)。这正式 QR 分解中的正交矩阵。 上三角矩阵中对角线上的元素 \(r_{ii}\) 则是向量 \(\boldsymbol{\beta}_i\) 的模长。其余元素 \(r_{ij} = \boldsymbol{\gamma}_i \cdot \boldsymbol{\alpha}_j\)。 实现参见例程QR_GramSchmidt。

Gram-Schmidt 算法十分简洁，但不幸的是它的数值稳定性不怎么样，在工程上很少使用。基于 Householder 变换的 QR 分解算法，能在计算过程中较好地控制舍入误差的传播， 非常适合用于大规模的矩阵分解。

2. Householder 矩阵

Householder 变换是一种正交变换，它通过构造一个镜像反射矩阵，也就是所谓的 Householder 矩阵 \(\boldsymbol{H}\)，在保持向量长度不变的同时， 将一个向量 \(\boldsymbol{x}\)映射到另一个指定方向的向量 \(\boldsymbol{y} = \boldsymbol{H}\boldsymbol{x}\)。通过合理的选择 Householder 向量 \(\boldsymbol{v}\)， 可以将原向量 \(\boldsymbol{x}\) 中的指定元素变成 0。

给定一个非零向量 \(\boldsymbol{v}\)，可以按照如下的方式构造 Householder 矩阵: $$ \boldsymbol{H} = \boldsymbol{I} - \frac{2 \boldsymbol{v}\boldsymbol{v}^T}{\| \boldsymbol{v} \|^2} $$ 其中，\(\boldsymbol{I}\) 是一个单位矩阵，\(\boldsymbol{v}\boldsymbol{v}^T\) 是向量 \(\boldsymbol{v}\) 的外积，\(\| \boldsymbol{v} \|\)则是 \(\boldsymbol{v}\) 的长度。 矩阵 \(\boldsymbol{H}\) 是正交的，即 \(\boldsymbol{H}^T = \boldsymbol{H}^{-1}\)；也是对称的，即 \(\boldsymbol{H}^T = \boldsymbol{H}\)；也是对合的， 即\(\boldsymbol{H}^{-1} = \boldsymbol{H}\)。Householder 矩阵的几何意义是，将向量 \(\boldsymbol{x}\) 关于一个垂直于 \(\boldsymbol{v}\) 的超平面 \(\text{span}(\boldsymbol{v})^{\perp}\) 的镜面反射\(\boldsymbol{H}\boldsymbol{x}\)。

由于 \(\boldsymbol{H}\) 本身就是正交的，所以我们可以不断的左乘不同的 Householder 矩阵， 将目标矩阵 \(\boldsymbol{A}\) 上三角化，即，\(\boldsymbol{A} = \boldsymbol{H}_k \cdots \boldsymbol{H}_2 \boldsymbol{H}_1 \boldsymbol{R}\)， 其中 \(\boldsymbol{H}_k \cdots \boldsymbol{H}_2 \boldsymbol{H}_1\) 是一系列的 Householder 矩阵，也就是 QR 分解中的那个矩阵 \(\boldsymbol{Q}\)，\(\boldsymbol{R}\) 则是那个上三角矩阵。

假设 \(\boldsymbol{x}, \boldsymbol{y}\) 是空间 \(\mathbb{F}^n\) 中两个非零向量，如果它们的长度相等，即 \(\|\boldsymbol{x}\| = \| \boldsymbol{y} \|\)， 那么就一定存在一个矩阵 \(\boldsymbol{H}\)，使得 \(\boldsymbol{y} = \boldsymbol{H}\boldsymbol{x}\)。因为两个向量的模长相等， 所以有 \(\boldsymbol{x}^T\boldsymbol{x} = \boldsymbol{y}^T\boldsymbol{y}\)，并且 \(\boldsymbol{x}^T\boldsymbol{y} = \boldsymbol{y}^T\boldsymbol{x}\)。因此，存在如下关系: $$ \| \boldsymbol{x} - \boldsymbol{y} \|^2 = (\boldsymbol{x} - \boldsymbol{y})^T(\boldsymbol{x} - \boldsymbol{y}) = \boldsymbol{x}^T\boldsymbol{x} + \boldsymbol{y}^T\boldsymbol{y} - (\boldsymbol{x}^T\boldsymbol{y} + \boldsymbol{y}^T\boldsymbol{x}) = 2(\boldsymbol{x}^T\boldsymbol{x} - \boldsymbol{y}^T\boldsymbol{x}) $$ 令 \(\boldsymbol{v} = \boldsymbol{x} - \boldsymbol{y}\)，根据 Householder 矩阵的定义有 $$ \boldsymbol{H}(\boldsymbol{v}) \boldsymbol{x} = \boldsymbol{x} - \frac{2 (\boldsymbol{x} - \boldsymbol{y})(\boldsymbol{x} - \boldsymbol{y})^T \boldsymbol{x}}{\| \boldsymbol{x} - \boldsymbol{y} \|^2} = \boldsymbol{x} - \frac{2 (\boldsymbol{x} - \boldsymbol{y})(\boldsymbol{x}^\boldsymbol{x} - \boldsymbol{y}^T\boldsymbol{x})} {2(\boldsymbol{x}^T\boldsymbol{x} - \boldsymbol{y}^T\boldsymbol{x})} = \boldsymbol{y} $$ 所以存在 Householder 矩阵 \(\boldsymbol{H}(\boldsymbol{v})\) 使得 \(\boldsymbol{y} = \boldsymbol{H}(\boldsymbol{v}) \boldsymbol{x}\)。 向量 \(\boldsymbol{v}\) 被称为 Householder 向量。假设 \(\boldsymbol{x} = \begin{bmatrix} x_1 & \cdots & x_n \end{bmatrix}^T\)， \(\boldsymbol{y} = \begin{bmatrix} \alpha & \cdots & 0 \end{bmatrix}\)，其中 \(\alpha = \pm \| \boldsymbol{x} \|\)。 我们可以直接构造出 Householder 向量 \(\boldsymbol{v} = \begin{bmatrix} x_1 - \alpha & \cdots & x_n \end{bmatrix}\)， 使得经 Householder 矩阵映射后的向量除第一个元素外的所有元素都变成零。这一特性可以拿来对矩阵进行 QR 分解。在工程应用中，为了控制舍入误差，我们应当尽量避免两个相近的数做减法运算， 所以，通常取 \(\alpha = -\text{sign}(x_1) \| \boldsymbol{x} \|\)。 我们在单元测试 TEST(QR, HouseholderMatrix) 中给出了一种只保留指定元素的实现。

3. 基于 Householder 的分解

基于 Householder 变换的 QR 分解的基本思想就是不断的构造 Householder 矩阵对目标矩阵进行消元，只保留上三角区域的元素。假设有一个 \(m \times n\) 的矩阵 \(\boldsymbol{A}\)， 对其进行 QR 分解之后的正交矩阵是 \(m \times m\) 的，上三角矩阵则是 \(m \times n\) 的。

我们将 \(\boldsymbol{A} = \begin{bmatrix} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \cdots & a_{mn} \end{bmatrix}\)，写成列向量的形式， \(\boldsymbol{A} = \begin{bmatrix} \boldsymbol{\alpha}_1 & \cdots & \boldsymbol{\alpha}_n \end{bmatrix}\)，其中 \(\boldsymbol{\alpha}_i = \begin{bmatrix} a_{1i} \\ \vdots \\ a_{mi} \end{bmatrix}\)。对于 \(\boldsymbol{\alpha}_1\) 有矩阵 \(\boldsymbol{H}_1\)使得 \(\boldsymbol{H}_1 \boldsymbol{\alpha}_1 = \begin{bmatrix} r_{11} \\ \vdots \\ 0 \end{bmatrix}\)。那么 $$ \boldsymbol{H}_1 \boldsymbol{A} = \left[ \begin{array}{c|c} r_{11} & \begin{matrix} r_{12} & \cdots & r_{1n} \end{matrix} \\ \hline \begin{matrix} 0 \\ \vdots \\ 0 \end{matrix} & \tilde{\boldsymbol{A}}_2 \end{array}\right] $$ 其中 \(\tilde{\boldsymbol{A}}_2\) 是一个 \((m-1) \times (n-1)\) 的矩阵，同样的方式，我们可以构造 \((m-1) \times (m-1)\) 的矩阵 \(\tilde{\boldsymbol{H}}_2\) 对 \(\tilde{\boldsymbol{A}}_2\) 消元仅保留第一列第一行的元素，即 $$ \tilde{\boldsymbol{H}}_2 \tilde{\boldsymbol{A}}_2 = \left[ \begin{array}{c|c} r_{22} & \begin{matrix} r_{23} & \cdots & r_{2n} \end{matrix} \\ \hline \begin{matrix} 0 \\ \vdots \\ 0 \end{matrix} & \tilde{\boldsymbol{A}}_3 \end{array}\right] $$ 构造 \(\boldsymbol{H}_2 = \begin{bmatrix} 1 & \boldsymbol{0} \\ \boldsymbol{0} & \tilde{\boldsymbol{H}}_2 \end{bmatrix}\)，那么有 $$ \boldsymbol{H}_2 \boldsymbol{H}_1 \boldsymbol{A} = \left[ \begin{array}{c|c} \begin{matrix} r_{11} & r_{12} \\ 0 & r_{22} \end{matrix} & \begin{matrix} r_{13} & \cdots & r_{1n} \\ r_{23} & \cdots & r_{2n} \end{matrix} \\ \hline \begin{matrix} 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{matrix} & \tilde{\boldsymbol{A}}_3 \end{array}\right] $$ 不断重复上述过程，我们需要构造出 \(m - 1\) 个 Householder 矩阵，才能最终将矩阵\(\boldsymbol{A}\) 上三角化，即 \(\boldsymbol{H}_{m-1} \cdots \boldsymbol{H}_2 \boldsymbol{H}_1\boldsymbol{A} = \boldsymbol{R} = \begin{bmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\ & r_{22} & \cdots & r_{2n} \\ & & \ddots & \vdots \\ & & & r_{mn} \end{bmatrix}\)。

令 \(\boldsymbol{Q} = (\boldsymbol{H}_{m-1} \cdots \boldsymbol{H}_2\boldsymbol{H}_1)^{-1} = \boldsymbol{H}_1\boldsymbol{H}_2 \cdots \boldsymbol{H}_{m-1}\)， 则有对矩阵 \(\boldsymbol{A}\) 的 QR 分解 \(\boldsymbol{A} = \boldsymbol{Q}\boldsymbol{R}\)。 实现参见例程QR_Householder。

4. 完