看其他篇章到目錄選擇。
補充上一次的矩陣知識，這次主要講講矩陣的一些分解運算——Matrix Decomposition：
矩陣分解主要有三種方式：LU分解，QR分解和奇異值分解。當然在Math的linear包中提供了對應的接口有CholeskyDecomposition、EigenDecomposition、LUDecomposition、QRDecomposition和SingularValueDecomposition這5種分解方式。
每一個接口定義對應著一個***DecompositionImpl的實現類。
先來看看LU分解。
LU分解是矩陣分解的一種，可以將一個矩陣分解為一個下三角矩陣和一個上三角矩陣的乘積（有時是它們和一個置換矩陣的乘積）。LU分解主要應用在數值分析中，用來解線性方程、求逆矩陣或計算行列式。
Math庫中的LU分解主要是LUP分解，即針對n*n方陣A，找出三個n*n的矩陣L、U和P，滿足PA=LU。其中L是一個單位下三角矩陣，U是一個上三角矩陣，P是一個置換矩陣。非奇異的矩陣（可逆）都有這樣一種分解（可證明）。LUP分解的計算方法就是高斯消元。具體的算法見算導第28章或者參看我列出的參考資料。
Commons Math中的實現LUP分解是這樣的，主要在LUDecompositionImpl的構造方法中。見源碼：
1
public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold)
2
throws NonSquareMatrixException 
{
3
4
if (!matrix.isSquare()) {
5 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
6
}
7
8 final int m = matrix.getColumnDimension();
9 lu = matrix.getData();
10 pivot = new int[m];
11 cachedL = null;
12 cachedU = null;
13 cachedP = null;
14
15 // Initialize permutation array and parity
16 for (int row = 0; row < m; row++) {
17 pivot[row] = row;
18 }
19 even = true;
20 singular = false;
21
22 // Loop over columns
23 for (int col = 0; col < m; col++) {
24
25 double sum = 0;
26
27 // upper
28 for (int row = 0; row < col; row++) {
29 final double[] luRow = lu[row];
30 sum = luRow[col];
31 for (int i = 0; i < row; i++) {
32 sum -= luRow[i] * lu[i][col];
33 }
34 luRow[col] = sum;
35 }
36
37 // lower
38 int max = col; // permutation row
39 double largest = Double.NEGATIVE_INFINITY;
40 for (int row = col; row < m; row++) {
41 final double[] luRow = lu[row];
42 sum = luRow[col];
43 for (int i = 0; i < col; i++) {
44 sum -= luRow[i] * lu[i][col];
45 }
46 luRow[col] = sum;
47
48 // maintain best permutation choice
49 if (Math.abs(sum) > largest) {
50 largest = Math.abs(sum);
51 max = row;
52 }
53 }
54
55 // Singularity check
56 if (Math.abs(lu[max][col]) < singularityThreshold) {
57 singular = true;
58 return;
59 }
60
61 // Pivot if necessary
62 if (max != col) {
63 double tmp = 0;
64 final double[] luMax = lu[max];
65 final double[] luCol = lu[col];
66 for (int i = 0; i < m; i++) {
67 tmp = luMax[i];
68 luMax[i] = luCol[i];
69 luCol[i] = tmp;
70 }
71 int temp = pivot[max];
72 pivot[max] = pivot[col];
73 pivot[col] = temp;
74 even = !even;
75 }
76
77 // Divide the lower elements by the "winning" diagonal elt.
78 final double luDiag = lu[col][col];
79 for (int row = col + 1; row < m; row++) {
80 lu[row][col] /= luDiag;
81 }
82 }
83
84
}
85
測試代碼示例：
1/** *//**
2 *
3 */
4package algorithm.math;
5
6import org.apache.commons.math.linear.ArrayRealVector;
7import org.apache.commons.math.linear.DecompositionSolver;
8import org.apache.commons.math.linear.LUDecomposition;
9import org.apache.commons.math.linear.LUDecompositionImpl;
10import org.apache.commons.math.linear.MatrixUtils;
11import org.apache.commons.math.linear.RealMatrix;
12
13/** *//**
14 * @author Jia Yu
15 * @date 2010-11-20
16 */
17public class MatrixDecompositionTest {
18
19 /** *//**
20 * @param args
21 */
22 public static void main(String[] args) {
23 // TODO Auto-generated method stub
24 matrixDecomposite();
25 }
26
27 private static void matrixDecomposite() {
28 // TODO Auto-generated method stub
29 double[][] testData = {
30 { 1.0, 2.0, 3.0},
31 { 2.0, 5.0, 3.0},
32 { 1.0, 0.0, 8.0}
33 };
34
35 RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
36 //LUP decomposition
37 LUDecomposition LU = new LUDecompositionImpl(matrix);
38 RealMatrix l = LU.getL();
39 RealMatrix u = LU.getU();
40 RealMatrix p = LU.getP();
41 System.out.println("L is : "+l);
42 System.out.println("U is : "+u);
43 System.out.println("P is : "+p);
44 System.out.println("PA is "+(p.multiply(matrix)));
45 System.out.println("LU is "+(l.multiply(u)));
46 System.out.println("PA = LU is "+p.multiply(matrix).equals(l.multiply(u)));
47 //matrix singular
48 System.out.println("LU is not singular : "+LU.getSolver().isNonSingular());
49 //matrix determinant
50 System.out.println("matrix determinant is : "+LU.getDeterminant());
51 //matrix solver
52 RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
53 { 1, 0 }, { 2, -5 }, { 3, 1 }
54 });
55 DecompositionSolver solver = LU.getSolver();
56 System.out.println("solve Ax = b (when b is matrix) is x = "+solver.solve(b));
57 System.out.println("solve Ax = b (when b is vector) is x = "+new ArrayRealVector(solver.solve(b.getColumn(0))));
58 //matrix inverse
59 System.out.println("matrix inverse is "+solver.getInverse());
60 }
61
62}
63
輸出結果：
L is : Array2DRowRealMatrix{{1.0,0.0,0.0},{0.5,1.0,0.0},{0.5,0.2,1.0}}
U is : Array2DRowRealMatrix{{2.0,5.0,3.0},{0.0,-2.5,6.5},{0.0,0.0,0.19999999999999996}}
P is : Array2DRowRealMatrix{{0.0,1.0,0.0},{0.0,0.0,1.0},{1.0,0.0,0.0}}
PA is Array2DRowRealMatrix{{2.0,5.0,3.0},{1.0,0.0,8.0},{1.0,2.0,3.0}}
LU is Array2DRowRealMatrix{{2.0,5.0,3.0},{1.0,0.0,8.0},{1.0,2.0,3.0}}
PA = LU is true
LU is not singular : true
matrix determinant is : -0.9999999999999998
solve Ax = b (when b is matrix) is x = Array2DRowRealMatrix{{19.000000000000004,-71.00000000000001},{-6.000000000000002,22.000000000000007},{-2.0000000000000004,9.000000000000002}}
solve Ax = b (when b is vector) is x = {19; -6; -2}
matrix inverse is Array2DRowRealMatrix{{-40.00000000000001,16.000000000000004,9.000000000000002},{13.000000000000004,-5.000000000000002,-3.000000000000001},{5.000000000000001,-2.0000000000000004,-1.0000000000000002}}
對于其中求解Ax=b的算法，主要是一個正向替換和逆向替換的過程。這里簡單講一下過程，要求Ax=b的解，對兩邊同時乘以P，得到等價的PAx=Pb，通過LUP分解即LUx=Pb。定義y=Ux，正向替換：Ly=Pb，得到y，再逆向替換Ux=y，得到x。過程其實就是Ax=(P^-1)LUx=(P^-1)Ly=(P^-1)Pb=b。
矩陣求逆的運算等價于Ax=I的計算，I是對角單位矩陣。
QR分解和其他的分解對應的接口定義與LU分解是類似的。測試的代碼就不多貼了。有興趣的同學可以翻看一下源代碼，對于理解這一塊，源代碼還是在算法上有很大幫助的。哦，對了，補充一點，QR分解的實現是利用Householder算法的，想用其他算法練手的同學完全可以實現QRDecomposition這個接口并做實驗比對。
相關資料：
矩陣知識：http://zh.wikipedia.org/zh/%E7%9F%A9%E9%98%B5
矩陣分解：http://zh.wikipedia.org/zh-cn/%E7%9F%A9%E9%98%B5%E5%88%86%E8%A7%A3
QR分解：http://zh.wikipedia.org/zh-cn/QR%E5%88%86%E8%A7%A3
LU分解：http://zh.wikipedia.org/zh-cn/LU%E5%88%86%E8%A7%A3
高斯消元法：http://zh.wikipedia.org/zh-cn/%E9%AB%98%E6%96%AF%E6%B6%88%E5%8E%BB%E6%B3%95
奇異值分解：http://zh.wikipedia.org/zh-cn/%E5%A5%87%E5%BC%82%E5%80%BC%E5%88%86%E8%A7%A3
Householder算法：http://zh.wikipedia.org/zh-cn/Householder%E5%8F%98%E6%8D%A2
Commons math包：http://commons.apache.org/math/index.html