本站原创文章，转载请说明来自《老饼讲解-机器学习》www.bbbdata.com
LU分解是矩阵常用的分解方法之一，它可以把方阵分解成上下三角矩阵的积
本文介绍LU分解中的Crout分解方法，包括它的推导及代码实现
本节介绍什么是LU分解、Crout分解，并讲解Crout分解的分解方法
LU分解
LU分解
LU分解是指将方阵A分解为下三角矩阵L与上三角矩阵U的乘积
即
其中
L是下三角矩阵
U是上三角矩阵
Crout分解
Crout分解则是LU分解方法中的一种，
它的特点是U的对角元素全为1
即Crout分解把A分解如下
其中
，
LU分解的类别
LU分解从L、U的特点可以分为三种
👉 Doolittle分解 ：L对角元全为1(即单位下三角矩阵)
👉 Crout分解 ：U对角元全为1(即单位上三角矩阵)
👉 Cholesky分解 ：U是L的转置(要求A必须是正定对称方阵)
Crout分解
Crout分解对L、U的求解是
通过对U逐行、L逐列进行迭代式计算求得
即算出U的第1行，再算L的第1列，
然后算U的第2行，再算L的第2列....
L的计算公式
在已算出L的前(i-1)列，U的(i-1)行时，
L的第ｉ列计算公式如下
特别地，L的第1列元素
U的计算公式
在已算出L的前ｉ列，U的前(i-1)行时，
U的第ｉ行计算公式如下
特别地，U的第1行元素
本节介绍上节Crout分解方法中的公式是如何推导出来的
L的求解
L第i列的求解，可由下式得到
A的第i列为 L乘以U 的i列
表示为连加的形式
U是上三角矩阵，k>i时为0
独立抽出第i项
这是因为Uii=1
根据上式可得到L第i列的求解公式：
因此，在L的前(i-1)列和U的前(i-1)行已求出
即可根据上式求出L第i列的所有元素，
👉补充：上述推导原理对于L的第一列同样适用
U的求解
U第i行的求解，可由下式得到
A的第i行为 L的i行乘以U
表示为连加的形式
这是因为L是下三角矩阵，k>i时为0
独立抽出第i项
根据上式可得到U第i行的求解公式：
因此，在L的前 i 列和U的前(i-1)行已求出
即可根据上式求出U第i行的所有元素
👉补充：上述推导原理对于U的第一行同样适用
本节根据上述原理编写Doolittle矩阵分解的实现代码
python代码实现Doolittle矩阵分解
# -*- coding: utf-8 -*-
"""
Crout矩阵分解(LU分解的一种)
"""
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)
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. 61. 0. ]
[8. -22. 129. -1.68852459]]
U=
[[1. 3. 6. 2. ]
[0. 1. 8. 1. ]
[0. 0. 1. 0.09836066]
[0. 0. 0. 1. ]]
---------结果验证--------
L*U=
[[1. 3. 6. 2.]
[2. 5. 4. 3.]
[5. 4. 3. 5.]
[8. 2. 1. 5.]]相关参考
【1】线性方程组的解法——Crout分解法 https://www.jianshu.com/p/6819775e80c5
End