LU分解的實現
阿新 • • 發佈:2019-01-07
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];
}
}
}