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>$$

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

#### 基本思想

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

設$\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)$$