【介绍】LU分解之：Cholesky分解

Cholesky分解属于LU分解方法中的一种，它将矩阵A分解为A=L*U，而它的特点是分解出来的U是L的转置，但它要求被分解的矩阵必须是正定矩阵，Cholesky分解过程相对简单，只需对L逐列计算就可以了，下面就让我们一起来看看Cholesky分解到底是怎么干的吧~！

一、Cholesky分解是什么

Cholesky分解则是LU分解方法中的一种，我们先来看看LU分解是什么。

LU分解是指将方阵A分解为下三角矩阵L与上三角矩阵U的乘积，即

 其中，L是下三角矩阵，U是上三角矩阵，L和U维度与A一致

而Cholesky分解是LU分解中的一种，分解出来的U是L的转置，但它要求被分解的矩阵必须是正定矩阵。

也就是，Cholesky分解把正定矩阵A分解如下：

其中 

   , 

总的来说， Cholesky分解就是将正定方阵A分解成下三角矩阵与它的转置的积。

二、Cholesky矩阵分解流程

Cholesky矩阵分解流程，就是L的求解流程，它对L进行逐列求解，即先求第1列，再求第2列，再求第3列....如此类推。那下面我们看看它每列是如何计算出来的。

每列的求解分为两步：

1. 先求对角元素

2. 再求非对角元素

假设前列已求出，按如下流程求第列：

Step-1. 先求解L第i列对角元素

 的求解公式如下：

 更形象地，分别如下：

可以看到，第个对角元素就是：根号(-L第i行前(i-1)个元素的平方和)。

Step-2. 再求解L第i列所有元素

在求得后，第 列的求解公式如下：

解说：即L第列为：(A的第i列- L前(i-1)列*各自第i行的元素)/ 

矩阵形式为：

特别地，第1列的计算公式为

公式就只有这么多，按照上面的公式，逐列计算L就可以了。

三、Cholesky矩阵分解-公式推导

好了， 上面一起来看看Cholesky分解的公式是怎么推导出来的。

3.1. Cholesky分解-对角元素的公式推导

对于L的对角元素，公式的推导过程如下：

              A_ii等于L的 i 行乘的 i 列          

         的 i 列就是L的第i行的转置 

                 i之后的元素是0,所以只要累加到i即可

根据上式，即可得到

3.2. Cholesky分解-非角元素的公式推导

假设前(i-1)列已求出，且也已经求得，现在求L第i列的所有元素，则有：

                                 A的i列为 L乘以的 i 列

                                          L与互为转置

                                A的i列为L每列乘以各自第i行元素后再相加

                                L的i行i列之后的元素为0，故不需i之后的列

                                        独立抽出第i列

根据上式，即可得到：

将上式写成矩阵形式，则为：

三、Cholesky矩阵分解-代码实现

好了，最后我们根据上面的原理，来用python实现一下Crout分解。

闲话少说，直接上代码，如下：

"""
 本代码用于展示如何实现Cholesky分解
 本代码来自《老饼讲解-机器学习》www.bbblearn.com
"""
import numpy as np

# Cholesky分解函数
def CholeskyDecompose(A):
    r,c = A.shape                                           # 获取矩阵的行数、列数
    L = np.zeros((r,r))                                     # 初始化L
    L[0,0] = np.sqrt(A[0,0])                                # 计算L11
    L[:,0] = A[:,0]/L[0,0]                                  # 计算L第1列的其余元素
    for i in range(1,r):                                    # 逐列计算
        L[i,i] = np.sqrt(A[i,i]-sum(L[i,:i+1]*L[i,:i+1]))   # 计算第i列的对角元素Lii
        L[:,i] = (A[:,i]-L[:,:i]@L[i,:i].T)/L[i,i]          # 计算第i列的非对角元素
    return L                                                # 返回结果

# 测试样例
if __name__ == "__main__":
    # 生成要分解的A
    np.random.seed(0)                                       # 设置随机种子
    rnd_mat = np.random.randint(0,10,(5,5))                 # 生成一个随机矩阵
    L_true = np.tril(rnd_mat,0)                             # 将上三角置0,即得到下三角矩阵
    A = L_true@L_true.T                                     # 将L*L^T组合成A
												            
    # Cholesky分解                                          
    L = CholeskyDecompose(A)                                # 调用函数对A进行Cholesky分解
    print('\n----生成的L-----:\n',L_true)                   # 打印生成的L
    print('\n----LL^T组合成的A-----:\n',A)                  # 打印由LL^T得到的A
    print('\n----从A分解得到的L-------:\n',L)               # 打印从A中分解得到的L

代码运行结果如下：

从结果中可以看到，我们先用合成了A，然后用Cholesky分解A，就能分解出原始的L，说明程序是有效的。

总结

总的来说， Cholesky分解可以将正定方阵A分解成，而L是下三角矩阵。分解过程相对简单，只需要逐列求解L就可以了。