[数值计算] 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分解 QR分解矩阵QR分解 QR分解与最小二乘【矩阵】RQ/QR 分解 QR分解与线性回归矩阵分解---QR正交分解，LU分解简单QR分解之Gram-Schmit正交化&&Householder变换&&Givens Rotation变换&&计算步骤 8.QR分解的python实现