3.3 QR分解

QR分解是将一个矩阵分解一个正交矩阵(酉矩阵)和一个三角矩阵的乘积. QR分解被广泛应用于线性最小二乘问题的求解和矩阵特征值的计算.

3.3.1 QR分解的存在性与唯一性

定理3.9 (QR分解) [73] 设 $A \in \mathbb{C}^{m \times n}$ ( $m \geq n$ ), 则存在一个单位列正交矩阵 $Q \in \mathbb{C}^{m \times n}$ (即 $Q^{*}Q = I_{n \times n}$ ) 和一个上三角矩阵 $R \in \mathbb{C}^{n \times n}$ , 使得

若 $A$ 列满秩, 则存在一个具有正对角线元素的上三角矩阵 $R$ 使得 (3.8) 成立, 且此时 QR 分解唯一, 即 $Q$ 和 $R$ 都唯一.

证明. 设 $A = [a_{1}, a_{2}, \ldots, a_{n}] \in \mathbb{C}^{m \times n}$ . 若 $A$ 列满秩, 即 $\operatorname{rank}(A) = n$ . 则 QR 分解 (3.8) 就是对 $A$ 的列向量组进行 Gram-Schmidt 正交化过程的矩阵描述 (见算法 3.3).

算法3.3.Gram-Schmidt正交化过程

1: $r_{11} = \left\| a_1\right\| _2$
2: $q_{1} = a_{1} / r_{11}$
3: for $j = 2$ to $n$ do
4: $q_{j} = a_{j}$
5: for $i = 1$ to $j - 1$ do
6: $r_{ij} = (a_j, q_i) \quad \%$ 计算内积, 用于正交化
7: $q_{j} = q_{j} - r_{ij}q_{i}$
8: end for
9: $r_{jj} = \| q_j\| _2$
10: $q_{j} = q_{j} / r_{jj}$
11: end for

由算法3.3可知

记 $Q = [q_{1}, q_{2}, \ldots, q_{n}]$ , $R = [r_{ij}]_{n \times n}$ , 其中

于是Gram-Schmidt正交化过程可表示为

如果 $A$ 不是列满秩, 我们可以通过下面的方式做类似的正交化过程:

• 如果 $a_1 = 0$ , 则令 $q_1 = 0$ ; 否则令 $q_1 = a_1 / \|a_1\|_2$ ;

• 对于 $j = 2,3,\ldots ,n$ ，计算 $\tilde{q}_j = a_j - \sum_{i = 1}^{j - 1}(q_i^* a_j)q_i$ 如果 $\tilde{q}_j = 0$ ，则表明 $a_{j}$ 可以由 $a_1,a_2,\dots ,a_{j - 1}$ 线性表出，令 $q_{j} = 0.$ 否则令 $q_{j} = \tilde{q}_{j} / \| \tilde{q}_{j}\|_{2}$ ：

于是我们有

其中 $Q = [q_{1}, q_{2}, \ldots, q_{n}]$ 列正交 (但不是单位列正交), 其列向量要么是单位向量, 要么就是零向量. $R = [r_{ij}]_{n \times n}$ 的定义同 (3.9). 需要注意的是, 如果 $Q$ 的某一列 $q_{k} = 0$ , 那么 $R$ 中对应的第 $k$ 行就全部为 0.

设 $\operatorname{rank}(A) = l < n$ ，则 $Q$ 有 $l$ 个非零列，设为 $q_{i_1}, q_{i_2}, \ldots, q_{i_l}$ . 它们形成 $\mathbb{C}^m$ 中的一个单位正交向量组，所以我们可以将其扩展成 $\mathbb{C}^m$ 中的一组标准正交基，即

然后我们用 $\tilde{q}_1$ 替换 $Q$ 中的第一个零列, 用 $\tilde{q}_2$ 替换 $Q$ 中的第二个零列, 依此类推, 将 $Q$ 中的所有零列都替换掉. 将最后得到的矩阵记为 $\tilde{Q}$ , 则 $\tilde{Q} \in \mathbb{C}^{m \times n}$ 单位列正交, 且

这是由于 $\tilde{Q}$ 中的新添加的列向量正好与 $R$ 中的零行相对应. 所以我们有 QR 分解

下面证明满秩矩阵QR分解的存在唯一性

存在性: 由于 $A$ 列满秩, 由 Gram-Schmidt 正交化过程 (算法 3.3) 可知, 存在上三角矩阵 $R = [r_{ij}]_{n \times n}$ 满足 $r_{jj} > 0$ , 使得 $A = QR$ , 其中 $Q$ 单位列正交.

唯一性: 假设 $A$ 存在 QR 分解

其中 $Q_{1}, Q_{2} \in \mathbb{C}^{m \times n}$ 单位列正交, $R_{1}, R_{2} \in \mathbb{C}^{n \times n}$ 为具有正对角元素的上三角矩阵. 则有

于是

又 $R_{1}, R_{2}$ 均为上三角矩阵, 所以 $R_{2}R_{1}^{-1}$ 也是上三角矩阵, 且其对角线元素为 $R_{2}(i,i) / R_{1}(i,i)$ ,

$i = 1,2,\dots ,n.$ 因此

同理可证 $R_{1}(i,i) / R_{2}(i,i)\leq 1.$ 所以

又 $\| Q_1\| _F^2 = \mathrm{tr}(Q_1^* Q_1) = n,$ 所以由(3.10)可知

由于 $R_{2}R_{1}^{-1}$ 的对角线元素都是1,所以 $R_{2}R_{1}^{-1}$ 只能是单位矩阵，即 $R_{2} = R_{1}$ 因此 $Q_{2} = AR_{2}^{-1} =$ $AR_1^{-1} = Q_1$ ，即 $A$ 的QR分解是唯一的. □

A 有时也将 QR 分解定义为: 存在酉矩阵 $Q \in \mathbb{C}^{m \times m}$ 使得

其中 $R = \left[ \begin{array}{c} R_{11} \\ 0 \end{array} \right] \in \mathbb{C}^{m \times n}$ 是上三角矩阵.

如果 $A$ 是实矩阵, 则上面证明中的运算都可以在实数下进行, 因此 $Q$ 和 $R$ 都可以是实矩阵.
如果 $A$ 是非奇异的方阵, 则 QR 分解也可以用来求解线性方程组 $Ax = b$ .

Gram-Schmidt正交化与正交投影

Gram-Schmidt 正交化过程的第 $j$ 步:

可以看作是 $a_{j}$ 减去其在 $\operatorname{span}\{q_1, q_2, \ldots, q_{j-1}\}$ 内的正交投影. 事实上, $r_{ij}q_i$ 就是 $a_{j}$ 在 $q_i$ 方向的正交投影, 即

下面给出QR分解的具体实现方法，分别基于MGS正交化过程, Householder变换和Givens变换.

3.3.2 基于MGS的QR分解

在证明QR分解的存在性时,我们利用了Gram-Schmidt正交化过程.但由于数值稳定性方面的原因,在实际计算中,我们一般不采用Gram-Schmidt正交化过程,取而代之的是修正的Gram-Schmidt正交化过程(modified Gram-Schmidt process,MGS),即对正交化过程做如下修改:

• Gram-Schmidt 正交化过程的第 $j$ 步:

(1) 计算 $r_{ij} = (a_j, q_i)$ , $i = 1, 2, \ldots, j - 1$ ;
(2) 计算 $\tilde{q}_j = a_j - r_{1j}q_1 - r_{2j}q_2 - \dots -r_{j - 1,j}q_{j - 1};$
(3) 计算 $r_{jj} = \|\tilde{q}_j\|$ , $q_j = \tilde{q}_j / r_{jj}$ ;

• MGS 正交化过程的第 $j$ 步:

(1) 令 $\tilde{q}_j = a_j$
(2) 计算 $r_{ij} = (\tilde{q}_j, q_i)$ , $\tilde{q}_j = \tilde{q}_j - r_{ij}q_i$ , $i = 1, 2, \ldots, j-1$ ;
(3) 计算 $r_{jj} = \|\tilde{q}_j\|$ , $q_j = \tilde{q}_j / r_{jj}$ ;

可以证明, 数学上这两个算法完全等价, 即 $r_{ij}$ 和 $q_j$ 都一样. 但在数值上, Gram-Schmidt 正交化过程是不稳定的, 而 MGS 是向后稳定的 [98].

算法3.4.基于MGS的QR分解

$\%$ Given $A \in \mathbb{C}^{m \times n}$ , compute $Q = [q_1, \ldots, q_n] \in \mathbb{R}^{m \times n}$ and $R \in \mathbb{R}^{n \times n}$ such that $A = QR$

1: Set $R = [r_{ij}] = 0_{n \times n}$ (the $n \times n$ zero matrix)
2: if $a_1 = 0$ then
3: $q_{1} = 0$
4: else
5: $r_{11} = \| a_1\| _2, \quad q_1 = a_1 / \| a_1\| _2$
6: end if
7: for $j = 2$ to $n$ do
8: $q_{j} = a_{j}$
9: for $i = 1$ to $j - 1$ do $\%$ MGS, 注意与 GS 的区别
10: $r_{ij} = (q_j,q_i),\quad q_j = q_j - r_{ij}q_i$
11: end for
12: if $q_{j}\neq 0$ then
13: $r_{jj} = \| q_j\| _2, \quad q_j = q_j / r_{jj}$
14: end if
15: end for

A 本算法的运算量大约为 $2mn^{2}$
算法中 $R$ 是一列一列计算的, 因此也称为按列MGS (column-oriented MGS). 按列MGS不适合列主元, 因此, 如果需要计算列主元QR时需采用按行MGS (row-oriented MGS), 可参见相关资料.
由MGS得到的QR分解中， $Q\in \mathbb{R}^{m\times n},R\in \mathbb{R}^{n\times n}$

3.3.3 基于Householder变换的QR分解

由定理3.5可知, 通过Householder变换, 我们可以将任何一个非零变量 $x \in \mathbb{R}^n$ 转化成 $\| x\|_2e_1$ , 即除第一个元素外, 其它都为零. 下面我们就考虑通过Householder变换来实现矩阵的QR分解.

设 $A \in \mathbb{R}^{m \times n}$ ( $m \geq n$ ), 构造 Householder 变换 $H_1 \in \mathbb{R}^{m \times m}$ , 使得

于是

其中 $\tilde{A}_2\in \mathbb{R}^{(m - 1)\times (n - 1)}$ .同样地,我们可以构造一个Householder变换 $\tilde{H}_2\in \mathbb{R}^{(m - 1)\times (m - 1)}$ ，将 $\tilde{A}_2$ 的第一列中除第一个元素外的所有元素都化为0，即

令 $H_{2} = \left[ \begin{array}{cc}1 & 0\\ 0 & \tilde{H}_{2} \end{array} \right]$ 则 $H_{2}\in \mathbb{R}^{m\times m}$ 也是Householder变换(见习题3.8)，且

不断重复上述过程, 我们就可以得到一系列的 Householder 变换

使得

由于Householder变换都是正交矩阵，令

则 $Q$ 也是正交矩阵，且

以上就是基于Householder变换的QR分解的具体实现过程.最后所得到的上三角矩阵 $R$ 就存放在 $A$ 的上三角部分.

如果不需要生成 $Q$ , 则基于 Householder 变换的 QR 分解的运算量大约为 $2mn^{2} - 2n^{3}/3$ .

矩阵 $Q$ 可通过下面的累积方法来计算：

如果保留了每一步的Householder向量, 则 $Q$ 也可以通过下面的向后累积方法来计算:

这样做的好处是一开始 $Q$ 会非常稀疏 (左上角为单位矩阵), 随着迭代的进行, $Q$ 才会慢慢变满. 而前面的计算方法, 第一步就将 $Q$ 变成了一个满矩阵.

采用向后累积方法计算 $Q$ 的运算量大约为 $4m^{2}n - 4mn^{2} + 4n^{3} / 3$ . 如果只需要计算 $Q$ 的前 $n$ 列, 则运算量大约为 $2mn^{2} - 2n^{3} / 3$ , 此时 QR 分解的总运算量为 $4mn^{2} - 4n^{3} / 3$ . 若 $m = n$ , 则为 $8n^{3} / 3$ .

算法3.5.基于Householder变换的QR分解

$\%$ Given $A \in \mathbb{R}^{m \times n}$ , compute $Q$ and $R$ such that $A = QR$ where $Q \in \mathbb{R}^{m \times m}$ and $R \in \mathbb{R}^{m \times n}$ .
$\%$ The upper triangular part of $R$ is stored in the upper triangular part of $A$

1: Set $Q = I_{m \times m}$
2: for $k = 1$ to $n$ do
3: $x = A(k:m,k)$
4: $[\beta_k, v_k] = \mathbf{House}(x)$
5: $A(k:m,k:n) = (I_{m - k + 1} - \beta_kv_kv_k^\top)A(k:m,k:n)$
6: $= A(k:m,k:n) - \beta_{k}v_{k}\big(v_{k}^{\mathsf{T}}A(k:m,k:n)\big)$
7: $Q(:, k : m) = Q(:, k : m)(I_{m-k+1} - \beta_k v_k v_k^\top)$
8: $= Q(:,k:m) - \beta_{k}\big(Q(:,k:m)v_{k}\big)v_{k}^{\mathsf{T}}$
9: end for

上面的算法只是关于Householder QR分解的一个简单描述，并没有考虑运算量问题.在实际计算时，我们通常会保留所有的Householder向量(用于向后累积法计算 $Q$ ).由于第 $k$ 步中 $\tilde{H}_k$ 所对应的Householder向量 $v_{k}$ 的长度为 $m - k + 1$ ，因此我们先把 $v_{k}$ 单位化，使得 $v_{k}$ 的第一元素为1,这样就只要存储 $v_{k}(2:end)$ ，共 $m - k$ 个元素.于是我们就可以把所有的

Householder 向量存放在 $A$ 的严格下三角部分. 而 $A$ 的上三角部分仍然存放 $R$ .

事实上, 如果没有明确要求生成 $Q$ , 通常只需存储 $v_{k}$ 和 $\beta_{k}$ , 即以因式分解 (Factored-Form) 形式存储 $Q$ [57].

我们也可以考虑块Householder QR分解，以便充分利用3级BLAS运算，提高计算效率.
算法3.5针对的是实矩阵，如果是复矩阵则要适当做些修改

3.3.4 列主元QR分解

如果 $A$ 不是满秩, 记 $l \triangleq \operatorname{rank}(A) < n$ , 则存在一个置换矩阵 $P$ , 使得 $AP$ 的前 $l$ 列是线性无关的. 因此我们可以对 $AP$ 进行 QR 分解, 于是我们可以得到下面的结论.

定理3.10 (列主元QR分解) 设 $A \in \mathbb{C}^{m \times n}$ ( $m \geq n$ ), 且 $\operatorname{rank}(A) = l < n$ . 则存在置换矩阵 $P$ 和正交矩阵 $Q \in \mathbb{C}^{m \times m}$ , 使得

其中 $R_{11}\in \mathbb{C}^{l\times l}$ 是非奇异上三角矩阵,且对角线元素满足 $r_{11}\geq r_{22}\geq \dots \geq r_{ll} > 0$

A 上述结论也可简化为

其中 $Q_{1}\in \mathbb{C}^{m\times l}$ 单位列正交(上述结论中 $Q$ 的前 $l$ 列).

列主元QR分解的实现过程与QR分解基本类似,只是在第 $k$ 步时，需要选列主元，同时可能需要做一个列交换.

假设经过 $k - 1$ 步后，我们得到下面的分解

其中 $P^{(k - 1)}$ 是置换矩阵, $Q^{(k - 1)}$ 是正交矩阵, $R_{11}^{(k - 1)}\in \mathbb{R}^{(k - 1)\times (k - 1)}$ 是非奇异上三角矩阵

下面考虑第 $k$ 步：

(1) 首先计算 $R_{22}^{(k - 1)}$ 的所有列的范数, 如果范数都为0, 则 $R_{22}^{(k - 1)} = 0$ , 此时必有 $k - 1 = l$ , 算法结束.
(2) 当 $k \leq l$ 时, $R_{22}^{(k-1)} \neq 0$ , 记其范数最大的列为第 $i_k$ 列 (如果有相等的, 取其中一个即可). 若 $i_k \neq 1$ , 则交换 $R^{(k-1)}$ 的第 $k$ 列与第 $i_k + k - 1$ 列, 并记相应的置换矩阵为 $P_k$ .
(3) 由于列交换不会影响到 $R^{(k - 1)}$ 的前 $k - 1$ 列, 因此列交换后的矩阵可记为

构造 $\tilde{R}_{22}^{(k - 1)}$ 的第1列所对应的Householder变换 $\tilde{H}_k$ ，并令 $H_{k} = \left[ \begin{array}{cc}I_{k - 1} & 0\\ 0 & \tilde{H}_{k} \end{array} \right],P^{(k)} = P^{(k - 1)}P_k,$ 于是

其中 $\tilde{H}_k\tilde{R}_{22}^{(k - 1)}$ 的第一列除第一个元素外, 其余都是零, 且该元素的值等于 $\tilde{R}_{22}^{(k - 1)}$ 的第1列的范数. 记 $Q^{(k)}\triangleq Q^{(k - 1)}H_k^{\mathsf{T}}$ , 则

其中 $R_{11}^{(k)}\in \mathbb{R}^{k\times k}$ 为非奇异上三角矩阵

依此类推, 直到第 $l$ 步, 我们就可以得到 $A$ 的列主元 QR 分解, 其中 $R_{11}$ 的对角线元素非负且按降序排列是由列主元的选取方法和 Householder 变换的性质得到的.

下面是列主元 QR 分解的一个应用.

推论3.11 (满秩分解) 设 $A \in \mathbb{C}^{m \times n}$ 且 $\operatorname{rank}(A) = r \leq \min\{m, n\}$ , 则存在满秩矩阵 $F \in \mathbb{C}^{m \times r}$ 和 $G \in \mathbb{C}^{r \times n}$ , 使得

3.3.5 基于Givens变换的QR分解

我们同样可以利用Givens变换来做QR分解

设 $A \in \mathbb{R}^{n \times n}$ , 首先构造一个Givens变换 $G_{21}$ , 作用在 $A$ 的最前面的两行上, 使得

由于 $G_{21}$ 只改变矩阵的第1行和第2行的值, 所以其它行保存不变. 然后再构造一个Givens变换 $G_{31}$ , 作用在 $G_{21}A$ 的第1行和第3行, 将其第一列的第三个元素化为零. 由于 $G_{31}$ 只改变矩阵的第1行和第3行的值, 所以第二行的零元素维持不变. 以此类推, 我们可以构造一系列的Givens变换 $G_{41}, G_{51}, \ldots, G_{n1}$ , 使得 $G_{n1} \cdots G_{21}A$ 的第一列中除第一个元素外, 其它元素都化为零, 即

下面我们可以对第二列进行类似的处理.构造Givens变换 $G_{32},G_{42},\ldots ,G_{n2}$ ，将第二列的第3至第 $_n$ 个元素全化为零，同时保持第一列不变

以此类推, 我们对其他列也做类似的处理. 最后, 通过构造 $\frac{1}{2} n(n - 1)$ 个Givens变换, 将 $A$ 转

化成一个上三角矩阵 $R$ ，即

令 $Q = (G_{n,n - 1}\dots G_{21})^{\mathsf{T}}$ .由于Givens变换是正交矩阵,所以 $Q$ 也是正交矩阵.于是，我们就得到矩阵 $A$ 的QR分解

与Householder变换一样,在进行Givens变换时，我们不需要显式地写出Givens矩阵.

对于稠密矩阵而言, 基于Givens变换的QR分解的运算量比Householder变换要多很多.
当需要连续应用一系列Givens变换时，我们可以使用快速Givens变换(见[16]).但即便如此,其速度仍然要慢于Householder变换.因此基于Givens变换的QR分解主要用于当矩阵的非零下三角元素相对较少时的情形(比如对上Hessenberg矩阵进行QR分解)，或者稀疏矩阵(可以尽可能地保留矩阵的稀疏性).
A 另外, 由于每次 Givens 变换只影响矩阵的两行或者两列, 因此非常适合并行计算.
如果 $A \in \mathbb{R}^{m \times n}$ , 其中 $m > n$ , 我们仍然可以通过Givens变换进行QR分解.

下面是基于Givens变换的QR分解的算法描述

算法3.6.基于Givens变换的QR分解

$\%$ Given $A \in \mathbb{R}^{m \times n}$ , compute $Q$ and $R$ such that $A = QR$ where $Q \in \mathbb{R}^{m \times m}$ and $R \in \mathbb{R}^{m \times n}$ .

$\%$ The upper triangular part of $R$ is stored in the upper triangular part of $A$

1: Set $Q = I_{m \times m}$
2: for $k = 1$ to $n$ do
3: for $i = k + 1$ to $m$ do
4: $[c, s] = \text{givens}(a_{kk}, a_{ik})$
5: $\begin{bmatrix} A(k,k:n)\\ A(i,k:n) \end{bmatrix} = G\begin{bmatrix} A(k,k:n)\\ A(i,k:n) \end{bmatrix}$ where $G = \begin{bmatrix} c & s\\ -s & c \end{bmatrix}$
6: $[Q(1:m,k),Q(1:m,i)] = [Q(1:m,k),Q(1:m,i)]G^{\mathsf{T}}$
7: end for
8: end for

3.3.6 QR分解的稳定性

基于Householder变换和Givens变换的QR分解都具有很好的数值稳定性.详细分析可以参考[135]和[70].基于MGS的QR分解也是向后稳定的,参见[98].

当需要计算矩阵 $Q$ 时, 基于MGS的QR分解的运算量相对较少. 因此当 $A$ 的列向量具有很好的线性无关性时, 我们可以使用MGS来计算QR分解. 但需要注意的是, MGS得到的 $Q$ 不是方阵, 除非 $A$ 是方阵.

但是, 由于舍入误差的原因, 最后得到的矩阵 $Q$ 会带有一定的误差, 可能会导致 $Q$ 失去正交性. Björck [15, 88] 证明了, 通过 MGS 计算的矩阵 $Q$ 满足

其中 $c_{1}(m,n)$ 是与 $m,n$ 相关的常数

而通过Householder变换计算的矩阵 $Q$ 满足

因此, 如果正交性至关重要, 则当 $A$ 的列向量接近线性相关时, 建议使用Householder变换.

另外, 基于MGS的QR分解所计算得到的 $Q$ (由于舍入误差的原因, $Q$ 不一定是列正交的)和 $R$ 满足[15,88]

而且存在单位列正交矩阵 $\tilde{Q}$ 使得 [57]

例3.3 编写程序, 分别用GS, MGS和Householder变换计算 $n$ 阶Hilbert矩阵 $H$ 的QR分解，并比较三种算法的稳定性，即观察 $\| \tilde{Q}\tilde{R} - H\|_2$ 和 $\|\tilde{Q}^{\mathsf{T}}\tilde{Q} - I\|_2$ 的值，其中 $\tilde{Q}$ 和 $\tilde{R}$ 是计算出来的QR分解矩阵因子. 试验结果如下. (QR_3methods.m)

也可以通过重正交来提升MGS的稳定性，即进行两次MGS.

例3.4 编写程序, 分别用 GS, MGS, 两次 MGS 和 Householder 变换计算 QR 分解. 下面是针对 256 的随机矩阵, 在不同条件数下各种方法的稳定性, 试验结果如下, 其中 MGS(double) 是指两次 MGS. (QR_3methods_MGS2.m, QR_3methods_256_64.m)