看其他篇章到目录选择。
补充上一次的矩阵知识,这次主要讲讲矩阵的一些分解运算——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![](/Images/OutliningIndicators/ExpandedBlockStart.gif)
throws NonSquareMatrixException
{
3![](/Images/OutliningIndicators/InBlock.gif)
4![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if (!matrix.isSquare())
{
5
throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
6
}
7![](/Images/OutliningIndicators/InBlock.gif)
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![](/Images/OutliningIndicators/InBlock.gif)
15
// Initialize permutation array and parity
16![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int row = 0; row < m; row++)
{
17
pivot[row] = row;
18
}
19
even = true;
20
singular = false;
21![](/Images/OutliningIndicators/InBlock.gif)
22
// Loop over columns
23![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int col = 0; col < m; col++)
{
24![](/Images/OutliningIndicators/InBlock.gif)
25
double sum = 0;
26![](/Images/OutliningIndicators/InBlock.gif)
27
// upper
28![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int row = 0; row < col; row++)
{
29
final double[] luRow = lu[row];
30
sum = luRow[col];
31![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int i = 0; i < row; i++)
{
32
sum -= luRow[i] * lu[i][col];
33
}
34
luRow[col] = sum;
35
}
36![](/Images/OutliningIndicators/InBlock.gif)
37
// lower
38
int max = col; // permutation row
39
double largest = Double.NEGATIVE_INFINITY;
40![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int row = col; row < m; row++)
{
41
final double[] luRow = lu[row];
42
sum = luRow[col];
43![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int i = 0; i < col; i++)
{
44
sum -= luRow[i] * lu[i][col];
45
}
46
luRow[col] = sum;
47![](/Images/OutliningIndicators/InBlock.gif)
48
// maintain best permutation choice
49![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if (Math.abs(sum) > largest)
{
50
largest = Math.abs(sum);
51
max = row;
52
}
53
}
54![](/Images/OutliningIndicators/InBlock.gif)
55
// Singularity check
56![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if (Math.abs(lu[max][col]) < singularityThreshold)
{
57
singular = true;
58
return;
59
}
60![](/Images/OutliningIndicators/InBlock.gif)
61
// Pivot if necessary
62![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if (max != col)
{
63
double tmp = 0;
64
final double[] luMax = lu[max];
65
final double[] luCol = lu[col];
66![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
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![](/Images/OutliningIndicators/InBlock.gif)
77
// Divide the lower elements by the "winning" diagonal elt.
78
final double luDiag = lu[col][col];
79![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
for (int row = col + 1; row < m; row++)
{
80
lu[row][col] /= luDiag;
81
}
82
}
83![](/Images/OutliningIndicators/InBlock.gif)
84
}
85![](/Images/OutliningIndicators/None.gif)
测试代码示例:
1![](/Images/OutliningIndicators/ExpandedBlockStart.gif)
/** *//**
2
*
3
*/
4
package algorithm.math;
5![](/Images/OutliningIndicators/None.gif)
6
import org.apache.commons.math.linear.ArrayRealVector;
7
import org.apache.commons.math.linear.DecompositionSolver;
8
import org.apache.commons.math.linear.LUDecomposition;
9
import org.apache.commons.math.linear.LUDecompositionImpl;
10
import org.apache.commons.math.linear.MatrixUtils;
11
import org.apache.commons.math.linear.RealMatrix;
12![](/Images/OutliningIndicators/None.gif)
13![](/Images/OutliningIndicators/ExpandedBlockStart.gif)
/** *//**
14
* @author Jia Yu
15
* @date 2010-11-20
16
*/
17![](/Images/OutliningIndicators/ExpandedBlockStart.gif)
public class MatrixDecompositionTest
{
18![](/Images/OutliningIndicators/InBlock.gif)
19![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
/** *//**
20
* @param args
21
*/
22![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
public static void main(String[] args)
{
23
// TODO Auto-generated method stub
24
matrixDecomposite();
25
}
26![](/Images/OutliningIndicators/InBlock.gif)
27![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
private static void matrixDecomposite()
{
28
// TODO Auto-generated method stub
29![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
double[][] testData =
{
30![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ 1.0, 2.0, 3.0},
31![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ 2.0, 5.0, 3.0},
32![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ 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![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
RealMatrix b = MatrixUtils.createRealMatrix(new double[][]
{
53![](/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ 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![](/Images/OutliningIndicators/InBlock.gif)
62
}
63![](/Images/OutliningIndicators/None.gif)
输出结果:
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