[數值計算] QR分解

本文轉載自查看原文 2021-06-02 11:03 340 數據算法

https://zhuanlan.zhihu.com/p/84415000

0. 為什么要用QR分解

$Ax=b,\ A\in\mathbb{R}^{m\times n}$ 的問題可以分成3類：

情況1：A是方陣，m=n
情況2：A是over-determined的，m>n
情況3：A是under-determined的，m<n

在[數值計算] 條件數的例子2里，遇到的情況1（A是方陣），通過構造拉格朗日插值來使得對A求逆足夠穩定。對於一般的情況下，解決思路是使用LU（LUP）分解來解決穩定性問題，在前一篇文中已經簡介過了[數值計算] LU分解、LUP分解、Cholesky分解。

對於后兩種情況， [數值計算] 數據擬合——線性最小二乘法分析了用正規方程組求解over-determined以及under-determined的問題。但在文中也提到了，對於over-determined的線性最小二乘問題，正規方程組是不穩定的，通常需要用QR分解來處理：

理論很美好，在小數據量的時候沒問題，然而直接使用正規方程組求解會在數據量大(e.g. data size > 100)的時候不穩定numerically unstable。原因是需要對 $A^TA$ 求逆，而A我們都知道是Vandermonde矩陣的一部分，本身就是poorly conditioned，而 $A^TA$ 只會更糟糕。解決的方法是使用QR分解，這也是Python MATLAB求解線性最小二乘問題的方法。

1. QR分解

1.1 定義

一個矩陣 $A \in \mathbb{R}^{m\times n},\ m\ge n$ 可以被分解成 $A = QR$ ，其中：

$Q\in \mathbb{R}^{m\times m}$ 是正交矩陣
$R \equiv \begin{bmatrix} \hat{R} \\ 0 \end{bmatrix} \in \mathbb{R}^{m\times n}$
$\hat{R} \in \mathbb{R}^{n\times n}$ 是上三角矩陣

1.2 正交矩陣的性質

$Q^T Q = QQ^T=I$
左乘一個正交矩陣對歐式范數的結果不影響（在下面證明eq.2的時候會用到）

$|| Qv ||_2^2 = v^TQ^TQv = v^Tv = || v ||_2^2 \tag{1}$

1.3 從QR分解角度看線性最小二乘

對於一個over-determined線性最小二乘問題 $Ax\simeq b$ ，其目標函數是 $\begin{align} \phi(x) = || r(x) ||_2^2 &= || b - Ax ||_2^2 = || b - Q\begin{bmatrix}\hat{R} \\ 0\end{bmatrix} x ||_2^2 \\ &= || Q^T (b - Q\begin{bmatrix}\hat{R} \\ 0\end{bmatrix} x) ||_2^2 \\ &= || Q^Tb - \begin{bmatrix}\hat{R} \\ 0\end{bmatrix} x ||_2^2 \end{align} \tag{2}$

這里 $Q\in \mathbb{R}^{m\times m}$ ， $Qb\in \mathbb{R}^{m\times 1}$ ， $\begin{bmatrix}\hat{R} \\ 0\end{bmatrix} \in \mathbb{R}^{m\times n}$ ， $\begin{bmatrix}\hat{R} \\ 0\end{bmatrix}x \in \mathbb{R}^{m\times 1}$ 。

如果把 $Q^Tb$ 拆分成上下兩部分，形式 $\begin{bmatrix}\hat{R} \\ 0\end{bmatrix}$ 類似， $Q^Tb = \begin{bmatrix} c_1\\c_2 \end{bmatrix},\ where\ c_1 \in \mathbb{R}^n,\ c_2\in\mathbb{R}^{m-n}$ 。那么目標函數可以寫成下面的形式：

$|| r(x) ||_2^2 = || c_1 - \hat{R}x ||_2^2 + || c_2||_2^2 \tag{3}$

可以看到，我們只能最小化前一部分 $|| c_1 - \hat{R}x ||_2^2$ 到0，即 $\hat{R}x = c_1$ ， $|| r(x) ||_2^2$ 的最小值為 $|| c_2||_2^2$ 。這樣處理之后就避免了求正規方程組中的 $(A^TA)^{-1}$ ，避免了條件數變成 $cond(A^TA) = cond(A)^2$ ，所以QR分解法更加數值穩定。

1.4 計算QR分解的方法

一共有三種：

Gram–Schmidt Orthogonalization
Householder Triangularization
Givens Rotations

1.5 Gram–Schmidt Orthogonalization

1.5.1 Reduced QR分解

GSO構建正交矩陣 $Q\in \mathbb{R}^{m\times m}$ 的方法是從A矩陣的n個列（ $A_{:,j} \in \mathbb{R}^{m\times 1}$ ）中構建互相正交的基，先選定 $A_{:,0}$ 為第一個基，然后把第二列 $A_{:,1}$ 減去平行於 $A_{:,0}$ 的部分，剩下的垂直於 $A_{:,0}$ 的部分作為下一個基，以此類推，直到生成了n個基。

$\begin{aligned} A_{:,0} &= r_{00}q_0 \\ A_{:,1} &= r_{0,1}q_0 + r_{1,1}q_1 \\ \vdots\\ A_{:,n-1} &= r_{0,n-1}q_0 + r_{1,n-1}q_1 + \dots + r_{n-1,n-1}q_{n-1}\\ A &= \hat{Q}\hat{R} \end{aligned} \tag{4}$

$\begin{bmatrix} A_{:,0} | A_{:,1} | \dots | A_{:,n-1} \end{bmatrix} = \begin{bmatrix} q_{0} | q_{1} | \dots | q_{n-1} \end{bmatrix} \begin{bmatrix} r_{0,0} & r_{0,1} & \dots & r_{0,n-1}\\ & r_{1,1} & \dots & r_{1,n-1}\\ & & \ddots & \vdots \\ & & & r_{n-1,n-1}\\ \end{bmatrix}\\ \tag{5}$

這個方法生成的 $\hat{Q}\in \mathbb{R}^{m\times n}$ ， $\hat{R} \in \mathbb{R}^{n\times n}$ ，和section1.1中定義的Q是方陣，R不是方陣有區別。這個結果被稱為Reduced QR分解，因為m>n，所以只滿足 $\hat{Q}^T\hat{Q} = I$ ，而不滿足 $\hat{Q}\hat{Q}^T = I$ 。

Credit to http://iacs-courses.seas.harvard.edu/courses/am205/schedule/

Reduced QR分解同樣可以求解over-determined線性最小二乘問題。形式類似Full QR分解：

$\hat{R}x = \hat{Q}^T b \tag{6}$

其中 $\hat{R}x \in \mathbb{R}^{n\times 1}$ ， $\hat{Q}^Tb \in \mathbb{R}^{n\times 1}$ 。

1.5.2 Full QR分解

為了實現定義中的完整的QR分解，需要把上面生成Q中的n個基拓展成m個互相正交的基。但此處並沒有對額外的m-n個基的順序有特殊要求，因此任意一種順序都可以。另外還需要把 $\hat{R}$ 下面加m-n行零矩陣。

在Python中，Reduced QR分解和Full QR分解對應於

q,r = np.linalg.qr(A) # reduced q,r = np.linalg.qr(A,mode="complete") # full

1.5.3 Classic Gram–Schmidt Orthogonalization算法 CGSO

觀察Eq.4可以發現，其實每一步迭代都只有一個 $r_{j,j}$ 未知：左邊 $A_{:,0}$ 已知，右邊 $q_0 \sim q_{j-1}$ 已知， q的系數們 $r_{0,j} \sim r_{j-1,j}$ 可以用公式 $r_{i,j} = q_i^T A_{:,j}$ 求得。把 $r_{i,j} = q_i^T A_{:,j}$ 代入Eq.4，並整理可得

$r_{j,j}q_{j} = v_j= A_{:,j} - (q_0^T A_{:,j})q_0 - (q_1^T A_{:,j})q_1 - \dots (q_{j-1}^T A_{:,j})q_{j-1} \tag{7}$

因此 $q_j = \frac{v_j}{||v_j||_2}$ ， $|r_{j,j}| = ||v_j||_2$ 。其中， $r_{j,j}$ 的符號不確定是因為，任意一個基方向反向之后，這個QR分解不會有任何問題，這個基仍然和其他基正交。為了計算方便，這里就規定 $r_{j,j} = ||v_j||_2$ 。

整理上面計算 $q_j$ 和 $r_{i,j}$ 的過程為算法的形式：

Credit to http://iacs-courses.seas.harvard.edu/courses/am205/schedule/

觀察算法過程，可以發現，唯一可能在理論上出問題的情況就是，出現某個 $r_{j,j}$ =0，導致在算法第8行出現0在分母上的情況。因此只要 $A\in\mathbb{R}^{m\times n}$ 是滿秩的，且每個 $r_{jj}$ 都>0，那么reduced QR分解的結果是唯一的。

1.5.4 Modified Gram–Schmidt Orthogonalization算法 MGSO

由於CGSO對舍入誤差很敏感，容易導致生成的基 $q_j$ 的正交性隨着迭代越來越弱，因此引入改進的GSO。核心思想是，在每個 $q_j$ 生成后，直接把A剩下的列（下面算法第7行）都去掉 $q_j$ 的成分（下面算法的第8-9行）。因為只是把計算的順序變了，所以理論上計算結果是一樣的。

Credit to http://iacs-courses.seas.harvard.edu/courses/am205/schedule/

但是改進之后穩定性會好很多。從實際計算步驟上來看，CGSO和MGSO的區別在於，CGSO中，每次迭代新的一列 $v_i$ ，計算每個 $r_{i,j} = q_i^Tv_j$ 都是用的同一個 $v_j$ ，而MGSO計算 $r_{i,j}$ 的時候用的 $v_j$ 是已經減去前面j-1個基的分量之后的 $v_j$ 。

這樣做的好處是：誤差的傳遞是局部的。比如計算 $q_1$ 是精確的，計算 $q_2$ 出現誤差，即， $q_2$ 在 $q_1$ 上存在一個微小分量。按照CGSO，接下來要分別計算 $v_3$ 在 $q_1$ 和 $q_2$ 的分量，最終 $q_1^Tv_3 \ne 0,\ q_2^Tv_3 \ne 0$ ；而MGSO則先計算 $v_3$ 在 $q_1$ 上的分量，去除掉這個分量之后成為 $v_3'$ ，再計算並去除 $v_3'$ 在 $q_2$ 上的分量得到最終的 $v_3''$ ，此時如果計算是精確的，那么至少可以保證 $v_3'' \perp q_2$ 。

直觀理解參考下面這張圖，在三維xyz坐標系里， $q_2'$ 是帶誤差的 $q_2$ 。用CGSO處理 $q_2'$ 的時候， $v_3$ 用的是初始值 $(1,1,1)$ ，包含了 $q_1$ 和 $q_2$ 兩個方向的誤差，而用MGSO處理 $q_2'$ 的時候， $v_3$ 用的是去掉 $q_1$ 分量之后的 $v_3' = (0,1,1)$ ，只有 $q_1$ 方向的誤差。

公式上計算這些誤差參考The modified Gram-Schmidt procedure：

Credit to https://www.math.uci.edu/~ttrogdon/105A/html/Lecture23.html

1.6 Givens Rotations

1.6.1 Givens Rotation Matrix

$G(i,j,\theta) = \begin{bmatrix} 1 & \dots & 0 & \dots & 0 & \dots & 0\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 0 & \dots & cos\theta & \dots & sin\theta & \dots & 0\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 0 & \dots & -sin\theta & \dots & cos\theta & \dots & 0\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots\\ 0 & \dots & 0 & \dots & 0 & \dots & 1\\ \end{bmatrix} \tag{8}$

$\begin{aligned} i &< j \\ G_{i,i} &= c\\ G_{j,j} &= c\\ G_{i,j} &= s\\ G_{j,i} &= -s\\ G_{k,k} &= 1,\ for\ k\ne i\ or\ j\\ G_{t,s} &= 0,\ otherwise \end{aligned} \tag{9}$

1.6.2 Givens Rotations的作用

對於一個矩陣 $A\in \mathbb{R}^{m\times n},\ m\ge n$ ，對於第i列的第j和k行 $A_{i,j},A_{i,k},\ j<k$ ，如果 $A_{i,k}$ 元素不為0，可以通過一個Givens Rotation把它轉換成0。

$\begin{bmatrix} c & s\\ -s & c \end{bmatrix} \begin{bmatrix} A_{i,j}\\ A_{i,k} \end{bmatrix} = \begin{bmatrix} \alpha\\ 0 \end{bmatrix}\\ \alpha = \sqrt{A_{i,j}^2 + A_{i,k}^2 }\\ c = \frac{A_{i,j}}{\alpha}\\ s = \frac{A_{i,k}}{\alpha}\\ \tag{10}$

當 $A_{i,j}$ 或者 $A_{i,k}$ 很小或者很大，且它們的平方不是用float表示的時候，對它們求平方會導致上溢出或者下溢出。因此更好的公式是：

如果 $|A_{i,j}| > | A_{i,k}|$ ，那么設 $t = A_{i,k}/A_{i,j}$

$c = \frac{1}{\sqrt{1+t^2}},\ s = ct$

如果 $|A_{i,j}| < | A_{i,k}|$ ，那么設 $\tau = A_{i,j}/A_{i,k}$

$s = \frac{1}{\sqrt{1+\tau^2}},\ c = s\tau$

不過這個問題基本只有在設計package造輪子的時候才會遇到，所以通常用Eq.10不會引起問題。詳見Scientific Computing - Heath的第128頁。

另外，在涉及反三角的數值運算的時候，建議使用atan2替代atan，范圍更大，更穩定。例如atan2(y,x)會返回一個(x,y)向量和正x軸的夾角。

the difference between atan and atan2 in C++?stackoverflow.com wikipedia Atan2en.wikipedia.org

1.6.3 Givens Rotations 算法

對於一個稠密的矩陣 $A\in \mathbb{R}^{m\times n},\ m\ge n$ ，逐漸把A消元成R（參考1.5.1的full QR的圖）。

Credit to http://iacs-courses.seas.harvard.edu/courses/am205/schedule/

注意第三行的循環，j是從大到小的迭代。

1.6.4 Givens Rotations 優勢

當A是稠密矩陣，Givens Rotations並沒有比另外兩種算法更高效，但如果A是稀疏矩陣，那么Givens Rotations大小為0的元素可以直接被忽略。另一個優勢是，Givens Rotations更容易並行化，因為Givens Rotations只對兩個元素進行操作，處理不同列的時候可以完全的獨立。

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 QR分解矩陣QR分解矩陣分解系列二：正交三角分解(UQ、QR分解) 矩陣QR分解的MATLAB與C++實現機器學習中的矩陣方法03：QR 分解數值計算（Python實現）（一） Matlab數值計算示例：牛頓插值法、LU分解法、拉格朗日插值法、牛頓插值法 QR分解迭代求特征值——原生python實現(不使用numpy) 數值分析筆記（3）——數值計算中的原則數值分析之奇異值分解(SVD)篇