QR分解

QR分解法是一种将矩阵分解的方式。这种方式，把矩阵分解成一个正交矩阵与一个上三角矩阵的积。QR分解经常用来解线性最小二乘法问题。QR分解也是特定特征值算法即QR算法的基础。

## 类别及定义

### 方阵

任何方块矩阵A都可以分解为

$$A = QR$$

其中Q是正交矩阵（意味着QTQ = I）而R是上三角矩阵。如果A是非奇异的，且限定R的对角线元素为正，则这个因数分解是唯一的。

更一般的说，我们可以因数分解复数$m$×$n$矩阵（有着m ≥ n）为$m$×$n$幺正矩阵（在Q ∗Q = I 的意义上，不需要是方阵）和$n$×$n$上三角矩阵的乘积。对m1是n×n上三角矩阵，0是(m − n)×n零矩阵，Q1是m×n，Q2是m×(m − n)，且Q1和Q2都是有正交列。

```
call ''Q''<sub>1</sub>''R''<sub>1</sub> the ''thin QR factorization'' of ''A''; Trefethen and Bau call this the ''reduced QR factorization''.<ref name="Trefethen"></ref> If ''A'' is of full [[matrix rank|rank]] ''n'' and we require that the diagonal elements of ''R''<sub>1</sub> are positive then ''R''<sub>1</sub> and ''Q''<sub>1</sub> are unique, but in general ''Q''<sub>2</sub> is not. ''R''<sub>1</sub> is then equal to the upper triangular factor of the [[Cholesky decomposition]] of ''A''<span style="vertical-align: 10%; margin-left: 0.1em">*</span> ''A'' (=&nbsp;''A''<sup>T</sup>''A'' if ''A'' is real).
```

### QL、RQ 和 LQ 分解

类似的，我们可以定义A的QL，RQ和LQ分解。其中L是下三角矩阵。

## QR分解的求法

QR分解的实际计算有很多方法，例如Givens旋转、Householder变换，以及Gram-Schmidt正交化等等。每一种方法都有其优点和不足。

### 使用Householder变换

#### Householder变换

Householder变换将一个向量关于某个平面或者超平面进行反射。我们可以利用这个操作对$m \times n ( m \geqq n)$的矩阵$A$进行QR分解。

矩阵$Q$可以被用于对一个向量以一种特定的方式进行反射变换，使得它除了一个维度以外的其他所有分量都化为0。

令$\mathbf{x}$为矩阵$A$的任一m维实列向量，且有$\|\mathbf{x}\| = |\alpha|$（其中$\alpha$为标量）。若该算法是通过浮点数实现的，则$\alpha$应当取和$\mathbf{x}$的第$k$维相反的符号（其中$x_k$是要保留不为0的项），这样做可以避免精度缺失。对于复数的情况，令

$$\alpha = - \mathrm{e}^{\mathrm{i} \arg x_k} \|\mathbf{x}\|$$

Stoer，并且在接下来矩阵$Q$的构造中要将矩阵转置替换为共轭转置。

接下来，设$\mathbf{e}_1$为单位向量$(1, 0, \cdots, 0)^T$。

·

为欧几里德范数，$I$为$m \times m$单位矩阵，令

$\mathbf{u} = \mathbf{x} - \alpha\mathbf{e}_1$。

$\mathbf{v} = {\mathbf{u}\over\|\mathbf{u}\|}$。

$Q = I - 2 \mathbf{v}\mathbf{v}^T$。

或者，若$A$为复矩阵，则

$Q = I - (1+w)\mathbf{v}\mathbf{v}^H$，其中$w = \mathbf{x}^H\mathbf{v}\mathbf{/}\mathbf{v}^H\mathbf{x}$

式中$\mathbf{x}^H$是$x$的共轭转置（亦称埃尔米特共轭或埃尔米特转置）。

则$Q$为一个$m \times m$的Householder矩阵，它满足

$$Q\mathbf{x} = (\alpha, 0, \cdots, 0)^T \ $$

利用Householder矩阵，可以将一个$m \times n$的矩阵$A'$变换为上三角矩阵。

首先，我们将A左乘通过选取矩阵的第一列得到列向量$x$的Householder矩阵$Q_1$。这样，我们得到的矩阵$Q_1 A$的第一列将全部为0（第一行除外）：

$$Q_1A = \begin{bmatrix} <pre> \alpha_1&\star&\dots&\star\\ 0 & & & \\ \vdots & & A' & \\ 0 & & & \end{bmatrix}$$

这个过程对于矩阵$A'$（即$Q_1 A$排除第一行和第一列之后剩下的方阵）还可以继续做下去，从而得到另一个Householder矩阵$Q_2$。注意到$Q_2$其实比$Q_1$要小，因为它是在$Q_1 A$而非$A$的基础上得到的。因此，我们需要在$Q_2$的左上角补上1，或者，更一般地来说：

$$Q_k = \begin{bmatrix} <pre> I_{k-1} & 0\\ 0 & Q_k'\end{bmatrix}$$

将这个迭代过程进行$t$次之后（$t = \min(m-1, n)$）,将有

$$R = Q_t \cdots Q_2Q_1A$$

其中R为一个上三角矩阵。因此，令

$$Q = Q_1^T Q_2^T \cdots Q_t^T,$$

则$A = QR$为矩阵$A$的一个QR分解。

相比与Gram-Schmidt正交化，使用Householder变换具有更好的数值稳定性。

#### 例子

现在要用Householder变换求解矩阵$A$的$QR$ 分解。

$$A = <pre> \begin{bmatrix} 0 & 3 & 1 \\ 0 & 4 & -2 \\ 2 & 1 & 1 \\ \end{bmatrix} </pre>$$

因为$\alpha_1 = [ 0,\ 0,\ 2 ]^T$, 令$a_1 = ||\alpha_1||_2 = 2$，则

$$\omega_1 = \frac{\alpha_1 - a_1e_1}{||\alpha_1 - a_1e_1||_2} = \frac{1}{\sqrt{2}}[ -1,\ 0,\ 1 ]^T$$

则有

$$H_1 = I - 2\omega_1\omega_1^H= <pre> \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \\ \end{bmatrix}$$

从而。

$$H_1A = <pre> \begin{bmatrix} 2 & 1 & 1 \\ 0 & 4 & -2 \\ 0 & 3 & 1 \\ \end{bmatrix}$$

记$\beta = [4,\ 3]^T$, 则$b_1 = ||\beta_2||_2 = 5$。令

$$\omega_2 = \frac{\beta_2 - b_1e_1}{||\beta_2 - b_1e_1||_2} = \frac{1}{\sqrt{10}}[ -1,\ 3 ]^T$$

$$\hat{H_2} = I - 2\omega_2\omega^H = \frac{1}{5} <pre> \begin{bmatrix} 4 & 3 \\ 3 & -4 \\ \end{bmatrix}$$

记。

$$H_2= <pre> \begin{bmatrix} 1 & 0^T \\ 0 & \hat{H_2} \\ \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \frac{4}{5} & \frac{3}{5} \\ 0 & \frac{3}{5} & -\frac{4}{5} \\ \end{bmatrix}$$

则。

$$R = H_2(H_1A) = <pre> \begin{bmatrix} 2 & 1 & 1 \\ 0 & 5 & -1 \\ 0 & 0 & -2 \\ \end{bmatrix}$$

那么

$$Q = H_1H_2 = \frac{1}{5} <pre> \begin{bmatrix} 0 & 3 & -4 \\ 0 & 4 & 3 \\ 5 & 0 & 0 \\ \end{bmatrix} </pre>$$

### 使用吉文斯旋转

#### 吉文斯旋转

吉文斯旋转表示为如下形式的矩阵

$$G(i, j, \theta) = <pre> \begin{bmatrix} 1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\ 0 & \cdots & c & \cdots & -s & \cdots & 0 \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ 0 & \cdots & s & \cdots & c & \cdots & 0 \\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1 \end{bmatrix}$$

这里的 c = cos(θ) 和 s = sin(θ) 出现在第 i 行和第 j 行与第 i 列和第 j 列的交叉点上。就是说，吉文斯旋转矩阵的所有非零元定义如下：:

$$\begin{align} <pre> g_{k\, k} &{}= 1 \qquad \text{for} \ k \ne i,\,j\\ g_{i\, i} &{}= c \\ g_{j\, j} &{}= c \\ g_{i\, j} &{}= s \\ g_{j\, i} &{}= -s </pre> \end{align}$$

乘积 G(i, j, θ)x 表示向量 x 在 (i,j)平面中的逆时针旋转 θ 弧度。

#### 吉文斯旋转作用于QR分解

对于一个向量

$$\begin{array}{lcl} A &=& \begin{bmatrix} <pre> a\\ b\\ </pre> \end{bmatrix} \\ \end{array}$$

如果，$r = \sqrt{a^2 + b^2}$ , $c = \frac{a}{r}$ , $s = - \frac{b}{r}$ ,

那么，就存在旋转矩阵G，使$A$ 底部转成0。

$$A_{2\_Sub} = <pre> \begin{bmatrix} c & -s \\ s & c \\ \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix} </pre> = \begin{bmatrix} r \\ <pre> 0 \\ \end{bmatrix} </pre>$$

每一次的旋转，吉文斯旋转都可以将一个元素化成0，直到将原始矩阵转成一个上三角矩阵，则完成分解。

$$A = QR$$

$$Q = G_1^TG_2^T \cdots G_k^T$$

#### 例子

$$A_1 = <pre> \begin{bmatrix} 6 & 5 & 0 \\ 5 & 1 & 4 \\ 0 & 4 & 3 \\ \end{bmatrix}$$

$$r = \sqrt{6 ^2 + 5^2} \approx 7.8102$$

$$c = 6 / r \approx 0.7682$$

$$s = -5 / r \approx -0.6402$$

$$A_2 = G_1A_1 = <pre> \begin{bmatrix} c & -s & 0 \\ s & c & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} 6 & 5 & 0 \\ 5 & 1 & 4 \\ 0 & 4 & 3 \\ \end{bmatrix} \approx \begin{bmatrix} 7.8102 & 4.4813 & 2.5607 \\ 0 & -2.4327& 3.0729 \\ 0 & 4 & 3 \\ \end{bmatrix} </pre>$$

对于:$A_2$子矩阵 :$A_{2\_Sub}$

$$A_{2\_Sub} = <pre> \begin{bmatrix} -2.4327& 3.0729 \\ 4 & 3 \\ \end{bmatrix} </pre>$$

$$r = \sqrt{(-2.4327)^2 + 4^2} \approx 4.6817$$

$$c = -2.4327 / r \approx -0.5196$$

$$s = - 5 / r \approx -0.8544$$

$$G_2A_2 = <pre> \begin{bmatrix} 1 & 0 & 0 \\ 0 & c & -s \\ 0 & s & c \\ \end{bmatrix} \begin{bmatrix} 7.8102 & 4.4813 & 2.5607 \\ 0 & -2.4327& 3.0729 \\ 0 & 4 & 3 \\ \end{bmatrix} \approx \begin{bmatrix} 7.8102 & 4.4813 & 2.5607 \\ 0 & 4.6817 & 0.9664 \\ 0 & 0 & -4.1843 \\ \end{bmatrix} </pre>$$

$$R = G_2A_2 = G_2 G_1 A_1$$

$$Q = G_1^T G_2^T = <pre> \begin{bmatrix} 0.7682 & 0.3327 & 0.5470 \\ 0.6402 & -0.3992 & -0.6564 \\ 0 & 0.8544 & -0.5196 \\ \end{bmatrix} </pre>$$

### 使用格拉姆-施密特正交化方法

#### 基本思想

<image>

格拉姆-施密特正交化的基本想法，是利用投影原理在已有正交基的基础上构造一个新的正交基。

设$\boldsymbol{v} \in \boldsymbol{V^n}$。$\boldsymbol{V}^k$是$\boldsymbol{V}^n$上的$k$维子空间，其标准正交基为$\{ \boldsymbol{\eta}_1,\ldots, \boldsymbol{\eta}_k \}$，且$\boldsymbol{v}$不在$\boldsymbol{V}^k$上。由投影原理知，$\boldsymbol{v}$与其在$\boldsymbol{V}^k$上的投影$\mathrm{proj}_{\boldsymbol{V^k}} \boldsymbol{v}$之差

$$<pre> \boldsymbol{\beta} = \boldsymbol{v} - \sum_{i=1}^{k}\mathrm{proj}_{\boldsymbol{\eta}_i}\,\boldsymbol{v} = \boldsymbol{v} - \sum_{i=1}^{k}\langle \boldsymbol{v}, \boldsymbol{\eta}_i \rangle \boldsymbol{\eta}_i$$

是正交于子空间$\boldsymbol{V}^k$的，亦即$\boldsymbol{\beta}$正交于$\boldsymbol{V}^k$的正交基$\boldsymbol{\eta}_i$。因此只要将$\boldsymbol{\beta}$单位化，即

$$\boldsymbol{\eta}_{k+1} = \frac{\boldsymbol{\beta}}{\|\boldsymbol{\beta}\|} = \frac{\boldsymbol{\beta}}{\sqrt{\langle \boldsymbol{\beta},\boldsymbol{\beta} \rangle }}$$

那么$\{ \boldsymbol{\eta}_1,\ldots, \boldsymbol{\eta}_{k}, \boldsymbol{\eta}_{k+1} \}$就是$\boldsymbol{V}^k$在$\boldsymbol{v}$上扩展的子空间$\mathrm{span}\{\boldsymbol{v},\boldsymbol{\eta}_1,...,\boldsymbol{\eta}_k\}$的标准正交基。

根据上述分析，对于向量组$\{ \boldsymbol{v}_1,\ldots, \boldsymbol{v}_{m} \}$张成的空间$\boldsymbol{V}^m$ ($m<n$)，只要从其中一个向量（不妨设为$\boldsymbol{v}_1$）所张成的一维子空间$\mathrm{span}\{\boldsymbol{v}_1\}$开始（注意到$\boldsymbol{v}_1$就是$\mathrm{span}\{\boldsymbol{v}_1\}$的正交基），重复上述扩展构造正交基的过程，就能够得到$\boldsymbol{V}^n$ 的一组正交基。这就是格拉姆-施密特正交化。

#### 格拉姆-施密特正交化算法

首先需要确定已有基底向量的顺序，不妨设为$\{ \boldsymbol{v}_1,\ldots, \boldsymbol{v}_{n} \}$。Gram-Schmidt正交化的过程如下：

$\boldsymbol{\beta}_1 = \boldsymbol{v}_1,$｜$\boldsymbol{\eta}_1 = {\boldsymbol{\beta}_1 \over \|\boldsymbol{\beta}_1\|}$

$\boldsymbol{\beta}_2 = \boldsymbol{v}_2-\langle \boldsymbol{v}_2, \boldsymbol{\eta}_1 \rangle \boldsymbol{\eta}_1,$｜$\boldsymbol{\eta}_2 = {\boldsymbol{\beta}_2 \over \|\boldsymbol{\beta}_2\|}$

$\boldsymbol{\beta}_3 = \boldsymbol{v}_3 - \langle \boldsymbol{v}_3, \boldsymbol{\eta}_1 \rangle \boldsymbol{\eta}_1 - \langle \boldsymbol{v}_3, \boldsymbol{\eta}_2 \rangle \boldsymbol{\eta}_2 ,$｜$\boldsymbol{\eta}_3 = {\boldsymbol{\beta}_3 \over \|\boldsymbol{\beta}_3\|}$

$\vdots$｜$\vdots$

$\boldsymbol{\beta}_n = \boldsymbol{v}_n-\sum_{i=1}^{n-1}\langle \boldsymbol{v}_n, \boldsymbol{\eta}_i \rangle \boldsymbol{\eta}_i,$｜$\boldsymbol{\eta}_n = {\boldsymbol{\beta}_n\over\|\boldsymbol{\beta}_n\|}$

这样就得到$\mathrm{span}\{ \boldsymbol{v}_1, \ldots , \boldsymbol{v}_n \}$上的一组正交基$\{ \boldsymbol{\beta}_1, \ldots , \boldsymbol{\beta}_n \}$，以及相应的标准正交基$\{ \boldsymbol{\eta}_1, \ldots , \boldsymbol{\eta}_n \}$。

#### 例子

现在要用格拉姆-施密特变换求解矩阵$A$的$QR$ 分解。

$$A = <pre> \begin{bmatrix} 1 & 2 & 4 \\ 0 & 0 & 5 \\ 0 & 3 & 6 \\ \end{bmatrix} </pre>$$

令， $a = [1 , 0 , 0]$

$$q_1 = \frac{a}{||a||} = <pre> \begin{bmatrix} 1 \\ 0 \\ 0 \\ \end{bmatrix} </pre>$$

$$\hat{q_2} = b - (b * q_1)q_1 = <pre> \begin{bmatrix} 2 \\ 0 \\ 3 \\ \end{bmatrix} - 2 \begin{bmatrix} 1 \\ 0 \\ 0 \\ \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 3 \\ \end{bmatrix} </pre>$$

$$q_2 = \frac{\hat{q_2}}{||\hat{q_2}||} = <pre> \begin{bmatrix} 0 \\ 0 \\ 1 \\ \end{bmatrix} </pre>$$

$$\hat{q_3} = c - (c * q_1)q_1 - (c * q_2)q_2 = <pre> \begin{bmatrix} 4 \\ 5 \\ 6 \\ \end{bmatrix} - 4 \begin{bmatrix} 1 \\ 0 \\ 0 \\ \end{bmatrix} - 6 \begin{bmatrix} 0 \\ 0 \\ 1 \\ \end{bmatrix} </pre> = <pre> \begin{bmatrix} 0 \\ 5 \\ 0 \\ \end{bmatrix} </pre>$$

$$q_3 = \frac{\hat{q_3}}{||\hat{q_3}||} = <pre> \begin{bmatrix} 0 \\ 1 \\ 0 \\ \end{bmatrix} </pre>$$

那么可知。

$$Q = <pre> \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \end{bmatrix} </pre>$$

由$A = QR$，可知。

$$R = <pre> \begin{bmatrix} 1 & 2 & 4 \\ 0 & 3 & 6 \\ 0 & 0 & 5 \\ \end{bmatrix} </pre>$$

## Matlab

MATLAB以qr函数来执行QR分解法，其语法为

$$[Q,R]=qr(A)$$

其中Q代表正规正交矩阵。

而R代表上三角形矩阵。

此外，原矩阵A不必为正方矩阵。

如果矩阵A大小为$m \times n$，则矩阵Q大小为$m \times m$，矩阵R大小为$m \times n$。

## 用途

### 解线性方程组

对于直接求解线性方程组的逆，用QR分解的方法求解会更具有数据的稳定性。

对于求解一个线性系统$Ax = b$, 这里$A$的维度是$m \times n$。

如果$m \leq n$, 那么$A^{T} =QR$,这里$Q^{T} = Q^{-1}$)。

$R$ 的形式是 $R = \begin{bmatrix} R_1 \\ 0 \end{bmatrix}$，$R_1$是$R$上不为0的部分。

那么对于

$$x = Q \begin{bmatrix} <pre> \left(R_1^\textsf{T}\right)^{-1}b \\ 0 \end{bmatrix} </pre>$$

如果$m > n$, 那么$A =QR$,这里$Q^{T} = Q^{-1}$)。本质是最小化$||A\hat{x} - b||$

$$\hat{x} = R_1^{-1} \left(Q_1^\textsf{T} b\right)$$