1. 程式人生 > >LU分解的實現

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];
				}
			}
		}