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

Crout分解(克劳特分解)是LU分解方法中的一种，它将矩阵A分解为A=L*U，而它的特点是分解出来的U，对角线元素全部都为1，Crout分解过程相对简单，只需对U逐行、对L逐列进行交替计算就可以了，下面就让我们一起来看看Crout分解到底是怎么干的吧~！

一、LU分解之Crout分解

1.1. Crout分解是什么

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

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

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

而Crout分解是LU分解中的一种，它的特点是U的对角元素全为1。

也就是，Crout分解把A分解如下：

其中 

 ，

总的来说， Crout分解就是将方阵A分解成下三角矩阵L和上三角矩阵U，而L的对角线元素全部为1。

1.2. 如何进行Crout分解

下面我们来说说Crout分解是如何求出L、U的。

它主要是通过"对U逐行、对L逐列"进行迭代式计算求得的：

如图，也就是先算出U的第1行，再算L的第1列，然后算U的第2行，再算L的第2列....如此类推。

好了，我们再来看看Crout分解中，U的每行、L的每列分别是怎么求出来的。

U的计算公式

在已算出L的前ｉ列，U的前(ｉ-1)行时，U的第ｉ行计算公式如下：

特别地，U的第1行元素 

L的计算公式

在已算出L的前(ｉ-1)列，U的(ｉ-1)行时， L的第ｉ列计算公式如下：

特别地，L的第1列元素

Crout分解计算起来还是挺方便的，按照上面的顺序、计算公式进行计算就可以了。

二、Crout分解-公式推导

上面说了Crout分解的计算方法了，那么为什么这样计算就能得到Crout分解呢？为了消除大家的好奇心，下面我们来说说它的原理和推导过程。

2.1. Crout分解-U的公式推导

U第ｉ行的求解，可由下式得到：

                      * A的第i行为 L的i行乘以U    

      * 表示为连加的形式 

                                * 这是因为L是下三角矩阵，k>i时为0 

                         * 独立抽出第i项

根据上式可得到U第ｉ行的求解公式为：

因此，在L的前ｉ列和U的前(ｉ-1)行已求出，即可根据上式求出U第i行的所有元素。

特别地，上述推导原理对于U的第一行同样适用。

2.2. Crout分解-L的公式推导

L第ｉ列的求解，可由下式得到： 

                      * A的第i列为 L乘以U 的i列 

       * 表示为连加的形式

                         * U是上三角矩阵，k>i时为0    

                             * 独立抽出第i项 

                                             * 这是因为Uii=1  

根据上式可得到L第ｉ列的求解公式为： 

因此，在L的前(ｉ-1)列和U的前(ｉ-1)行已求出，即可根据上式求出L第ｉ列的所有元素。

特别地，上述推导原理对于L的第一列同样适用。

三、Crout分解-代码实现

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

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

"""
 本代码用于展示如何实现Crout矩阵分解(LU分解的一种)
 本代码来自《老饼讲解-机器学习》www.bbblearn.com
"""
import numpy as np 

# Crout矩阵分解函数
def CroutDecompose(A):                                               
    size   = A.shape[0]                                                # 方阵的大小
    L      = np.zeros((size,size))                                     # 初始化L全为0
    U      = np.identity(size)                                         # 初始化U为单位矩阵
    L[:,0] = A[:,0].copy()                                             # 计算L的第一列
    U[0,:] = A[0,:]/L[0,0]                                             # 计算U的第一行
												                      
    for i in range(1,size):                                            # 逐列、行计算L和U
       L[:,i] =  A[:,i] - L[:,:i]@U[:i,i]                              # 计算L的第i列
       U[i,:] = (A[i,:] - L[i,:i]@U[:i,:])/L[i,i]                      # 计算U的第i行
	   
	# 为避免产生数值问题,将L上三角元素、U的下三角元素置0
    L=np.tril(L)                                                       # L只取下三角
    U=np.triu(U)                                                       # U只取上三角
    return L,U

# 测试样例
if __name__ == "__main__":
	A = np.array([[1.,2.,5,8],[3.,5.,4,2],[6.,4,3,1],[2.,3,5,5]]).T    # 生成需要分解的矩阵A
	L,U= CroutDecompose(A)                                             # Crout矩阵分解
	print('\n-----需要分解的矩阵A-----')                               # 打印标题 
	print('A=\n',A)                                                    # 打印矩阵A
	print('\n--------分解的结果------')                                # 打印标题
	print('L=\n',L)                                                    # 打印分解的L
	print('\nU=\n',U)                                                  # 打印分解的U
	print('\n---------结果验证--------')                               # 打印标题
	print('L*U=\n',L@U)                                                # 打印分解的LU乘积

代码运行结果如下：

从运行结果中可以看到，从A中分解出下三角矩阵L和上三角矩阵U，将L与U相乘后，的确是等于A，说明程序进行的Crout分解是无误的。

总结

总的来说， Crout分解就是将方阵A分解成下三角方阵L和上三角方阵U，而U的对角线元素全部为1。分解过程相对简单，只需要对U逐行、L逐列来进行计算就可以了。