基于Matlab 的 n阶非奇异方阵的LU分解实现

基于Matlab 的n阶方阵的LDU分解实现
1.引言
矩阵的LDU分解是“矩阵理论与方法”课程中非常重要的一部分。

LDU分解在实际工程应用中也非常广泛。

LDU分解可以将一个矩阵分解为一个下三角矩阵和一个对角矩阵和一个上三角矩阵的乘积。

LDU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式。

将系数矩阵A转变成等价两个矩阵L和U和对角矩阵的乘积，其中L和U 分别是下三角和上三角矩阵，D为对角矩阵。

当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LDU。

即：
Matlab 是很好的处理矩阵的工具。

它的功能非常强大，包括创建矩阵，对矩阵求逆，转置等操作非常简单，使其成为图像处理，信号分析等领域常用的工具。

Matlab 官方已经包括了对非奇异矩阵的LU分解函数[L,U]=lu(A)，为了加深对矩阵分解的理解，本文不采用Matlab 官方的LU分解函数对矩阵A进行LDU 分解，而是根据理论推导和编程实现LDU分解。

2.程序设计
2.1.输入合法检验
LU分解需要被分解矩阵A满足如下条件：
1)矩阵A为方阵
2)A的顺序主子式
全故LU分解需先检验A为n阶方阵，然后检验A的n-1个顺序主子式D
k 不为0，才可进行LU分解。

而检验主子式D
可以在n-1次循环LU分解中
k
进行，故先检验矩阵是否为方阵。

代码如下
2.2. n-1次循环LDU 分解
LDU 分解本质上是高斯消元法。

实质上是将A 通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。

从下至上地对矩阵A 做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L 矩阵，它也是一个单位下三角矩阵。

LDU 分解主要分为两步：
1根据高斯消元法对A i ()消元，消元矩阵为L i +1-1；2计算L i +1
-1A i ()=A i +1()以产生下一步迭代的A i ()。

2.2.1. 根据
构造L j ,L j -1(j =i +1)
高斯消元A i ()，使A i ()第i+1列从第i+2行至n 行都为0。

构造消元矩阵L j ,L j -1。

首先判断是否为0，为0则无法继续分解，退出；否则继续。

代码如下
2.2.2. 计算L i +1-1A i ()=A i +1()
D
k
计算L
i+1
-1A i()=A i+1()并储存覆盖A i()，并计算下次循环的主子式。

代码如下
2.3.计算并返回L和A n-1()为结果
通过n-1次循环累乘L
i 得到L矩阵，并且根据最后得到的A n-1
()矩阵分解出
D和U矩阵。

代码如下3.验证结果
3.1.可LDU分解的矩阵验证
构造可以LDU分解的矩阵A D
k
调用MyLU函数，分解A矩阵。

计算L*D*U验证是否为原矩阵A（方法下同）
3.2.输入矩阵A是否为空验证
输入矩阵A 中没有元素，输入不合法。

3.3. 输入矩阵A 为方阵验证
输入矩阵A 维度不合法。

3.4. 顺序主子式出现0的错误验证
因为计算()0
1
1*L A -时出现()1
22
0a =，导致顺序主子式()()0
1
211220a a ∆==。

该情况为不合法输入。

至此，函数功能和合法性检查全部验证完毕。

4.心得体会
通过“矩阵理论与方法”的理论指导和Matlab编程的实践经验，我基本掌握了矩阵分解中的LDU分解的推导过程和算法步骤。

熟练掌握LDU分解，对今后研究LU分解、Doolittle分解、Crout分解、QR分解等矩阵分解方法的实现有非常大的帮助。

对LDU分解的推导过程进行步骤分解和归纳，我将N维矩阵的LDU 分解归纳总结为n-1次循环，每步循环进行n-i次元素除法（计算高斯消元系数），2次N维矩阵乘法（计算L和()1i A+）和1次元素乘法（计算顺序主子式的值）。

按乘法计算，即时间复杂度为
()()
()()
34
121
O n n i n O n
--++≈。

因为计算过程中
需要
1
,
j j
L L-
辅助计算矩阵和L i+1-1A i
()=A i+1()
，即空间复杂度为()2
O n。

该LDU分解设计还不够快速，占用空间相对较多，是以后改进的方向。

5.附录：程序源码
function [L,D,U] = MyLU(A)
%check validity
if(isempty(A))%check A if is empty
error('A is empty!');
end
[N,D] = size(A);
if(N ~= D)%check A if is square
error('A is not a square!');
end
%LDU decomposing
L = eye(N);
det_k = A(1,1);
for step_n = 1:N-1
if (det_k == 0)% check if Sequence principal minor appears 0
error('Sequence principal minor is 0');
end
Li = eye(N);
Li_inv = eye(N);
for step_row = step_n+1:N
mod=A(step_row,step_n)/A(step_n,step_n);
Li(step_row,step_n) = mod;
Li_inv(step_row,step_n) = -mod;
end
A = Li_inv * A;
det_k = det_k * A(step_n+1,step_n+1);
L = L * Li;
end
L = L;
D = eye(N);
U = zeros(N);
for i=1:N
if(A(i,i) ~= 0)
D(i,i) = A(i,i);
U(i,:) = A(i,:)/D(i,i);
else
D(i,i) = 0;
U(i,i) = 1;
end
end
end。