[矩阵的三角分解系列四] 乔累斯基(Cholesky)分解公式

矩阵的三角分解是求解线性方程组常用的方法，包括LU分解，LDU分解，杜利特(Doolittle)分解，克劳特(Crout)分解，LLT(乔累斯基Cholesky)分解，LDLT(不带平方根乔累斯基)分解等，以及为了满足分解条件又加入行列变换的LPU分解，PLU分解，LUP分解，LDPU分解等。这里矩阵的三角分解系列教程主要是针对在学习三角分解时候的涉及到的一些细节，包括很多方法的来源和证明等，以及其中用到的一些矩阵操作的基础知识，主要包括：

• [矩阵的三角分解系列一] 高斯消元法
• [矩阵的三角分解系列二] LDU基本定理
• [矩阵的三角分解系列三] 杜利特/克劳特分解公式
• [矩阵的三角分解系列四] 乔累斯基(Cholesky)分解公式
• [矩阵的三角分解系列五] 三角分解中的行列变换
• [矩阵的三角分解系列六] Eigen中的三角分解

这个系列后面文章会用到前面文章的理论和技术，所以建议按照顺序查看。

简介

之前[矩阵的三角分解系列三] 杜利特/克劳特分解公式文章介绍了针对一般的$n$阶方阵$\mathbf{A}$的三角分解的方法，但对于一些特殊的公式在分解时候会满足某些特定的性质可以减少分解的计算，类似乔累斯基(Cholesky)分解这种。

LLT分解

如果$\mathbf{A}$是$n$阶对称正定矩阵，则存在
$$\mathbf{A}=\mathbf{L}{\mathbf{L}}^{\mathrm{T}},\text{\hspace{1em}(1)}$$
如果限定$\mathbf{L}$对角元素为正，则这个分解是唯一的，其中$\mathbf{L}$是下三角矩阵。

证明

因为$\mathbf{A}$是对称正定的，所以根据性质可知它的各阶顺序主子式均为正${\Delta }_k>0(k=1,2,\cdots ,n)$。所以由前面的LDU基本定理可知$\mathbf{A}$可以被唯一的分解为
$$\mathbf{A}=\tilde{\mathbf{L}}\mathbf{D}\tilde{\mathbf{U}}$$
又知$\tilde{\mathbf{L}}$是单位上三角矩阵，$\tilde{\mathbf{U}}$是单位下三角矩阵，$\mathbf{D}$是对角矩阵$\mathbf{D}=\mathrm{diag}(d_1,d_2,\cdots ,d_n)$。
又因为$\mathbf{A}$是对称矩阵，所以有
$$\mathbf{A}=\tilde{\mathbf{L}}\mathbf{D}\tilde{\mathbf{U}}={\mathbf{A}}^{\mathrm{T}}=(\tilde{\mathbf{L}}\mathbf{D}\tilde{\mathbf{U}}{)}^{\mathrm{T}}={\tilde{\mathbf{U}}}^{\mathrm{T}}{\mathbf{D}}^{\mathrm{T}}{\tilde{\mathbf{L}}}^{\mathrm{T}}$$
又知道分解是唯一的，所以
$$\tilde{\mathbf{L}}={\tilde{\mathbf{U}}}^{\mathrm{T}},\tilde{\mathbf{U}}={\tilde{\mathbf{L}}}^{\mathrm{T}},$$
所以$\mathbf{A}$的三角分解又可以写成
$$\mathbf{A}=\tilde{\mathbf{L}}\mathbf{D}{\tilde{\mathbf{L}}}^{\mathrm{T}}={\tilde{\mathbf{U}}}^{\mathrm{T}}\mathbf{D}\tilde{\mathbf{U}}$$
后面讲解以$\mathbf{A}={\tilde{\mathbf{U}}}^{\mathrm{T}}\mathbf{D}\tilde{\mathbf{U}}$为例，因为$|\tilde{\mathbf{U}}|=1\ne 0$，所以$\tilde{\mathbf{U}}$是可逆的，可得
$$\mathbf{D}=({\tilde{\mathbf{U}}}^{\mathrm{T}}{)}^{-1}\mathbf{A}{\tilde{\mathbf{U}}}^{-1}=({\tilde{\mathbf{U}}}^{-1}{)}^{\mathrm{T}}\mathbf{A}{\tilde{\mathbf{U}}}^{-1}$$
又因为$\mathbf{A}$是正定矩阵，又当$\mathbf{x}\ne 0$，明显${\tilde{\mathbf{U}}}^{-1}\mathbf{x}\ne 0$，所以
$${\mathbf{x}}^{\mathrm{T}}\mathbf{D}\mathbf{x}={\mathbf{x}}^{\mathrm{T}}({\tilde{\mathbf{U}}}^{-1}{)}^{\mathrm{T}}\mathbf{A}{\tilde{\mathbf{U}}}^{-1}\mathbf{x}=({\tilde{\mathbf{U}}}^{-1}\mathbf{x}{)}^{\mathrm{T}}\mathbf{A}{\tilde{\mathbf{U}}}^{-1}\mathbf{x}>0$$
所以$\mathbf{D}$也是正定矩阵，那么所有的一阶主子式(对角元素)$d_i>0(i=1,2,\cdots ,n)$，所以令
$${\mathbf{D}}^{\frac{1}{2}}=\mathrm{diag}(\sqrt{d_1},\sqrt{d_2},\cdots ,\sqrt{d_n}),$$
则$\mathbf{A}$又可以被唯一的分解为
$$\mathbf{A}=\tilde{\mathbf{L}}\mathbf{D}\tilde{\mathbf{U}}=\tilde{\mathbf{L}}{\mathbf{D}}^{\frac{1}{2}}{\mathbf{D}}^{\frac{1}{2}}{\tilde{\mathbf{L}}}^{\mathrm{T}}=(\tilde{\mathbf{L}}{\mathbf{D}}^{\frac{1}{2}})(\tilde{\mathbf{L}}{\mathbf{D}}^{\frac{1}{2}}{)}^{\mathrm{T}}=\mathbf{L}{\mathbf{L}}^{\mathrm{T}}$$
其中$\mathbf{L}=\tilde{\mathbf{L}}{\mathbf{D}}^{\frac{1}{2}}$对对角线元素全为正数$\sqrt{d_1},\sqrt{d_2},\cdots ,d_n$的下三角矩阵。

公式$(1)$证明完成！！!

公式$(1)$是对称正定矩阵$\mathbf{A}$的乔累斯基(Cholesky)分解，亦称为平方根分解。

具体解法

下面给出进行乔累斯基(Cholesky)分解具体分解公式和步骤，对于$n$阶对称正定矩阵$\mathbf{A}$，有分解式
$$\mathbf{A}=\mathbf{L}{\mathbf{L}}^{\mathrm{T}}$$
即
$$\left[\begin{array}{ccccc}{a}_{11}& \cdots & a_{1j}& \cdots & a_{1n}\\ ⋮& & ⋮& & ⋮\\ a_{i1}& \cdots & a_{ij}& \cdots & a_{in}\\ ⋮& & ⋮& & ⋮\\ a_{n1}& \cdots & a_{nj}& \cdots & a_{nn}\end{array}\right]=\left[\begin{array}{ccccc}{l}_{11}& & & & \\ ⋮& \ddots & & & \\ l_{i1}& \cdots & l_{ii}& & \\ ⋮& & & \ddots & \\ l_{n1}& \cdots & l_{ni}& \cdots & l_{nn}\end{array}\right]\left[\begin{array}{ccccc}{l}_{11}& \cdots & l_{j1}& \cdots & l_{n1}\\ & \ddots & & & ⋮\\ & & l_{jj}& \cdots & l_{nj}\\ & & & \ddots & ⋮\\ & & & & l_{nn}\end{array}\right]$$
由于$\mathbf{A}$对称，所以只考虑下三角元素，即$i\ge j$的情况，有
$$a_{ij}=\sum _{k=1}^j{l}_{ik}{l}_{jk}=\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk}+l_{ij}{l}_{jj},$$
即
$$l_{ij}=\left(a_{ij}-\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk}\right)/l_{jj},i\ge j,\text{\hspace{1em}(2)}$$
而且当$i=j$时，有
$$l_{ii}=\sqrt{a_{ii}-\sum _{k=1}^{i-1}{l}_{ik}^2}.\text{\hspace{1em}(3)}$$
这里与克劳特三角分解公式不同的是矩阵$\mathbf{U}={\mathbf{L}}^{\mathrm{T}}$，所以在求得$\mathbf{L}$的第$i$行元素之后，${\mathbf{L}}^{\mathrm{T}}$的第$i$列元素也已求出，所以计算量相当于克劳特分解的一半左右(当然对角线元素是都需要求一次)。

稳定性

同时又有
$$a_{ii}=\sum _{k=1}^i{l}_{ik}^2,i=1,2,\cdots ,n$$
所以
l i k 2 ≤ a i i ≤ max ⁡ 1 ≤ i ≤ n { a i i } l_{ik}^2 \leq a_{i i} \leq \max_{1 \leq i \leq n}\{ a_{i i} \} lik2​≤aii​≤1≤i≤nmax​{aii​}
于是
max ⁡ i , k { l i k 2 } ≤ max ⁡ 1 ≤ i ≤ n { a i i } , \max_{i,k}\{ l_{ik}^2 \} \leq \max_{1 \leq i \leq n}\{ a_{i i} \}, i,kmax​{lik2​}≤1≤i≤nmax​{aii​},
以上分析说明，分解过程中元素$l_{ik}$的数量级完全可以控制，从而计算过程是稳定的。

LDLT分解

然而，公式$(1)$这种$\mathbf{L}{\mathbf{L}}^{\mathrm{T}}$分解也还是存在很大的缺陷的，因为计算这种分解要进行$n$次开方运算，而绝大多数计算机上的开方运算是用子程序实现的(转化成对数计算)，这样不仅增加了计算量，而且还会扩大误差，甚至计算过程中有平方根号下出现负数的风险。为了避免上面的平方根运算，只需要在公式$(1)$中的$\mathbf{L}$和${\mathbf{L}}^{\mathrm{T}}$之间插入一个特殊的对角矩阵$\mathbf{D}$，就能达到预定的效果。
如果$\mathbf{A}$是$n$阶对称正定矩阵，则$\mathbf{A}$可以唯一的分解为
$$\mathbf{A}=\mathbf{L}\mathbf{D}{\mathbf{L}}^{\mathrm{T}},\text{\hspace{1em}(4)}$$
其中$\mathbf{L}$是下三角矩阵，其中$\mathbf{D}$是对角矩阵，且对角线元素是$\mathbf{L}$对角线元素的倒数，即
$$\mathbf{L}=\left[\begin{array}{cccc}{l}_{11}& & & \\ l_{21}& l_{22}& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & l_{nn}\end{array}\right],\mathbf{D}=\left[\begin{array}{cccc}\frac{1}{l_{11}}& & & \\ & \frac{1}{l_{22}}& & \\ & & \ddots & \\ & & & \frac{1}{l_{22}}\end{array}\right].$$

证明

利用特劳特分解公式以及矩阵对称性$a_{ij}=a=ji$就能推出公式$(4)$了。证明如下
$$\mathbf{A}=\left[\begin{array}{cccc}{l}_{11}& & & \\ l_{21}& l_{21}& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & l_{nn}\end{array}\right]\left[\begin{array}{cccc}1& u_{12}& \cdots & u_{1n}\\ & 1& \cdots & u_{2n}\\ & & \ddots & ⋮\\ & & & 1\end{array}\right]$$
其中
$$\begin{array}{rl}& \{\begin{array}{rl}{u}_{12}=& a_{12}/l_{11}=a_{21}/l_{11}=l_{21}/l_{11}\\ ⋮& \\ u_{1n}=& l_{n1}/l_{11}\end{array}\\ & \{\begin{array}{rl}{u}_{23}=& \left(a_{23}-l_{21}{u}_{13}\right)/l_{22}=\left(a_{32}-l_{21}{l}_{31}/l_{11}\right)/l_{22}\\ =& \left(a_{32}-l_{31}{u}_{12}\right)/l_{22}=l_{32}/l_{22}\\ ⋮& \\ u_{2n}=& l_{n2}/l_{22}\end{array}\end{array}$$
以此类推，$\mathbf{U}$的第$j$行元素有
$$\{\begin{array}{rl}& u_{j,j+1}=l_{j+1,j}/l_{jj},j=1,\cdots ,n-1,\\ & ⋮\\ & u_{jn}=l_{nj}/l_{jj}\end{array}$$
于是可以推导出$\mathbf{U}$为
$$\begin{array}{rl}\mathbf{U}& =\left[\begin{array}{ccccc}1& \cdots & \frac{l_{j1}}{l_{11}}& \cdots & \frac{l_{n1}}{l_{11}}\\ & \ddots & ⋮& & ⋮\\ & & 1& \cdots & \frac{l_{nj}}{l_{jj}}\\ & & & \ddots & ⋮\\ & & & & 1\end{array}\right]\\ & =\left[\begin{array}{ccccc}\frac{1}{l_{11}}& & & & \\ & \ddots & & & \\ & & \frac{1}{l_{jj}}& & \\ & & & \ddots & \\ & & & & \frac{1}{l_{nn}}\end{array}\right]\left[\begin{array}{ccccc}{l}_{11}& \cdots & l_{j1}& \cdots & l_{n1}\\ & \ddots & ⋮& & ⋮\\ & & l_{jj}& \cdots & l_{nj}\\ & & & \ddots & ⋮\\ & & & & l_{nn}\end{array}\right]=\mathbf{D}{\mathbf{L}}^{\mathrm{T}}\end{array}$$
所以有
$$\mathbf{L}\mathbf{D}{\mathbf{L}}^{\mathrm{T}}$$
公式$(15)$是对称正定矩阵$\mathbf{A}$的不带平方根乔累斯基(Cholesky)分解，亦称乔累斯基(Cholesky)分解的变形。

具体解法

下面给出进行不带平方根乔累斯基(Cholesky)分解具体分解公式和步骤，对于$n$阶对称正定矩阵$\mathbf{A}$，有分解式
$$\mathbf{L}\mathbf{D}{\mathbf{L}}^{\mathrm{T}}$$
即
$$\left[\begin{array}{ccccc}{a}_{11}& \cdots & a_{i1}& \cdots & a_{n1}\\ ⋮& & ⋮& & ⋮\\ a_{i1}& \cdots & a_{ii}& \cdots & a_{ni}\\ ⋮& & ⋮& & ⋮\\ a_{n1}& \cdots & a_{ni}& \cdots & a_{nn}\end{array}\right]=\left[\begin{array}{ccccc}{l}_{11}& & & & \\ ⋮& \ddots & & & \\ l_{i1}& \cdots & l_{ii}& & \\ ⋮& & ⋮& \ddots & \\ l_{n1}& \cdots & l_{ni}& \cdots & l_{nn}\end{array}\right]\left[\begin{array}{ccccc}1& \cdots & \frac{l_{j1}}{l_{11}}& \cdots & \frac{l_{n1}}{l_{11}}\\ & \ddots & ⋮& & ⋮\\ & & 1& \cdots & \frac{l_{nj}}{l_{jj}}\\ & & & \ddots & ⋮\\ & & & & 1\end{array}\right]$$
由于$\mathbf{A}$对称，所以只考虑下三角元素，即$i\ge j$的情况，有
$$a_{ij}=\sum _{k=1}^j{l}_{ik}(l_{jk}/l_{kk})=\sum _{k=1}^{j-1}{l}_{ik}(l_{jk}/l_{kk})+l_{ij},$$
即
$$l_{ij}=a_{ij}-\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk}/l_{kk},i\ge j,\text{\hspace{1em}(5)}$$
而且当$i=j$时，有
$$l_{ii}=a_{ii}-\sum _{k=1}^{j-1}{l}_{ik}^2/l_{kk},\text{\hspace{1em}(6)}$$
这样做分解计算就不会出现开方的运算，弥补了公式$(1)$平方根分解的缺陷。因此，$\mathbf{L}\mathbf{D}\mathbf{L}\mathbf{T}$是求解对称正定线性方程组最常用的一个分解公式。

例子

LLT分解

将下列矩阵进行乔累斯基(Cholesky)分解
$$\mathbf{A}=\left[\begin{array}{ccc}2& 1& 1\\ 1& 3& 2\\ 1& 2& 2\end{array}\right].$$
解：
根据公式$(13)$和公式$(14)$可得
$$\begin{array}{rl}{l}_{11}& =\sqrt{a_{11}}=\sqrt{2}\\ l_{21}& =a_{21}/l_{11}=\frac{\sqrt{2}}{2}\\ l_{22}& =\sqrt{a_{22}-l_{21}^2}=\frac{\sqrt{10}}{2}\\ l_{31}& =a_{31}/l_{11}=\frac{\sqrt{2}}{2}\\ l_{32}& =\left(a_{32}-l_{31}{l}_{21}\right)/l_{22}=\frac{3\sqrt{10}}{10}\\ l_{33}& =\sqrt{a_{33}-l_{31}^2-l_{32}^2}==\frac{\sqrt{15}}{5}\end{array}$$
则
$$\mathbf{A}=\mathbf{L}{\mathbf{L}}^{\mathrm{T}}=\left[\begin{array}{ccc}\sqrt{2}& & \\ \frac{\sqrt{2}}{2}& \frac{\sqrt{10}}{2}& \\ \frac{\sqrt{2}}{2}& \frac{3\sqrt{10}}{10}& \frac{\sqrt{15}}{5}\end{array}\right]\left[\begin{array}{ccc}\sqrt{2}& \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\\ & \frac{\sqrt{10}}{2}& \frac{3\sqrt{10}}{10}\\ & & \frac{\sqrt{15}}{5}\end{array}\right]$$

LDLT分解

将下列矩阵进行不带平方根乔累斯基(Cholesky)分解
$$\mathbf{A}=\left[\begin{array}{ccc}2& 1& 1\\ 1& 3& 2\\ 1& 2& 2\end{array}\right].$$
解：
根据公式$(16)$和公式$(17)$可得
$$\begin{array}{rl}{l}_{11}& =a_{11}=2\\ l_{21}& =a_{21}=1\\ l_{22}& =a_{22}-l_{21}^2/l_{11}=\frac{5}{2}\\ l_{31}& =a_{31}=1\\ l_{32}& =a_{32}-l_{31}{l}_{21}/l_{11}=\frac{3}{2}\\ l_{33}& =a_{33}-l_{31}^2/l_{11}-l_{32}^2/l_{22}==\frac{3}{5}\end{array}$$
则
$$\mathbf{A}=\mathbf{L}\mathbf{D}{\mathbf{L}}^{\mathrm{T}}=\left[\begin{array}{ccc}2& & \\ 1& \frac{5}{2}& \\ 1& \frac{3}{2}& \frac{3}{5}\end{array}\right]\left[\begin{array}{ccc}\frac{1}{2}& & \\ & \frac{2}{5}& \\ & & \frac{5}{3}\end{array}\right]\left[\begin{array}{ccc}2& 1& 1\\ & \frac{5}{2}& \frac{3}{2}\\ & & \frac{3}{5}\end{array}\right]$$

引用

【1】 矩阵论(第二版)