一种古老的矩阵分解:QR分解

      最近的研究涉及到了矩阵的QR分解和Rank Revealing QR分解,因此看了一些相关文献。本来想把QR分解和RRQR分解的相关的内容都整理下的。但是后面发现,RRQR实在是太复杂和琐碎,一时半会讲不清楚,只能把QR分解相关内容整理下了。

      什么是QR分解
      QR分解是一种矩阵分解。QR的分解公式如下,

(1)   \begin{eqnarray*} M = Q R \end{eqnarray*}

其中M是m \times n的矩阵,并且m\gen, rank(M) = n。Q是m \times m的正交单位矩阵,R是上三角矩阵。

      QR分解的物理含义: 找到矩阵M的列空间的一组单位正交基(Q的列向量),使得M的前k个列向量,恰好在前k个基向量所形成的子空间中。下面是QR分解的一个例子。

(2)   \begin{eqnarray*} \underset{M}{ \begin{bmatrix} 1 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 1  \end{bmatrix}} =  \underset{Q}{ \begin{bmatrix} -0.7071 & 0  &  0.7071 \\ -0.7071 & 0  & -0.7071 \\ 0       & 1 & 0  \end{bmatrix}} \underset{R}{ \begin{bmatrix} -1.414 & -1.414 & -0.7071 \\ 0      & 1      & 1 \\ 0      & 0      & -0.7071  \end{bmatrix}} \nonumber \\ \end{eqnarray*}

分别用Q_jM_j表示Q和M的第j列,我们发现M_1 = -1.414 Q_1, M_2= -1.414 Q_1 - Q_2, M_3 = -0.7074 Q_1 - Q_2 -0.7071Q_3,验证了M的前k个列向量,能用Q的前k个列向量表示,即M的前k个列向量恰好在前k个基向量所形成的子空间中。

QR算法
      我们有两种算法求解一个m \times n(m \ge n)的列满秩矩阵M的QR分解: 1)Gram-Schmidt(格拉姆-施密特正交法), 2)HouseHolder Transformation算法。其中Gram-Schmidt给出了QR分解存在性的构造性证明。

      一. Gram-Schmidt算法
      运用Gram-Schmidt正交法, 求m \times n(m \ge n)的列满秩矩阵M = [M_1,...,M_n]列向量的一组正交基:

(3)   \begin{eqnarray*} O_1 &=& M_1 \nonumber \\ ... \nonumber\\ O_k &=& M_{k} - \frac{M_{k}^TO_1}{O_1^tO_1}O_1 - ... - \frac{M_{k}^TO_{k-1}}{O_{k-1}^tO_{k-1}}O_{k-1} \nonumber \\ ... \nonumber  \end{eqnarray*}

转换一下

(4)   \begin{eqnarray*} M_1 &=& O_1 \nonumber \\ ... \nonumber\\ M_k &=& O_{k} + \frac{M_{k}^TO_1}{O_1^tO_1}O_1 + ... + \frac{M_{k}^TO_{k-1}}{O_{k-1}^tO_{k-1}}O_{k-1} \nonumber \\ ... \nonumber  \end{eqnarray*}

      根据上面的公式,易推出

(5)   \begin{eqnarray*} [M_1,M_2,...] &=& [O_1,O_2,...]\begin{bmatrix} 1 &  \frac{M_{2}^TO_1}{O_1^tO_1} & \frac{M_{3}^TO_1}{O_1^tO_1} & ...\\ 0 &  1                           & \frac{M_{3}^TO_2}{O_1^tO_2} & ...\\ 0 &  0                           &  1                          & ...\\ ...  &  ...                      &  ...                        & ...\end{bmatrix} \nonumber \\  &=&  \underset{\boldsymbol{Q}}{[\frac{O_1}{||O_1||},\frac{O_2}{||O_2||},...]}       \underset{\boldsymbol{R}}{\begin{bmatrix} 1 &  \frac{M_{2}^TO_1}{O_1^tO_1||O_1||} & \frac{M_{3}^TO_1}{O_1^tO_1||O_1||} & ...\\ 0 &  1                           & \frac{M_{3}^TO_2}{O_1^tO_2||O_2||} & ...\\ 0 &  0                           &  1                          & ...\\ ...  &  ...                      &  ...                        & ...\end{bmatrix}} \nonumber  \end{eqnarray*}

      由于对于n个线性不相关的向量[M_1,...,M_n] (即M列满秩),Gram-Schmidt正交法总能找到一组正交基[O_1,O_2,...], 因此M必定存在QR分解。Gram-Schmidt给出了QR分解存在性的构造性证明。

      二, HouseHolder Transformation算法

      HouseHolder算法是迭代算法。算法从R_0=M开始, 迭代地产生R_1,R_2,...,R_kR_n就是要求的上三角矩阵R。迭代公式如下:

(6)   \begin{eqnarray*} R_{k+1} = (\pmb{I} - \beta_k \pmb{u}_k\pmb{u}_k^{T})R_k = R_k - \beta_k \pmb{u}_k\pmb{u}_k^{T}R_k \end{eqnarray*}

其中\beta_k\pmb{u}的计算公式如下所示。\pmb{u}是m元向量,其中第1到第k-1元素为0,而第k+1到最后一个元素是R_{k}的第k列向量后面的元素。

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
\sigma_k = (\sum_{i=k}^{m} a_{i,k}^2)^{1/2} &\quad& \beta_k  = \frac{1}{\sigma_k*(\sigma_k + |a_{k,k}|)} \nonumber \\
\pmb{u}_k[i]    = 0   (k \lt i \le m) \quad 
\pmb{u}_k[k]    &=& sign(a_{k,k})(\sigma_k + |a_{k,k}|) \quad
\pmb{u}_k[i]    = a_{i,k}   (1 \le i \lt k) \nonumber 
\end{eqnarray*}

*** Error message:
Undefined control sequence \lt.
leading text: \pmb{u}_k[i]    = 0   (k \lt

如果要知道HouseHolder算法中Q的值,可以利用以下公式计算。

(8)   \begin{eqnarray*} Q = \prod_{k=1}^{m}(\pmb{I} - \beta_k \pmb{u}_k\pmb{u}_k^{T}) \nonumber \end{eqnarray*}

      让我们举个例子, 假设我们要求解的矩阵为R_0

(9)   \begin{eqnarray*} R_0 = \begin{bmatrix} 1 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 1  \end{bmatrix} \nonumber \\ \end{eqnarray*}

根据上面的公式,我们可以计算出\beta_0,\pmb{u}_0,继而计算出R_1

(10)   \begin{eqnarray*} \beta_0 = 0.292893218  \quad \pmb{u}_0 = \begin{bmatrix} 2.414 \\ 1\\ 0 \end{bmatrix}  \quad R_1 = (\pmb{I} - \beta_0 \pmb{u}_0\pmb{u}_0^{T})R_0 = \begin{bmatrix} -1.414 & -1.414 & -0.707 \\ 0      &  0     & 0.707  \\ 0      &  1     & 1       \end{bmatrix} \nonumber \end{eqnarray*}

继续迭代,可以计算出R_2R_3

(11)   \begin{eqnarray*} R_2 = \begin{bmatrix} -1.414 & -1.414 & -0.707 \\ 0      &  1     & 1  \\ 0      &  0     & 0.7071    \end{bmatrix} \quad  R_3 = \begin{bmatrix} -1.414 & -1.414 & -0.707 \\ 0      &  1     & 1  \\ 0      &  0     & -0.7071  \end{bmatrix}  \nonumber\\ \end{eqnarray*}

最终的结果为R = R_3

      这个例子可以体现HouseHolder算法的如下两个性质。这两个性质的证明太繁琐,就不列上了。
      1,算法的第k次迭代,在保证R_k的前k列向量不变的情况下,将R_k的第k+1列向量转变成,只有前k个元素不为0的列向量。即R_{k+1}的前k+1列满足上三角矩阵的形式。
      2,算法的第k次迭代之后, |R_{k+1}[k,k]| = \sigma_k

      QR分解的应用
      说了那么久的QR分解,那QR分解有什么用处呢?QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交相似变化成为Hessenberg矩阵,然后再应用QR方法求特征值和特征向量。运用QR分解求特征值的算法流程如下:令M_1 = M,然后交替进行如下步骤

*** QuickLaTeX cannot compile formula:
\begin{eqnarray*}
1.& M_k     &=& QR (QR分解) \nonumber \\
2.& M_{k+1} &=& RQ \nonumber
\end{eqnarray*}

*** Error message:
Too many columns in eqnarray environment.
leading text: 2

因为M_{k+1} = RQ = Q^TM_kQ,所以M_{k+1}M_{k}有相同的特征值。随着上述两个步骤交替进行, M的特征值就在M_k对角线上。

此条目发表在算法荟萃分类目录。将固定链接加入收藏夹。

一种古老的矩阵分解:QR分解》有 4 条评论

  1. 就冲这打Tex公式的精神也要赞一下!

  2. 匿名说:

    在讲GS算法的时候,在“根据上面的公式,易推出…”下面的推到的第二个等号右边的第二个矩阵,它的对角线元素应该是 || O_i || ?对角线上方||O_i||似乎应该在分子上而不在分母上?

发表评论

电子邮件地址不会被公开。