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

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

      什么是QR分解
      QR分解是一种矩阵分解。QR的分解公式如下,
\begin{eqnarray}
M = Q R
\end{eqnarray}
其中M是 m \times n 的矩阵,并且m \ge n, rank(M) = n。Q是 m \times m 的正交单位矩阵,R是上三角矩阵。

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

\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] 列向量的一组正交基:
\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}
转换一下
\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}

      根据上面的公式,易推出
\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 。迭代公式如下:
\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列向量后面的元素。
\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}
如果要知道HouseHolder算法中Q的值,可以利用以下公式计算。
\begin{eqnarray}
Q = \prod_{k=1}^{m}(\pmb{I} - \beta_k \pmb{u}_k\pmb{u}_k^{T}) \nonumber
\end{eqnarray}

      让我们举个例子, 假设我们要求解的矩阵为 R_0
\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
\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
\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 ,然后交替进行如下步骤
\begin{eqnarray}
1.& M_k &=& QR (QR分解) \nonumber \\
2.& M_{k+1} &=& RQ \nonumber
\end{eqnarray}
因为 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||似乎应该在分子上而不在分母上?

发表评论

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