看其他篇章到目录选择。
补充上一次的矩阵知识，这次主要讲讲矩阵的一些分解运算——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