三十分钟理解:矩阵Cholesky分解,及其在求解线性方程组、矩阵逆的应用
写一篇关于Cholesky分解的文章,作为学习笔记,尽量一文看懂矩阵Cholesky分解,以及用Cholesky分解来求解对称正定线性方程组,以及求对称正定矩阵的逆的应用。
先简单理解下正定矩阵和半正定矩阵的定义[1][2][3]:
- 给定一个 n × n n\times n n×n 的实对称矩阵 A A A ,若对于任意长度为 n 的非零向量 X,有 X T A X > 0 X^TAX \gt 0 XTAX>0 恒成立,则矩阵 A A A 是一个正定矩阵。
- 给定一个
n
×
n
n\times n
n×n 的实对称矩阵
A
A
A ,若对于任意长度为 n 的非零向量 X,有
X
T
A
X
≥
0
X^TAX \geq 0
XTAX≥0 恒成立,则矩阵
A
A
A 是一个半正定矩阵。
直接Cholesky分解
定理——若 A ∈ R n × n A \in R^{n\times n} A∈Rn×n是对称正定矩阵,则存在一个对角元全为正数的下三角矩阵 L ∈ R n × n L \in R^{n\times n} L∈Rn×n,使得 A = L L T A=LL^T A=LLT成立。
L
L
L是一个下三角形,形式是这样的:
推导 A = L L T A=LL^T A=LLT:我们先令
A = [ a 11 A 21 T A 21 A 22 ] , L = [ l 11 0 L 21 L 22 ] , L T = [ l 11 L 21 T 0 L 22 T ] A = \left[ \begin{matrix} a_{11}&A_{21}^{T}\\ A_{21}&A_{22}\\ \end{matrix} \right], \quad L = \left[ \begin{matrix} l_{11}&0\\ L_{21}&L_{22}\\ \end{matrix} \right], \quad L^{T} = \left[ \begin{matrix} l_{11}&L_{21}^{T}\\ 0&L_{22}^{T}\\ \end{matrix} \right] A=[a11A21A21TA22],L=[l11L210L22],LT=[l110L21TL22T]
其中
a
11
a_{11}
a11和
l
11
l_{11}
l11是一个标量,
A
21
A_{21}
A21和
L
21
L_{21}
L21是一个列向量,
A
22
A_{22}
A22是一个n-1阶的方阵,而
L
22
L_{22}
L22是一个n-1阶的下三角形。那么:
[
a
11
A
21
T
A
21
A
22
]
=
[
l
11
0
L
21
L
22
]
[
l
11
L
21
T
0
L
22
T
]
=
[
l
11
2
l
11
L
21
T
l
11
L
21
L
21
L
21
T
+
L
22
L
22
T
]
\left[ \begin{matrix} a_{11}&A_{21}^{T}\\ A_{21}&A_{22}\\ \end{matrix} \right] = \left[ \begin{matrix} l_{11}&0\\ L_{21}&L_{22}\\ \end{matrix} \right] \left[ \begin{matrix} l_{11}&L_{21}^{T}\\ 0&L_{22}^{T}\\ \end{matrix} \right]= \left[ \begin{matrix} l_{11}^{2}&l_{11}L_{21}^{T}\\ l_{11}L_{21}&L_{21}L_{21}^{T}+L_{22}L_{22}^{T}\\ \end{matrix} \right]
[a11A21A21TA22]=[l11L210L22][l110L21TL22T]=[l112l11L21l11L21TL21L21T+L22L22T]
未知量只有标量
l
11
l_{11}
l11,列向量
L
21
L_{21}
L21,和下三角形
L
22
L_{22}
L22,也是我们要求的。很容易得到:
l
11
=
a
11
l_{11} = \sqrt {a_{11}}
l11=a11
L 21 = 1 l 11 A 21 L_{21} = \frac {1}{l_{11}}A_{21} L21=l111A21
L 22 L 22 T = A 22 − L 21 L 21 T L_{22}L_{22}^{T} = A_{22} - L_{21}L_{21}^{T} L22L22T=A22−L21L21T
其中 l 11 l_{11} l11, L 21 L_{21} L21我们直接可以求出来了,并且可以求出 A 22 ′ = A 22 − L 21 L 21 T A_{22}' = A_{22} - L_{21}L_{21}^{T} A22′=A22−L21L21T。
而 A 22 ′ = L 22 L 22 T A_{22}' = L_{22}L_{22}^{T} A22′=L22L22T又是一个Cholesky分解!被分解的矩阵是一个n-1阶方阵 A 22 ′ A_{22}' A22′。因此,Cholesky分解算法具有递归性质,每一轮可以求出 L L L的一列,依次往下求,就可以把整个L求出来。
一个例子[4]:
另外,上述的方法需要进行开方,这有可能损失精度和增加运算量,为了避免开方,Cholesky分解有个改进的版本。将对称正定矩阵通过分解成
A
=
L
D
L
T
A=LDL^T
A=LDLT,其中L是单位下三角矩阵(单位下三角矩阵的对角线右上方的系数全部为零,左下方的系数全为一),D是对角均为正数的对角矩阵。把这一分解叫做LDL分解,是Cholesky分解的变形。具体先不展开了,可以参考[5],以及其中的参考代码。
分块Cholesky分解
这一部分其实原理和上面是一样的,只是
l
11
l_{11}
l11也是一个矩阵块来做,这样算法就是在一个矩阵块粒度递归往下做,效率上快很多,可以比较容易利用到GPU计算加速。我主要参考了文献[6]的内容,只做一个笔记用。
首先,把
A
A
A矩阵分块,其中
A
11
A_{11}
A11是一个
r
×
r
r\times r
r×r方阵,
r
r
r是我们设定的可以采用直接cholesky分解算法求解的块大小。
类似的,位面可以得到:
其中公式(8)其实我们一般计算的是线性方程组:
S ⋅ L 11 T = B S\cdot L_{11}^T = B S⋅L11T=B
其中 L 11 T L_{11}^T L11T是一个上三角形,因此,我们可以比较容易求出S(先直接可以求出S的第一列,然后是第二列,以此类推)。可以直接调用各种BLAS库的trsm函数求解[7]。
公式(9)又是一个Cholesky分解,一直递归下去采用分块Cholesky分解,直到 A ^ − S ⋅ S T \hat{A}-S\cdot S^T A^−S⋅ST的size小于等于 r × r r\times r r×r就不在分解了,最后采用一次直接Cholesky分解。
示意图以及流程总结如下:
到这里就把Cholesky分解的计算方法讲清楚了。
Cholesky分解应用于求解线性方程组,以及矩阵求逆
如果
A
A
A为对称正定矩阵,现在要求解线性方程组
A
X
=
B
AX=B
AX=B:
步骤:
- 求 A A A的Cholesky分解,得到 A = L L T A=LL^T A=LLT
- 把 A X AX AX看成 L ( L T X ) = L Y L(L^T X)=LY L(LTX)=LY,求 L Y = B LY = B LY=B,得到 Y Y Y
- 求 L T X = Y L^T X = Y LTX=Y,得到 X X X
这样就求出了线性方程组 A X = B AX=B AX=B的解 X X X。
类似的,如果
A
A
A为对称正定矩阵,我们要求
A
−
1
A^{-1}
A−1,我们实际上只要求解线性方程组
A
X
=
I
AX=I
AX=I:
步骤:
- 求 A A A的Cholesky分解,得到 A = L L T A=LL^T A=LLT
- 把 A X AX AX看成 L ( L T X ) = L Y L(L^T X)=LY L(LTX)=LY,求 L Y = I LY = I LY=I,得到 Y Y Y
- 求 L T X = Y L^T X = Y LTX=Y,得到 X X X
这样就求出了逆矩阵 X = A − 1 X= A^{-1} X=A−1。
参考资料
[1] https://www.jiqizhixin.com/articles/2019-03-05-8
[2] https://www.cnblogs.com/marsggbo/p/11461155.html
[3] https://baike.baidu.com/item/%E6%AD%A3%E5%AE%9A%E7%9F%A9%E9%98%B5/11030459?fr=aladdin
[4] https://www.qiujiawei.com/linear-algebra-11/
[5] https://blog.csdn.net/ACdreamers/article/details/44656847
[6] 使用GPU加速计算矩阵的Cholesky分解,沈 聪,高火涛
[7] https://blog.csdn.net/zb1165048017/article/details/70207812
更多推荐
所有评论(0)