LU分解的实现

LU分解是将矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积。矩阵可以不是NxN的矩阵

一个可逆矩阵可以进行LU分解当且仅当它的所有子式都非零。如果要求其中的L矩阵（或U矩阵）为单位三角矩阵，那么分解是唯一的。同理可知，矩阵的LDU可分解条件也相同，并且总是唯一的。

即使矩阵不可逆，LU仍然可能存在。实际上，如果一个秩为k的矩阵的前k个顺序主子式不为零，那么它就可以进行LU分解，但反之则不然。

目前，在任意域上一个方块矩阵可进行LU分解的充要条件已经被发现，这些充要条件可以用某些特定子矩阵的秩表示。用高斯消元法来得到LU分解的算法也可以扩张到任意域上。

任意矩阵A（不仅仅是方块矩阵）都可以进行LUP分解。其中的L和U矩阵是方阵，P矩阵则与A形状一样。

这里给出高斯消元法的思想

Matrix A (M x N)
for(column index i from 1 to N){
select max(A[i][i...M]), swap to row i //第i列中从第i行到第N行绝对值最大值元素的行作为第i行
if(A[i][i] is not zero){
for(row index j from i + 1 to M){
A[j][i] /= A[i][i]
for(column index k from i + 1 to N){
A[j][k] -= A[j][i] * A[i][k]
}
}
}
}

L如下

1, 	0, 	0, ... 	, 	0
A[21], 	1,	0, ...	,	0
A[31],	A[32],	1, ...	,	0
A[41],	A[42],	A[43], ...,	0
...	...
A[m1],	A[m2],	A[m3], ...,	1

U如下

A[11], 	A[12], 	A[13], ..., 	A[1n]
0, 	A[22],	A[23], ...,	A[2n]
0,	0,	A[33], ...,	A[3n]
0,	0,	0, ...	,	A[4n]
...	...
0,	0,	0, ...	,	A[mn]

下面给出实现，摘自JAMA

// Initialize.
LU = A.getArrayCopy();
m = A.getRowDimension();
n = A.getColumnDimension();
piv = new int[m];
for (int i = 0; i < m; i++) {
piv[i] = i;
}
pivsign = 1;
// Main loop.
for (int k = 0; k < n; k++) {
// Find pivot.
int p = k;
for (int i = k+1; i < m; i++) {
if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) {
p = i;
}
}
// Exchange if necessary.
if (p != k) {
for (int j = 0; j < n; j++) {
double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t;
}
int t = piv[p]; piv[p] = piv[k]; piv[k] = t;
pivsign = -pivsign;
}
// Compute multipliers and eliminate k-th column.
if (LU[k][k] != 0.0) {
for (int i = k+1; i < m; i++) {
LU[i][k] /= LU[k][k];
for (int j = k+1; j < n; j++) {
LU[i][j] -= LU[i][k]*LU[k][j];
}
}
}
}

由于Java对点乘支持的问题，上述方法在JAMA中并没有使用，而是使用了Crout/Doolittle算法，描述如下：







根据上述写的代码如下

for(int a = 0; a < length; a++){
for(int b = 0; b < length; b++){
double sum = 0.0;
if(a <= b){
for(int i = 0; i < a; i++){
sum += l[a][i] * u[i][b];
}
u[a][b] = matrix[a][b] - sum;
}else{
for(int i = 0; i < b; i++){
sum += l[a][i] * u[i][b];
}
l[a][b] = (matrix[a][b] - sum) / u[b][b];
}
}
}