老饼讲解-机器学习 机器学习 神经网络 深度学习
相关基础数学

矩阵LU分解(一):Doolittle分解

作者 : 老饼 发表日期 : 2023-02-06 05:59:28 更新日期 : 2023-02-07 08:53:55
本站原创文章,转载请说明来自《老饼讲解-机器学习》www.bbbdata.com




LU分解是矩阵常用的分解方法之一,它可以把方阵分解成上下三角矩阵的积

本文介绍LU分解中的Doolittle分解方法,包括它的推导及代码实现




   01. LU分解之Doolittle分解方法   



本节介绍什么是LU分解、Doolittle分解,并讲解Doolittle分解的分解方法



    LU分解与Doolittle分解    


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

 
 其中
L是下三角矩阵
U是上三角矩阵

Doolittle分解
Doolittle分解则是LU分解方法中的一种,
它的特点是L的对角元素全为1
即Doolittle分解把A分解如下
 
其中 
 
✍️补充:关于LU分解的类别
 LU分解从L、U的特点可以分为三种
 👉 Doolittle分解   :L对角元全为1(即单位下三角矩阵)                
  👉 Crout分解        :U对角元全为1(即单位上三角矩阵)                
  👉 Cholesky分解  :U是L的转置(要求A必须是正定对称方阵)       





    Doolittle分解的计算公式    


Doolittle分解对L、U的求解是
通过对U逐行、L逐列进行迭代式计算求得
 即算出U的第1行,再算L的第1列,  
然后算U的第2行,再算L的第2列....


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

 
 
 
特别地,U的第1行元素
  
 L的计算公式
在已算出U的前i行,L的(i-1)列时,
 L的第i列计算公式如下

 
 
 
特别地,L的第1列元素







   02. Doolittle分解原理与推导   


本节介绍上节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的第一行同样适用




    U的求解公式推导    


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

                        A的第i列为 L乘以U 的i列            
                            表示为连加的形式         
               U是上三角矩阵,k>i时为0     
 独立抽出第i项    

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

  

因此,在L的前(i-1)列和U的前 i 行已求出
即可根据上式求出L第i列的所有元素
👉补充:上述推导原理对于L的第一列同样适用






   03. Doolittle矩阵分解的python代码实现   



本节根据上述原理编写Doolittle矩阵分解的实现代码



      python代码实现Doolittle矩阵分解     


# -*- coding: utf-8 -*-
"""
Doolittle矩阵分解(LU分解的一种)
"""
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)
	print('\n-----需要分解的矩阵A-----')
	print('A=\n',A)

	print('\n--------分解的结果------')
	print('L=\n',L)
	print('\nU=\n',U)

	print('\n---------结果验证--------')
	print('\nL*U=\n',L@U)



     代码运行结果    


-----需要分解的矩阵A-----
A=
 [[1. 3. 6. 2.]
  [2. 5. 4. 3.]
  [5. 4. 3. 5.]
  [8. 2. 1. 5.]]

--------分解的结果------
L=
 [[ 1.     0.    0.         0.    ]
  [ 2.     1.    0.         0.    ]
  [ 5.    11.    1.         0.    ]
  [ 8.    22.    2.1147541  1.    ]

U=
 [[ 1.     3.     6.     2.        ]
  [ 0.    -1.    -8.    -1.        ]
  [ 0.     0.    61.     6.        ]
  [ 0.     0.     0.    -1.68852459]]

---------结果验证--------

L*U=
 [[1. 3. 6. 2.]
  [2. 5. 4. 3.]
  [5. 4. 3. 5.]
  [8. 2. 1. 5.]]





    相关参考    


 【1】 Doolittle分解法(LU分解)详细分析以及matlab的实现 https://blog.csdn.net/lol_IP/article/details/78491457 





 End 








联系老饼