本文详细介绍了线性代数中常见的矩阵分解方法：Cholesky分解、QR分解和LU分解。内容包括它们的定义、数值方法、稳定性分析以及算法复杂度。Cholesky分解用于正定矩阵，QR分解适用于任何实数矩阵，LU分解则关注矩阵的LU因子是否存在。文章还讨论了这些分解在数值稳定性上的挑战和应用，例如在凸优化问题和线性系统求解中的重要性。
本篇尝试介绍几个常见的矩阵分解及其数值方法，探究其算法复杂度及不稳定性来源，能够让我们更好的理解工具，使用工具。常见的矩阵分解可以在维基中找到
Cholesky
定义
正定的矩阵（如果不能分解的对称矩阵，一定是有非正的特征值），分解为对称的up down矩阵的乘积
即对于
A
∈
S
(
R
)
N
×
N
+
+
A \in \mathcal{S}(\mathbb{R})^{++}_{N\times N}
A∈S(R)N×N++
一定存在 L一个下三角矩阵（left bottom) 矩阵，能够使得有
A
=
L
L
T
A = L L^T
A=LLT
对于复数依旧成立，但是会要求A是一个Hermitian矩阵（即共轭转置等于自己）
数值方法
对于Cholesky求解有两种方法，参考维基的内容，第一种是利用高斯消元法，这种方法称为Chelosky 算法，第二种主要精神是首先假设这样的L存在，并反解，称为The Cholesky–Banachiewicz and Cholesky–Crout算法。
-
Chelosky算法
对于一个矩阵A，我们希望通过迭代的方式逐渐消掉A的元素，假设在第i步迭代我们已经有了如下矩阵（当然左上角的那个Id矩阵可以是size 0，这就使得这个矩阵和原矩阵A完全一样）。
A ( i ) = ( I i − 1 0 0 0 a i , i b i ∗ 0 b i B ( i ) ) , \mathbf{A}^{(i)}=\begin{pmatrix} \mathbf{I}_{i-1} & 0 & 0 \\ 0 & a_{i,i} & \mathbf{b}_{i}^{*} \\ 0 & \mathbf{b}_{i} & \mathbf{B}^{(i)} \end{pmatrix}, A(i)=⎝⎛Ii−1000ai,ibi0bi∗B(i)⎠⎞,
因为A是一个 S ( R ) N × N + + \mathcal{S}(\mathbb{R})^{++}_{N\times N} S(R)N×N++, 因此 a i , i a_{i,i} ai,i一定是一个严格大于零的数字（如果非严格大于零，我们可以找到一个基底函数 x = ( 0 , . . , 1 , . . 0 ) T x=(0,..,1,..0)^T x=(0,..,1,..0)T，只在i行有数值，这样 x T A x = a i , i x^T A x = a_{i,i} xTAx=ai,i，根据A的正定定义，这个数值必须大于零）。
如果我们定义一个矩阵 L i L_i Li有如下形式
L i : = ( I i − 1 0 0 0 a i , i 0 0 1 a i , i b i I n − i ) \mathbf{L}_{i}:= \begin{pmatrix} \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{pmatrix} Li:=⎝⎛Ii−1000ai,iai,i1bi00In−i⎠⎞
那么就会有如下关系
A ( i ) = L i A ( i + 1 ) L i ∗ \mathbf{A}^{(i)} = \mathbf{L}_{i} \mathbf{A}^{(i+1)} \mathbf{L}_{i}^{*} A(i)=LiA(i+1)Li∗
其中
A ( i + 1 ) = ( I i − 1 0 0 0 1 0 0 0 B ( i ) − 1 a i , i b i b i ∗ ) \mathbf{A}^{(i+1)}= \begin{pmatrix} \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}^{*} \end{pmatrix} A(i+1)=⎝⎛Ii−10001000B(i)−ai,i1bibi∗⎠⎞
因此我们又有了下一步需要优化的部分，即 B ( i ) − 1 a i , i b i b i ∗ \mathbf{B}^{(i)} - \frac{1}{a_{i,i}} \mathbf{b}_{i} \mathbf{b}_{i}^{*} B(i)−ai,i1bibi∗。（存疑，为什么新的形式下这个矩阵的最左上角元素依然是大于零的？）。最后只需要让 L L L有如下形式即可。
L : = L 1 L 2 … L n \mathbf{L} := \mathbf{L}_{1} \mathbf{L}_{2} \dots \mathbf{L}_{n} L:=L1L2…Ln -
Cholesky–Banachiewicz and Cholesky–Crout算法
假设存在一个矩阵L
A
=
L
L
T
=
(
L
11
0
0
L
21
L
22
0
L
31
L
32
L
33
)
(
L
11
L
21
L
31
0
L
22
L
32
0
0
L
33
)
=
(
L
11
2
(
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
)
\begin{aligned} \mathbf{A} = \mathbf{LL}^T & = \begin{pmatrix} L_{11} & 0 & 0 \\ L_{21} & L_{22} & 0 \\ L_{31} & L_{32} & L_{33}\\ \end{pmatrix} \begin{pmatrix} L_{11} & L_{21} & L_{31} \\ 0 & L_{22} & L_{32} \\ 0 & 0 & L_{33} \end{pmatrix} \\[8pt] & = \begin{pmatrix} 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{pmatrix} \end{aligned}
A=LLT=⎝⎛L11L21L310L22L3200L33⎠⎞⎝⎛L1100L21L220L31L32L33⎠⎞=⎝⎛L112L21L11L31L11L212+L222L31L21+L32L22(symmetric)L312+L322+L332⎠⎞
于是便有
L
=
(
A
11
0
0
A
21
/
L
11
A
22
−
L
21
2
0
A
31
/
L
11
(
A
32
−
L
31
L
21
)
/
L
22
A
33
−
L
31
2
−
L
32
2
)
\begin{aligned} \mathbf{L} = \begin{pmatrix} \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{pmatrix} \end{aligned}
L=⎝⎛A11A21/L11A31/L110A22−L212(A32−L31L21)/L2200A33−L312−L322⎠⎞
只要逐行（或者逐渐列）处理，就可以算出每一个L中的元素
L
j
,
j
=
(
±
)
A
j
,
j
−
∑
k
=
1
j
−
1
L
j
,
k
2
L
i
,
j
=
1
L
j
,
j
(
A
i
,
j
−
∑
k
=
1
j
−
1
L
i
,
k
L
j
,
k
)
for
i
>
j
\begin{aligned} 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) \quad \text{for } i>j \end{aligned}
Lj,jLi,j=(±)Aj,j−k=1∑j−1Lj,k2=Lj,j1(Ai,j−k=1∑j−1Li,kLj,k)for i>j
数值方法的数值稳定性
这个数值方法的不稳定性来源：（误差来源？）如果对角线的数量级相差很大，那么有可能造成四舍五入的时候累计误差会（累计在之前计算Ljk的时候）让最后的Ljj变成了对于一个负数求解。如果矩阵非常的不规则(ill conditioned, 一定程度可以用最大最小的特征值比值作为参考)，有可能会造成无解或者实数矩阵变成复空间的解。
数值方法的复杂度
两个算法的复杂度大约为为 1 / 3 n 3 1/3 n^3 1/3n3，第二种相对会更好一点，因为对于元素的iteration相对更加规律，对于内存使用率会更好一些。这个算法相比于LU要快一倍，LU目前最好的approximation大概是 2 / 3 n 3 2/3 n^3 2/3n3。 同样，QR分解的最好办法是Gram-Schmidt，复杂度也是 2 / 3 n 3 2/3 n^3 2/3n3。
应用
对于凸优化问题非常有用，而且相对最快。如果有一个covariance matrix，就能够很好的将其转换成二次扁平化的一个乘积。
x
T
A
x
=
∣
∣
L
T
x
∣
∣
2
x^TAx = ||L^Tx||^2
xTAx=∣∣LTx∣∣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。
其他性质：
- 如果A可逆，那么可以找到唯一的R，使得R的的对角线全是正数。
- Q的前k列是一个orthonormal的base，正好能够span生成A前k列的空间
- 因为R是一个upper triangle matrix，A的前K列只与Q的前K列有线性关系
数值方法
- Gram-Schmidt
使用迭代的方法，一列一列的去掉当前列与前面的垂直部分。因为在第k步的时候已经求出了能够span前k-1列的垂直归一向量，那么在第k步的时候去掉第k列在之前k-1个基向量上面构成超平面的投影向量，再去做归一处理即可。具体可以参考维基
这个方法已经几乎被弃用，虽然逻辑清晰，但是数值不稳定。可以看到下一步的垂直向量是收到之前所有误差的累加，误差随矩阵规模（应该是线性）增加，因此数量级差的大的话会有问题；而下面这种方法，误差并不会累加，可能会加加减减，保持同一数量级或者累乘（累乘的话误差越来越小，数值稳定） - Householder reflections
巧妙构造一个orthonormal的矩阵，且这个矩阵能够QA之后去掉第一列除了第一个元素其他的数值，逐步实现一列一列消除。
对于任意的向量x，使得 α \alpha α = x的模，于是就有这样的构造使得Q是orthonormal (可以验证 Q Q T = I d ) QQ^T=Id) QQT=Id)。
u = x − α e 1 , v = u ∥ u ∥ , Q = I − 2 v v T . \begin{aligned} \mathbf{u} &= \mathbf{x} - \alpha\mathbf{e}_1, \\ \mathbf{v} &= \frac{\mathbf{u}}{\|\mathbf{u}\|}, \\ Q &= I - 2 \mathbf{v}\mathbf{v}^\textsf{T}. \end{aligned} uvQ=x−αe1,=∥u∥u,=I−2vvT.
如果A是一个复空间矩阵，那么有 Q = I − 2 v v ∗ Q = I - 2\mathbf{v}\mathbf{v}^* Q=I−2vv∗。可以验证，
Q x = [ α 0 ⋮ 0 ] Q\mathbf{x} = \begin{bmatrix} \alpha \\ 0 \\ \vdots \\ 0 \end{bmatrix} Qx=⎣⎢⎢⎢⎡α0⋮0⎦⎥⎥⎥⎤
如果我们让之前的x取矩阵A的第一列，记 Q 1 = Q Q_1 = Q Q1=Q，便有了：
Q 1 A = [ α 1 ⋆ ⋯ ⋆ 0 ⋮ A ′ 0 ] Q_1A = \begin{bmatrix} \alpha_1 & \star & \cdots & \star \\ 0 & & & \\ \vdots & & A' & \\ 0 & & & \end{bmatrix} Q1A=⎣⎢⎢⎢⎡α10⋮0⋆⋯A′⋆⎦⎥⎥⎥⎤
对A’进行同样的操作，便会消掉A的第一列（注意这里的A’尺寸更小了），只要在左上方补上Id就可以保证Q的那部分效果不变。于是在消元第k列的时候，只需要让新的 Q k Q_k Qk有如下形式便能够让 A k ′ A_k' Ak′消掉第一列除了最上角元素的其他数值
Q k = [ I k − 1 0 0 Q k ′ ] Q_k = \begin{bmatrix} I_{k-1} & 0 \\ 0 & Q_k' \end{bmatrix} Qk=[Ik−100Qk′]
如果把所求到的所有 Q i Q_i Qi累乘起来，就能够得到一个右上角矩阵R
R = Q t ⋯ Q 2 Q 1 A R = Q_t \cdots Q_2 Q_1 A R=Qt⋯Q2Q1A
又因为这些Q都是orthonormal，因此只需要两边分别乘以这些Q的转置就有了最终的结果
Q = Q 1 T Q 2 T ⋯ Q t T = Q 1 Q 2 ⋯ Q t A = Q R \begin{aligned} Q &= Q_1^\textsf{T} Q_2^\textsf{T} \cdots Q_t^\textsf{T} = Q_1 Q_2 \cdots Q_t \\ A & =QR \end{aligned} QA=Q1TQ2T⋯QtT=Q1Q2⋯Qt=QR
算法复杂度约为 2 / 3 n 3 2/3 n^3 2/3n3。算法稳定性也比较好， 这个算法的问题是无法并行，一次性需要的内存较多。 - Given Rotation
利用四角rotation矩阵
G 1 = [ cos ( θ ) 0 − sin ( θ ) 0 1 0 sin ( θ ) 0 cos ( θ ) ] G_1 = \begin{bmatrix} \cos(\theta) & 0 & -\sin(\theta) \\ 0 & 1 & 0 \\ \sin(\theta) & 0 & \cos(\theta) \end{bmatrix} G1=⎣⎡cos(θ)0sin(θ)010−sin(θ)0cos(θ)⎦⎤
逐步去掉左下角的元素，最后将各种G累乘起来，这个算法求解路径比较复杂，而且求三角函数的反函数也会有比较大的误差和计算量，虽然巧妙但是不实用。
应用
- 求解determinant，因为Q的det是1，因此只需要把R的对角乘积求出来就可以了
- 线性问题求解，这种方法比直接求逆来的更快速且数值更稳定
LU
定义
所谓LU就是存在一个两个triangle矩阵，对于一个矩阵A，存在只有对角线下面/上面数值非零，有形如：
A
=
L
U
A = LU
A=LU
但是这个分解的存在性要取决于A矩阵的arrangement，如果前面存在leading principal minor为零的情况，可能会导致无解。当然这个分解并不一定唯一。
条件
- 任何NxN方阵都有LUP和PLU分解
- 当A可逆时，有LU分解 等价于 任意 leading principal minors（即去掉后面的n-k行/列后剩下的方阵的determinant）非零。
- 如果A不可逆但是前k 个leading principal minor非零，那么他也有LU分解
因此对于一个方阵存在三种可能
- 唯一的LU分解（当可逆矩阵存在LDU分解， A = L D U A = LDU A=LDU 此时L和U都是unitriangular的矩阵，D是对角线矩阵，即对角线全是1，此时找到的L或U就是LDU分解中的（LD）或者（DU））
- 存在无穷多个LU分解形式（当前n-1列中存在两个或以上的列是线性相关，或者前n-1列中存在0向量）
- 不存在LU分解（当前n-1列中无非零向量，且线性不想管，并且存在一个leading principal minor值为零时）（对于这种情况，可以对A矩阵中某一个对角线上的值加一个误差 ϵ \epsilon ϵ来近似LU分解的值
如果这个矩阵恰好是 H e r m i t i a n Hermitian Hermitian，那么LU分解就成了Cholesky分解
数值方法
- 高斯消元法，复杂度大概在 2 / 3 n 3 2/3 n^3 2/3n3（由于消元过程中误差会累计，如果对角线数值比较极端，那么有可能造成数值不稳定）
- recursion
- randomized
- 据维基描述，这个算法的复杂度大约和矩阵相乘差不多 O ( n 2.376 ) O(n^{2.376}) O(n2.376)
应用
- 求解线性系统， 矩阵求逆（为什么不用QR？）
- 计算Determinant （为什么不用QR？）
扫一扫
8546
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
