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

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

一、LU分解之Doolittle分解

1.1. Doolittle分解是什么

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

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

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

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

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

其中 

 ，

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

1.2. 如何进行Doolittle分解

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

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

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

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

U的计算公式

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

特别地，Ｕ的第1行元素

L的计算公式

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

 特别地，L的第1列元素

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

二、Doolittle分解-公式推导

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

2.1. Doolittle分解-U的求解公式推导

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

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

       * 表示为连加的形式         

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

                        * 独立抽出第i个 

                                         * 这是因为Lii=1  

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

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

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

2.2. Doolittle分解-L的求解公式推导

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

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

         * 表示为连加的形式                    

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

      * 独立抽出第i项

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

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

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

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

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

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

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

# Doolittle矩阵分解函数
def DoolittleDecompose(A):
    size   = A.shape[0]                                                   # 方阵的大小
    L      = np.identity(size)                                            # 初始化L为单位矩阵
    U      = np.zeros((size,size))                                        # 初始化U全为0
    U[0,:] = A[0,:].copy()                                                # 计算U的第一行
    L[:,0] = A[:,0]/U[0,0]                                                # 计算L的第一列
												                         
    for i in range(1,size):                                               # 逐行、列计算U和L
       U[i,:] =  A[i,:] - L[i,:i]@U[:i,:]                                 # 计算U的第i行
       L[:,i] = (A[:,i] - L[:,:i]@U[:i,i])/U[i,i]                         # 计算L的第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 = DoolittleDecompose(A)                                           # 将A进行Doolittle分解
	print('\n-----需要分解的矩阵A-----')                                  # 打印标题
	print('A=\n',A)                                                       # 打印矩阵A

	print('\n--------分解的结果------')                                   # 打印标题
	print('L=\n',L)                                                       # 打印分解后得到的L
	print('\nU=\n',U)                                                     # 打印分解后得到的U

	print('\n---------结果验证--------')                                  # 打印标题
	print('\nL*U=\n',L@U)                                                 # 打印验证结果

代码运行结果如下：

可以看到，Doolittle分解将矩阵A分解为了上三角矩阵L和下三角矩阵U，两者的积就是A。

总结

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