LU三角分解
LU分解(LU Factorization)是矩阵分解的一种,可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积(有时是它们和一个置换矩阵的乘积)。LU分解主要应用在数值分析中,用来解线性方程、求反矩阵或计算行列式。
简介
将系数矩阵A转变成等价两个矩阵L和U的乘积 ,其中L和U分别是单位下三角矩阵和上三角矩阵。当A的所有顺序主子式都不为0时,矩阵A可以分解为A=LU。其中L是下三角矩阵,U是上三角矩阵。
算法
LU分解在本质上是高斯消元法的一种表达形式。**实质上是将A通过初等行变换变成一个上三角矩阵,其变换矩阵就是一个单位下三角矩阵。**这正是所谓的杜尔里特算法(Doolittle algorithm):从下至上地对矩阵A做初等行变换,将对角线左下方的元素变成零,然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵,这一系列单位下三角矩阵的乘积的逆就是L矩阵,它也是一个单位下三角矩阵。这类算法的复杂度一般在(三分之二的n三次方) 左右
Java算法:
import java.util.Arrays;
/**
* 矩阵的直接三角分解 ,调用示例:
*
* DirectDecomposition dd = new DirectDecomposition(data);//data为一个二维double数组,代替一个矩阵
*
* double[][] l = dd.getL();//获取L
*
* double[][] u = dd.getU();//获取U
*
* @author jungle
*/
public class DoolittleDecomposition {
private double[][] data;
private double[][] l;
private double[][] u;
private int n;
/**
* 创建一个n阶的矩阵
*
* @param n
*/
public DoolittleDecomposition(double[][] data) {
if (data == null || data.length == 0 || data.length != data[0].length) {
throw new RuntimeException("不是一个方阵");
}
this.data = data;
n = data.length;
l = new double[n][n];
u = new double[n][n];
countLU();
}
protected void countLU() {
for (int i = 0; i < n; i++) {// 第一步,计算L的第一列和U的第一行:U1i=A1i,Li1=Ai1/U11
u[0][i] = data[0][i];
l[i][0] = data[i][0] / u[0][0];
}
//计算U的第r行,L的第r列元素
//uri=ari-Σ(k=1->r-1)lrkuki (i=r,r+1,...,n)
//lir=air-Σ(k=1->r-1)likukr (i=r+1,r+2,...,n且r≠n)
for (int r = 1; r < n; r++) {
for (int i = r; i < n; i++) {
u[r][i] = data[r][i] - sumLrkUki(r, i);
if(i==r) l[r][r]=1;
else if(r==n) l[n][n]=1;
else l[i][r] = (data[i][r] - sumLikUkr(r, i)) / u[r][r];
}
}
}
/**
* 求和:Lrk*Uki 对k求和:1<=k<=r-1
*
* @param r
* @param i
* @return
*/
private double sumLrkUki(int r, int i) {
double re = 0.0;
for (int k = 0; k < r; k++) {
re += l[r][k] * u[k][i];
}
return re;
}
private double sumLikUkr(int r, int i) {
double re = 0.0;
for (int k = 0; k < r; k++) {
re += l[i][k] * u[k][r];
}
return re;
}
public double[][] getData() {
return data;
}
public double[][] getL() {
return l;
}
public double[][] getU() {
return u;
}
public static void main(String[] args) {
double[][] data= {
{1,2,6},
{2,5,15},
{6,15,46},
};
DoolittleDecomposition dd = new DoolittleDecomposition(data);
double[][] l = dd.getL();
double[][] u = dd.getU();
int n = l.length;
System.out.println("L阵:");
for (int i = 0; i < n; i++) {
System.out.println(Arrays.toString(u[i]));
}
System.out.println("---------------------");
System.out.println("U阵:");
for (int i = 0; i < n; i++) {
System.out.println(Arrays.toString(l[i]));
}
}
}
MATLAB算法:
%LU分解法求解Ax=b,假定A矩阵可进行LU分解以及对角线元素均不为0
function [x] = Dool(A,b)
n=length(A);A(2:n,1)=A(2:n,1)/A(1,1);
for t=2:n-1 %进行LU分解
A(t,t:n)=A(t,t:n)-A(t,1:t-1)*A(1:t-1,t:n);
A(t+1:n,t)=(A(t+1:n,t)-A(t+1:n,1:t-1)*A(1:t-1,t))/A(t,t);
end
A(n,n)=A(n,n)-A(n,1:n-1)*A(1:n-1,n);
L=tril(A,-1)+eye(n);U=triu(A); %矩阵A分解出L和U
for t=2:n %解Lx=b
b(t)=b(t)-L(t,1:t-1)*b(1:t-1);
end
b(n)=b(n)/U(n,n); %解Ux=b
for t=1:n-1;
k=n-t;b(k)=(b(k)-U(k,k+1:n)*b(k+1:n))/U(k,k);
end
x=b; %方程Ax=b的解即为x