[数值计算] QR分解

[数值计算] QR分解
https://zhuanlan.zhihu.com/p/84415000

0. 为什么要用QR分解

$Ax=b, Ainmathbb{R}^{m imes 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 imes n}, mge n$ 可以被分解成 $A = QR$ ，其中：
- $Qin mathbb{R}^{m imes m}$ 是正交矩阵
- $R equiv egin{bmatrix} hat{R} \ 0 end{bmatrix} in mathbb{R}^{m imes n}$
- $hat{R} in mathbb{R}^{n imes n}$ 是上三角矩阵
1.2 正交矩阵的性质
- $Q^T Q = QQ^T=I$
- 左乘一个正交矩阵对欧式范数的结果不影响（在下面证明eq.2的时候会用到）
$|| Qv ||_2^2 = v^TQ^TQv = v^Tv = || v ||_2^2 ag{1}$

1.3 从QR分解角度看线性最小二乘

对于一个over-determined线性最小二乘问题 $Axsimeq b$ ，其目标函数是 $egin{align} phi(x) = || r(x) ||_2^2 &= || b - Ax ||_2^2 = || b - Qegin{bmatrix}hat{R} \ 0end{bmatrix} x ||_2^2 \ &= || Q^T (b - Qegin{bmatrix}hat{R} \ 0end{bmatrix} x) ||_2^2 \ &= || Q^Tb - egin{bmatrix}hat{R} \ 0end{bmatrix} x ||_2^2 end{align} ag{2}$

这里 $Qin mathbb{R}^{m imes m}$ ， $Qbin mathbb{R}^{m imes 1}$ ， $egin{bmatrix}hat{R} \ 0end{bmatrix} in mathbb{R}^{m imes n}$ ， $egin{bmatrix}hat{R} \ 0end{bmatrix}x in mathbb{R}^{m imes 1}$ 。

如果把 $Q^Tb$ 拆分成上下两部分，形式 $egin{bmatrix}hat{R} \ 0end{bmatrix}$ 类似， $Q^Tb = egin{bmatrix} c_1\c_2 end{bmatrix}, where c_1 in mathbb{R}^n, c_2inmathbb{R}^{m-n}$ 。那么目标函数可以写成下面的形式：

$|| r(x) ||_2^2 = || c_1 - hat{R}x ||_2^2 + || c_2||_2^2 ag{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构建正交矩阵 $Qin mathbb{R}^{m imes m}$ 的方法是从A矩阵的n个列（ $A_{:,j} in mathbb{R}^{m imes 1}$ ）中构建互相正交的基，先选定 $A_{:,0}$ 为第一个基，然后把第二列 $A_{:,1}$ 减去平行于 $A_{:,0}$ 的部分，剩下的垂直于 $A_{:,0}$ 的部分作为下一个基，以此类推，直到生成了n个基。

$egin{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} ag{4}$

$egin{bmatrix} A_{:,0} | A_{:,1} | dots | A_{:,n-1} end{bmatrix} = egin{bmatrix} q_{0} | q_{1} | dots | q_{n-1} end{bmatrix} egin{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}\ ag{5}$

这个方法生成的 $hat{Q}in mathbb{R}^{m imes n}$ ， $hat{R} in mathbb{R}^{n imes n}$ ，和section1.1中定义的Q是方阵，R不是方阵有区别。这个结果被称为Reduced QR分解，因为m>n，所以只满足 $hat{Q}^That{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 ag{6}$

其中 $hat{R}x in mathbb{R}^{n imes 1}$ ， $hat{Q}^Tb in mathbb{R}^{n imes 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} ag{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在分母上的情况。因此只要 $Ainmathbb{R}^{m imes 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 e 0, q_2^Tv_3 e 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, heta) = egin{bmatrix} 1 & dots & 0 & dots & 0 & dots & 0\ vdots & ddots & vdots & ddots & vdots & ddots & vdots\ 0 & dots & cos heta & dots & sin heta & dots & 0\ vdots & ddots & vdots & ddots & vdots & ddots & vdots\ 0 & dots & -sin heta & dots & cos heta & dots & 0\ vdots & ddots & vdots & ddots & vdots & ddots & vdots\ 0 & dots & 0 & dots & 0 & dots & 1\ end{bmatrix} ag{8}$

$egin{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 e i or j\ G_{t,s} &= 0, otherwise end{aligned} ag{9}$

1.6.2 Givens Rotations的作用

对于一个矩阵 $Ain mathbb{R}^{m imes n}, mge n$ ，对于第i列的第j和k行 $A_{i,j},A_{i,k}, j<k$ ，如果 $A_{i,k}$ 元素不为0，可以通过一个Givens Rotation把它转换成0。

$egin{bmatrix} c & s\ -s & c end{bmatrix} egin{bmatrix} A_{i,j}\ A_{i,k} end{bmatrix} = egin{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}\ ag{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}|$ ，那么设 $au = A_{i,j}/A_{i,k}$
$s = frac{1}{sqrt{1+ au^2}}, c = s au$

不过这个问题基本只有在设计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 算法

对于一个稠密的矩阵 $Ain mathbb{R}^{m imes n}, mge 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只对两个元素进行操作，处理不同列的时候可以完全的独立。
相关阅读:
win8.1下解决Visual C++不兼容的方法
 Java文件File类学习总结
 Java可视化日历（Date类、DATe Format类、Calendar类综合运用），开发可视化日历小程序
 java时间日期类（Date、DateFormat、Calendar）学习
 使用SimpleDateFormat类来实现时间跟字符串的转化
 Java中运行时异常和非运行时异常什么鬼？
Java构造器练习题
 Word中一条删除不掉的单或双横线的解决办法
 Java中常用类（包装类扩展知识）
Java异常学习总结二
原文地址：https://www.cnblogs.com/dhcn/p/14840472.html

[数值计算] QR分解

0. 为什么要用QR分解

1. QR分解

1.1 定义

1.2 正交矩阵的性质

1.3 从QR分解角度看线性最小二乘

1.4 计算QR分解的方法

1.5 Gram–Schmidt Orthogonalization

1.5.1 Reduced QR分解

1.5.2 Full QR分解

1.5.3 Classic Gram–Schmidt Orthogonalization算法 CGSO

1.5.4 Modified Gram–Schmidt Orthogonalization算法 MGSO

1.6 Givens Rotations

1.6.1 Givens Rotation Matrix

1.6.2 Givens Rotations的作用

1.6.3 Givens Rotations 算法

1.6.4 Givens Rotations 优势