CS205A Lecture 5 Column spaces and QR

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第五讲，这一讲主要介绍了$QR$分解。

Chapter 4 列空间和$QR$分解

正规方程的结构

上一讲主要介绍了对称矩阵的一些性质，现在考虑$A^TA$的条件数：

$$\begin{aligned} \text{cond }A^T A &=||A^T A||.||(A^T A)^{-1}||\\ &\approx ||A^T||.||A||.||A^{-1}||.||(A^T)^{-1}||\\ &=||A||^2 ||A^{-1}||^2\\ &=(\text{cond }A)^2 \end{aligned}$$

这说明如下方程更加不稳定：

$$A^TA \vec x =\vec b$$

注意到

$$\text{cond } I_{n} =1$$

所以如果我们有

$$A^T A=I_n$$

那么就不会有不稳定的情形发生，这就引入了正交性的概念。

正交性

对于矩阵

$$Q= \left( \begin{matrix} \mid & \mid & & \mid \\ \vec q_1 &\vec q_2 & \ldots & \vec q_n \\ \mid & \mid & & \mid \end{matrix} \right)$$

我们计算$Q^T Q$

$$Q^TQ= \left( \begin{matrix} — (q_1)^T— \\ — (q_2)^T— \\ \vdots\\ — (q_n)^T— \end{matrix} \right)\left( \begin{matrix} \mid & \mid & & \mid \\ \vec q_1 &\vec q_2 & \ldots & \vec q_n \\ \mid & \mid & & \mid \end{matrix} \right) =\left( \begin{matrix} \vec q_1.\vec q_1 &\vec q_1.\vec q_2 &\ldots& \vec q_1.\vec q_n \\ \vec q_2.\vec q_1 &\vec q_2.\vec q_2 &\ldots &\vec q_2.\vec q_n \\ \vdots &\vdots &\ldots &\vdots\\ \vec q_n.\vec q_1 &\vec q_n.\vec q_2 &\ldots& \vec q_n.\vec q_n \end{matrix} \right)$$

如果

$$Q^T Q=I_n$$

那么

$$\vec q_i .\vec q_j=\begin{cases} 1 & i=j\\ 0 &i\neq j \end{cases}$$

我们给出如下定义：

定义 4.1 (正交性，正交矩阵)

一组向量$\{\vec v_1,…,\vec v_k \}$正交当且仅当

$$\vec v_i .\vec v_j=\begin{cases} 1 & i=j\\ 0 &i\neq j \end{cases}$$

如果一个方阵的列向量正交，那么称其为正交矩阵。

性质

正交矩阵有如下一些性质：

1. $$Q^{-1}=Q^T$$

   这一点是显然的，由逆矩阵的定义即可得到。这个性质告诉我们正交矩阵的逆可以直接得到，大大减少了运算量。

2. $$(Q\vec x ).(Q\vec y ) =\vec x^T Q^T Q \vec y =\vec x^T \vec y =\vec x .\vec y$$

   这个性质表明正交矩阵保持内积不变。

3. 对性质2中取$\vec x =\vec y$，那么

   $$||Q\vec x ||^2 =\vec x^T Q^T Q \vec x =\vec x^T \vec x =||\vec x ||^2$$

   这个性质表明正交矩阵保持模长不变。

由性质2，3，我们知道$Q$对应于几何变换，即保持长度和角度不变。

对非正交矩阵的处理

首先回顾我们一直讨论的问题：

$$\min_{\vec x} ||A\vec x -\vec b||$$

从几何角度来看，可以理解为上述问题是在将$\vec b$投影到$A$的列空间，所以这一部分主要从列空间的角度讨论，我们先给出如下引理：

引理 4.1 (列空间的不变性)

对于任意$A\in \mathbb R^{m\times n}$和可逆矩阵$B\in \mathbb R^{n\times n}$，我们有

$$\text{col }A =\text{col }AB$$

证明：$\forall \vec b \in \text{col }A​$，存在$\vec x​$，使得

$$\vec b = A\vec x$$

因为$B​$可逆，所以

$$\vec b = (AB) (B^{-1}\vec x) \in \text{col }AB$$

所以

$$\text{col }A \subset \text{col }AB$$

$\forall \vec c \in \text{col }AB​$，存在$\vec y​$，使得

$$\vec c = (AB)\vec y =A(B\vec y )\in \text{col }A$$

所以

$$\text{col }AB \subset \text{col }A$$

因此

$$\text{col }A =\text{col }AB$$

上述性质提醒我们可以使用初等列变换对$A$进行操作，直至其变成正交矩阵，即

$$Q= AE_1E_2...E_k$$

如果$E_i$可逆，那么上述引理告诉我们

$$\text{col }Q =\text{col }A$$

那么

$$A=QR ,R=E^{-1}_k E^{-1}_{k-1}... E^{-1}_{1}$$

如果得到上述分解，考虑正规方程

$$A^TA\vec x = A^T\vec b$$

那么我们有

$$\begin{aligned} \vec x &= (A^T A)^{-1} A^T\vec b\\ &=(R^T Q^T Q R)^{-1}R^T Q^T \vec b\\ &=(R^TR)^{-1}R^T Q^T \vec b\\ &=R^{-1}(R^T)^{-1}R^T Q^T\vec b \\ &=R^{-1}Q^T\vec b \end{aligned}$$

或者等价的

$$R\vec x = Q^T \vec b$$

如果我们设计$R$为三角阵，那么求解$R\vec x = Q^T \vec b$只是在执行back-substitution。关于$QR$分解，常见的有两种方法，分别在后面介绍。

Gram-Schmidt正交化

介绍Gram-Schmidt正交化之前，首先考虑投影问题。

投影

投影问题的形式如下：对于向量$\vec a ,\vec b$，下式何时最小

$$\min_c ||c\vec a -\vec b||^2$$

这是$A\vec x =\vec b$的特殊形式，我们知道上述问题的解为

$$c = \frac{\vec a^T \vec b}{\vec a^T \vec a} =\frac{\vec a .\vec b}{||\vec a||^2}$$

我们称如下量为$\vec b $在$\vec a$上的投影：

$$\text{proj}_{\vec a} \vec b =c\vec a =\frac{\vec a .\vec b}{||\vec a||^2} \vec a$$

接着，我们考虑残差$\vec b - \text{proj}_{\vec a} \vec b$和$\vec a ​$的关系：

$$\begin{aligned} \vec a . (\vec b - \text{proj}_{\vec a} \vec b) &=\vec a .\vec b - \vec a .\Big(\frac{\vec a .\vec b}{||\vec a||^2} \vec a \Big)\\ &=\vec a .\vec b- \frac{\vec a .\vec b}{||\vec a||^2} (\vec a .\vec a )\\ &=\vec a .\vec b-\vec a .\vec b\\ &=0 \end{aligned}$$

上述事实告诉我们，我们将$\vec b$分解为平行于$\vec a $的部分和正交于$\vec a​$的部分。

现在，假设向量$\hat a_1 ,\hat a_2,…,\hat a_k​$正交（$||\hat a_i ||=1​$），那么对于每个$i​$，我们有

$$\text{proj}_{\vec a_i} \vec b =({\hat a_i .\vec b}) \hat a_i$$

进一步，我们可以将$\vec b$投影到$\{\hat a_1 ,\hat a_2,…,\hat a_k\}$张成的子空间，通过对下式关于$c_1,…,c_k \in \mathbb R$最小化，

$$\begin{aligned} || c_1 \hat a_1 +c_2 \hat a_2+...+c_k \hat a_k -\vec b||^2 &=\Big(\sum_{i=1}^k\sum_{j=1}^k c_i c_j(\hat a_i.\hat a_j) \Big) -2\vec b .\Big(\sum_{i=1}^k c_i \hat a_i \Big)+\vec b.\vec b\\ &= \sum_{i=1}^k (c_i^2 -2c_i \vec b. \hat a_i )+||\vec b||^2 \end{aligned}$$

最后一步是因为

$$\vec a_i .\vec a_j=\begin{cases} 1 & i=j\\ 0 &i\neq j \end{cases}$$

对上式关于$c_i$求偏导，不难得到

$$c_i =\hat a_i .\vec b$$

所以

$$\text{proj}_{\text{span}\{\hat a_1 ,\hat a_2,...,\hat a_k\}} \vec b = \sum_{i=1}^k({\hat a_i .\vec b}) \hat a_i$$

并且不难得到

$$\hat a_i .(\vec b -\text{proj}_{\text{span}\{\hat a_1 ,\hat a_2,...,\hat a_k\}} \vec b)=0$$

将上述内容整理一下即得到Gram-Schmidt正交化方法：

输入：线性独立的向量$\{\vec v_1 ,…\vec v_k\}​$

输入：正交基$\{\hat v_1 ,…\hat v_k\}$

1. 令

   $$\hat a_1 =\frac{\vec v_1}{||\vec v_1||}$$

2. 对$i=2…k$

   1. 计算

      $$\vec p_i =\text{proj}_{\text{span}\{\hat a_1 ,\hat a_2,...,\hat a_{i-1}\}} \vec v_i$$
   $$\hat a_i =\frac{\vec v_i -\vec p_i }{||\vec v_i -\vec p_i ||}$$

不难发现，对所有的$i$，我们有

$${\text{span}\{\hat v_1 ,\hat v_2,...,\hat v_i\}}={\text{span}\{\hat a_1 ,\hat a_2,...,\hat a_i\}}$$

Gram-Schmidt正交化可以解决很多问题，但有时候会产生数值不稳定性，考虑如下列子

$$\vec v_1 = \left( \begin{matrix} 1\\ 1 \end{matrix} \right), \vec v_2=\left( \begin{matrix} 1+\epsilon\\ 1 \end{matrix} \right)$$

此时显然有正交基

$$(1,0),(0,1)$$

但是如果使用Gram-Schmidt正交化，则有

$$\begin{aligned} \hat a_1 & =\frac{\vec v_1}{||\vec v_1||} =\frac 1 {\sqrt 2} \left( \begin{matrix} 1\\ 1 \end{matrix} \right) \\ \vec p_2 &=\frac {2+\epsilon}{2}\left( \begin{matrix} 1\\ 1 \end{matrix} \right) \\ \vec v_2 -\vec p_2 &= \left( \begin{matrix} 1+\epsilon\\ 1 \end{matrix} \right) -\frac {2+\epsilon}{2}\left( \begin{matrix} 1\\ 1 \end{matrix} \right) \\ &=\frac 1 2 \left( \begin{matrix} \epsilon\\ - \epsilon \end{matrix} \right) \end{aligned}$$

注意到

$$|| \vec v_2 -\vec p_2|| =\frac {\sqrt 2 } 2 \epsilon$$

如果$\epsilon $很小，那么除以上述量会产生计算问题。

Householder变换

Householder变换采取高斯消元法的思路，我们不断左乘正交矩阵$Q_i$，最终将$A$变为上三角矩阵，具体形式如下

$$Q_k...Q_1 A= R$$

所以

$$A= Q_1^T ...Q_k^T R =Q R\\ Q= Q_1^T ...Q_k^T$$

注意正交矩阵有无数种，Householder变换选择特殊的正交矩阵——反射矩阵，这是因为反射矩阵可以很方便的用投影表示出来。事实上，$2\text{proj}_{\vec v} \vec b -\vec b$为$\vec b$关于$\vec v $的反射向量，并且

$$\begin{aligned} 2\text{proj}_{\vec v} \vec b -\vec b &= 2\frac{\vec v .\vec b}{\vec v .\vec v}\vec v - \vec b \\ &=2\vec v.\frac{\vec v^T\vec b}{\vec v^T\vec v} - \vec b \\ &=\Big(\frac{2\vec v\vec v^T}{\vec v^T\vec v}-I_n \Big)\vec b\\ &\triangleq -H_{\vec v}\vec b \end{aligned}$$

所以反射可以看成对$\vec b​$做线性变换，事实上，我们还有

$$\begin{aligned} H_{\vec v}^T H_{\vec v} &=\Big(I_n -\frac{2\vec v\vec v^T}{\vec v^T\vec v}\Big)^T\Big(I_n -\frac{2\vec v\vec v^T}{\vec v^T\vec v} \Big)\\ &=\Big(I_n -\frac{2\vec v\vec v^T}{\vec v^T\vec v}\Big)\Big(I_n -\frac{2\vec v\vec v^T}{\vec v^T\vec v}\Big)\\ &=I_n+\frac{2\vec v\vec v^T}{\vec v^T\vec v}\frac{2\vec v\vec v^T}{\vec v^T\vec v}-\frac{2\vec v\vec v^T}{\vec v^T\vec v}-\frac{2\vec v\vec v^T}{\vec v^T\vec v} \\ &=I_n+\frac{4(\vec v\vec v^T)(\vec v^T\vec v)}{(\vec v^T\vec v)(\vec v^T\vec v)} -\frac{4\vec v\vec v^T}{\vec v^T\vec v} \\ &=I_n+\frac{4\vec v\vec v^T}{\vec v^T\vec v}-\frac{4\vec v\vec v^T}{\vec v^T\vec v}\\ &=I_n \end{aligned}$$

从而$H_{\vec v}$为正交矩阵，我们利用这一点来介绍Householder变换。

假设我们在做高斯消元法forward substitution的第一步，记矩阵$A$的第一列为$\vec a $，第一个单位向向量为$\vec e_1$，那么第一步消元的效果为除了第一行第一列元素以外都变成$0​$，即

$$H_{\vec v} A= \left( \begin{matrix} c & \times & \times &\times \\ 0 & \times & \times &\times \\ \vdots & \vdots & \vdots &\vdots \\ 0 & \times & \times &\times \end{matrix} \right)$$

所以

$$\begin{aligned} c\vec e_1 &= H_{\vec v} \vec a \\ &=\Big(I_n-\frac{2\vec v\vec v^T}{\vec v^T\vec v} \Big)\vec a \\ &=\vec a - 2\vec v \frac{\vec v^T\vec a}{\vec v^T\vec v} \end{aligned}$$

移项可得

$$\vec v =(\vec a -c\vec e_{1}) .\frac{\vec v^T \vec v}{2 \vec v^T \vec a }$$

所以$\vec v$平行于$\vec a -c\vec e_1$，因为$\vec v$的模长不改变$H_{\vec v}$，所以不妨选择

$$\vec v =a-c\vec e_1$$

从而

$$\begin{aligned} 1 &=\frac{\vec v^T \vec v}{2 \vec v^T \vec a } \\ &=\frac{(a-c\vec e_1)^T (a-c\vec e_1)}{2 (a-c\vec e_1)^T \vec a}\\ &=\frac{||\vec a||^2 -2c\vec e_1 .\vec a +c^2}{2(\vec a .\vec a -c\vec e_1 .\vec a )} \end{aligned}$$

即

$$\begin{aligned} ||\vec a||^2 -2c\vec e_1 .\vec a +c^2 &=2(\vec a .\vec a -c\vec e_1 .\vec a ) \\ c^2 &=||\vec a ||^2 \\ c&=\pm ||\vec a || \end{aligned}$$

现在将上述步骤推广，假设进行到第$k$步，记$A$的第$k​$列为

$$\vec a = \left( \begin{matrix} \vec a_1 \\ \vec a_{2} \end{matrix} \right)$$

其中$\vec a _1 \in \mathbb R^{k-1},\vec a _2 \in \mathbb R^{n-k+1}$，$a_0$是实数，消元以后的结果为

$$H_{\vec v} \vec a = \left( \begin{matrix} \vec a_1 \\ c \\ 0 \end{matrix} \right)$$

所以

$$\begin{aligned} \Big(I_n-\frac{2\vec v\vec v^T}{\vec v^T\vec v} \Big)\vec a &= \left( \begin{matrix} \vec a_1 \\ c \\0 \end{matrix} \right) \\ \vec a -2\vec v \frac{\vec v^T \vec a}{\vec v^T\vec v} &=\left( \begin{matrix} \vec a_1 \\ c \\0 \end{matrix} \right) \\ 2\vec v \frac{\vec v^T\vec a}{\vec v^T\vec v} & = \vec a -\left( \begin{matrix} \vec a_1 \\ c \\0 \end{matrix} \right)\\ &=\left( \begin{matrix} \vec a_1 \\ \vec a_{2} \end{matrix} \right) -\left( \begin{matrix} \vec a_1 \\ c \\0 \end{matrix} \right) \\ &= \left( \begin{matrix} 0\\ \vec a_{2} \end{matrix} \right) -c e_k \end{aligned}$$

同之前讨论，我们取

$$\vec v = \left( \begin{matrix} 0 \\ \vec a_{2} \end{matrix} \right) -c e_k$$

所以

$$\begin{aligned} \frac{2\vec v^T\vec a}{\vec v^T\vec v} &=1 \\ 2\vec v^T\vec a&=\vec v^T\vec v \\ 2\Big(\left( \begin{matrix} 0 \\ \vec a_{2} \end{matrix} \right) -c e_k\Big)^T \left( \begin{matrix} \vec a_1 \\ \vec a_2 \end{matrix} \right) &=\Big(\left( \begin{matrix} 0 \\ \vec a_{2} \end{matrix} \right) -c e_k\Big)^T\Big(\left( \begin{matrix} 0 \\ \vec a_{2} \end{matrix} \right) -c e_k\Big)\\ 2(|| \vec a_2||^2-c e_k^T \vec a )&= ||\vec a_2||^2 +c^2 -2c e_k^T \vec a \\ c^2 &=||\vec a _2 ||^2 \end{aligned}$$

所以最终结果为

$$\vec v = \left( \begin{matrix} 0 \\ \vec a_{2} \end{matrix} \right) -c e_k \\ c=\pm ||\vec a_2||$$

约化的$QR$分解

之前讨论的操作都是在方阵上，对于非方阵$A\in \mathbb R^{m\times n }$也可以进行同样的操作，但是Gram-Schmidt正交化和Householder变换的结果有一些不同：

• 如果使用Gram-Schmidt正交化，因为使用的是列变换，所以最终结果为

  $$Q\in \mathbb R^{m\times n}, R\in \mathbb R^{n\times n}$$

• 如果使用Householder变换，因为是使用正交矩阵进行消元，所以最终结果为

  $$Q\in \mathbb R^{m\times m}, R\in \mathbb R^{m\times n}$$

假设现在有$m\gg n$，因为数值稳定性，我们更倾向于使用Householder变换，但是$m\times m$的矩阵可能非常大，幸运的是，我们知道$R$为上三角矩阵，所以

$$R= \left( \begin{matrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & 0& \times \\ 0&0&0\\ 0&0&0 \end{matrix} \right)$$

因此

$$A=QR= \left( \begin{matrix} Q_1 & Q_2 \end{matrix} \right)\left( \begin{matrix} R_1 \\ 0 \end{matrix} \right)=Q_1 R_1$$

这里

$$Q_1\in \mathbb R^{m\times n}, R_1\in \mathbb R^{n\times n}$$

这种形式被称为$A$的约化的$QR$分解。