写一篇关于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 XTAX0 恒成立,则矩阵 A A A 是一个半正定矩阵。
    在这里插入图片描述

直接Cholesky分解

定理——若 A ∈ R n × n A \in R^{n\times n} ARn×n是对称正定矩阵,则存在一个对角元全为正数的下三角矩阵 L ∈ R n × n L \in R^{n\times n} LRn×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=A22L21L21T

其中 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=A22L21L21T

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 SL11T=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^SST的size小于等于 r × r r\times r r×r就不在分解了,最后采用一次直接Cholesky分解。

示意图以及流程总结如下:
在这里插入图片描述
在这里插入图片描述

到这里就把Cholesky分解的计算方法讲清楚了。

Cholesky分解应用于求解线性方程组,以及矩阵求逆

如果 A A A为对称正定矩阵,现在要求解线性方程组 A X = B AX=B AX=B
步骤:

  1. A A A的Cholesky分解,得到 A = L L T A=LL^T A=LLT
  2. 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
  3. 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} A1,我们实际上只要求解线性方程组 A X = I AX=I AX=I
步骤:

  1. A A A的Cholesky分解,得到 A = L L T A=LL^T A=LLT
  2. 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
  3. L T X = Y L^T X = Y LTX=Y,得到 X X X

这样就求出了逆矩阵 X = A − 1 X= A^{-1} X=A1

参考资料

[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

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐