[Math]常见矩阵分解及复杂度 Cholesky QR LU

本篇尝试介绍几个常见的矩阵分解及其数值方法，探究其算法复杂度及不稳定性来源，能够让我们更好的理解工具，使用工具。常见的矩阵分解可以在维基中找到

Cholesky

定义

正定的矩阵（如果不能分解的对称矩阵，一定是有非正的特征值），分解为对称的up down矩阵的乘积
即对于
$$A\in \mathcal{S}(\mathbb{R}{)}_{N\times N}^{++}$$
一定存在 L一个下三角矩阵（left bottom) 矩阵，能够使得有
$$A=L{L}^T$$
对于复数依旧成立，但是会要求A是一个Hermitian矩阵（即共轭转置等于自己）

数值方法

对于Cholesky求解有两种方法，参考维基的内容，第一种是利用高斯消元法，这种方法称为Chelosky 算法，第二种主要精神是首先假设这样的L存在，并反解，称为The Cholesky–Banachiewicz and Cholesky–Crout算法。

1. Chelosky算法
   对于一个矩阵A，我们希望通过迭代的方式逐渐消掉A的元素，假设在第i步迭代我们已经有了如下矩阵（当然左上角的那个Id矩阵可以是size 0，这就使得这个矩阵和原矩阵A完全一样）。
   $${\mathbf{A}}^{(i)}=\left(\begin{array}{ccc}{\mathbf{I}}_{i-1}& 0& 0\\ 0& a_{i,i}& {\mathbf{b}}_i^{\ast }\\ 0& {\mathbf{b}}_i& {\mathbf{B}}^{(i)}\end{array}\right),$$
   因为A是一个 $\mathcal{S}(\mathbb{R}{)}_{N\times N}^{++}$, 因此$a_{i,i}$一定是一个严格大于零的数字（如果非严格大于零，我们可以找到一个基底函数$x=(0,..,1,..0{)}^T$，只在i行有数值，这样$x^{T}Ax=a_{i,i}$，根据A的正定定义，这个数值必须大于零）。
   如果我们定义一个矩阵$L_i$有如下形式
   $${\mathbf{L}}_i:=\left(\begin{array}{ccc}{\mathbf{I}}_{i-1}& 0& 0\\ 0& \sqrt{a_{i,i}}& 0\\ 0& \frac{1}{\sqrt{a_{i,i}}}{\mathbf{b}}_i& {\mathbf{I}}_{n-i}\end{array}\right)$$
   那么就会有如下关系
   $${\mathbf{A}}^{(i)}={\mathbf{L}}_i{\mathbf{A}}^{(i+1)}{\mathbf{L}}_i^{\ast }$$
   其中
   $${\mathbf{A}}^{(i+1)}=\left(\begin{array}{ccc}{\mathbf{I}}_{i-1}& 0& 0\\ 0& 1& 0\\ 0& 0& {\mathbf{B}}^{(i)}-\frac{1}{a_{i,i}}{\mathbf{b}}_i{\mathbf{b}}_i^{\ast }\end{array}\right)$$
   因此我们又有了下一步需要优化的部分，即${\mathbf{B}}^{(i)}-\frac{1}{a_{i,i}}{\mathbf{b}}_i{\mathbf{b}}_i^{\ast }$。（存疑，为什么新的形式下这个矩阵的最左上角元素依然是大于零的？）。最后只需要让$L$有如下形式即可。
   $$\mathbf{L}:={\mathbf{L}}_1{\mathbf{L}}_2\dots {\mathbf{L}}_n$$

2. Cholesky–Banachiewicz and Cholesky–Crout算法
   假设存在一个矩阵L

$$\begin{array}{rl}\mathbf{A}={\mathbf{L}\mathbf{L}}^T& =\left(\begin{array}{ccc}{L}_{11}& 0& 0\\ L_{21}& L_{22}& 0\\ L_{31}& L_{32}& L_{33}\end{array}\right)\left(\begin{array}{ccc}{L}_{11}& L_{21}& L_{31}\\ 0& L_{22}& L_{32}\\ 0& 0& L_{33}\end{array}\right)\\ & =\left(\begin{array}{ccc}{L}_{11}^2& & (\text{symmetric})\\ L_{21}{L}_{11}& L_{21}^2+L_{22}^2& \\ L_{31}{L}_{11}& L_{31}{L}_{21}+L_{32}{L}_{22}& L_{31}^2+L_{32}^2+L_{33}^2\end{array}\right)\end{array}$$
于是便有
$$\begin{array}{r}\mathbf{L}=\left(\begin{array}{ccc}\sqrt{A_{11}}& 0& 0\\ A_{21}/L_{11}& \sqrt{A_{22}-L_{21}^2}& 0\\ A_{31}/L_{11}& \left(A_{32}-L_{31}{L}_{21}\right)/L_{22}& \sqrt{A_{33}-L_{31}^2-L_{32}^2}\end{array}\right)\end{array}$$
只要逐行（或者逐渐列）处理，就可以算出每一个L中的元素
$$\begin{array}{rl}{L}_{j,j}& =(\pm )\sqrt{A_{j,j}-\sum _{k=1}^{j-1}{L}_{j,k}^2}\\ L_{i,j}& =\frac{1}{L_{j,j}}\left(A_{i,j}-\sum _{k=1}^{j-1}{L}_{i,k}{L}_{j,k}\right)\text{for }i>j\end{array}$$

数值方法的数值稳定性

这个数值方法的不稳定性来源：（误差来源？）如果对角线的数量级相差很大，那么有可能造成四舍五入的时候累计误差会（累计在之前计算Ljk的时候）让最后的Ljj变成了对于一个负数求解。如果矩阵非常的不规则(ill conditioned, 一定程度可以用最大最小的特征值比值作为参考)，有可能会造成无解或者实数矩阵变成复空间的解。

数值方法的复杂度

两个算法的复杂度大约为为$1/3n^3$，第二种相对会更好一点，因为对于元素的iteration相对更加规律，对于内存使用率会更好一些。这个算法相比于LU要快一倍，LU目前最好的approximation大概是$2/3n^3$。 同样，QR分解的最好办法是Gram-Schmidt，复杂度也是$2/3n^3$。

应用

对于凸优化问题非常有用，而且相对最快。如果有一个covariance matrix，就能够很好的将其转换成二次扁平化的一个乘积。
$$x^{T}Ax=||L^{T}x|{|}^2$$

QR Decomposition

定义

任何一个实数方阵A都可以被分解为QR（复数也可以，Q就是一个unitary matrix)，Q是orthogonal，R是一个upper triangle； 如果A是个MxN，的矩形矩阵，且要求M>N，那么存在Q属于MxM方阵，A=QR，其中R只只有前M行有数值，且也是一个MxM的upper triangle。

其他性质：

1. 如果A可逆，那么可以找到唯一的R，使得R的的对角线全是正数。
2. Q的前k列是一个orthonormal的base，正好能够span生成A前k列的空间
3. 因为R是一个upper triangle matrix，A的前K列只与Q的前K列有线性关系

数值方法

1. Gram-Schmidt
   使用迭代的方法，一列一列的去掉当前列与前面的垂直部分。因为在第k步的时候已经求出了能够span前k-1列的垂直归一向量，那么在第k步的时候去掉第k列在之前k-1个基向量上面构成超平面的投影向量，再去做归一处理即可。具体可以参考维基
   这个方法已经几乎被弃用，虽然逻辑清晰，但是数值不稳定。可以看到下一步的垂直向量是收到之前所有误差的累加，误差随矩阵规模（应该是线性）增加，因此数量级差的大的话会有问题；而下面这种方法，误差并不会累加，可能会加加减减，保持同一数量级或者累乘（累乘的话误差越来越小，数值稳定）
2. Householder reflections
   巧妙构造一个orthonormal的矩阵，且这个矩阵能够QA之后去掉第一列除了第一个元素其他的数值，逐步实现一列一列消除。
   对于任意的向量x，使得$\alpha$ = x的模，于是就有这样的构造使得Q是orthonormal (可以验证$Q{Q}^T=Id)$。
   $$\begin{array}{rl}\mathbf{u}& =\mathbf{x}-\alpha {\mathbf{e}}_1,\\ \mathbf{v}& =\frac{\mathbf{u}}{\Vert \mathbf{u}\Vert },\\ Q& =I-2\mathbf{v}{\mathbf{v}}^{\text{T}}.\end{array}$$
   如果A是一个复空间矩阵，那么有$Q=I-2\mathbf{v}{\mathbf{v}}^{\ast }$。可以验证，
   $$Q\mathbf{x}=\left[\begin{array}{c}\alpha \\ 0\\ ⋮\\ 0\end{array}\right]$$
   如果我们让之前的x取矩阵A的第一列，记$Q_1=Q$，便有了：
   $$Q_{1}A=\left[\begin{array}{cccc}{\alpha }_1& \star & \cdots & \star \\ 0& & & \\ ⋮& & A^{\prime }& \\ 0& & & \end{array}\right]$$
   对A’进行同样的操作，便会消掉A的第一列（注意这里的A’尺寸更小了），只要在左上方补上Id就可以保证Q的那部分效果不变。于是在消元第k列的时候，只需要让新的$Q_k$有如下形式便能够让$A_k^{\prime }$消掉第一列除了最上角元素的其他数值
   $$Q_k=\left[\begin{array}{cc}{I}_{k-1}& 0\\ 0& Q_k^{\prime }\end{array}\right]$$
   如果把所求到的所有$Q_i$累乘起来，就能够得到一个右上角矩阵R
   $$R=Q_t\cdots Q_2{Q}_{1}A$$
   又因为这些Q都是orthonormal，因此只需要两边分别乘以这些Q的转置就有了最终的结果
   $$\begin{array}{rl}Q& =Q_1^{\text{T}}{Q}_2^{\text{T}}\cdots Q_t^{\text{T}}=Q_1{Q}_2\cdots Q_t\\ A& =QR\end{array}$$
   算法复杂度约为$2/3n^3$。算法稳定性也比较好， 这个算法的问题是无法并行，一次性需要的内存较多。
3. Given Rotation
   利用四角rotation矩阵
   $$G_1=\left[\begin{array}{ccc}\cos (\theta )& 0& -\sin (\theta )\\ 0& 1& 0\\ \sin (\theta )& 0& \cos (\theta )\end{array}\right]$$
   逐步去掉左下角的元素，最后将各种G累乘起来，这个算法求解路径比较复杂，而且求三角函数的反函数也会有比较大的误差和计算量，虽然巧妙但是不实用。

应用

• 求解determinant，因为Q的det是1，因此只需要把R的对角乘积求出来就可以了
• 线性问题求解，这种方法比直接求逆来的更快速且数值更稳定

LU

定义

所谓LU就是存在一个两个triangle矩阵，对于一个矩阵A，存在只有对角线下面/上面数值非零，有形如：
$$A=LU$$
但是这个分解的存在性要取决于A矩阵的arrangement，如果前面存在leading principal minor为零的情况，可能会导致无解。当然这个分解并不一定唯一。

条件

• 任何NxN方阵都有LUP和PLU分解
• 当A可逆时，有LU分解 等价于 任意 leading principal minors（即去掉后面的n-k行/列后剩下的方阵的determinant）非零。
• 如果A不可逆但是前k 个leading principal minor非零，那么他也有LU分解

因此对于一个方阵存在三种可能

1. 唯一的LU分解（当可逆矩阵存在LDU分解，$A=LDU$ 此时L和U都是unitriangular的矩阵，D是对角线矩阵，即对角线全是1，此时找到的L或U就是LDU分解中的（LD）或者（DU））
2. 存在无穷多个LU分解形式（当前n-1列中存在两个或以上的列是线性相关，或者前n-1列中存在0向量）
3. 不存在LU分解（当前n-1列中无非零向量，且线性不想管，并且存在一个leading principal minor值为零时）（对于这种情况，可以对A矩阵中某一个对角线上的值加一个误差$ϵ$来近似LU分解的值

如果这个矩阵恰好是$Hermitian$，那么LU分解就成了Cholesky分解

数值方法

• 高斯消元法，复杂度大概在$2/3n^3$（由于消元过程中误差会累计，如果对角线数值比较极端，那么有可能造成数值不稳定）
• recursion
• randomized
• 据维基描述，这个算法的复杂度大约和矩阵相乘差不多$O(n^{2.376})$

应用

• 求解线性系统， 矩阵求逆（为什么不用QR？）
• 计算Determinant （为什么不用QR？）