Change Dir

先知cd——热爱生活是一切艺术的开始

统计

留言簿(18)

积分与排名

“牛”们的博客

各个公司技术

我的链接

淘宝技术

阅读排行榜

评论排行榜

Commons Math学习笔记——矩阵分解

 

看其他篇章到目录选择。

补充上一次的矩阵知识,这次主要讲讲矩阵的一些分解运算——Matrix Decomposition

矩阵分解主要有三种方式:LU分解,QR分解和奇异值分解。当然在Mathlinear包中提供了对应的接口有CholeskyDecompositionEigenDecompositionLUDecompositionQRDecompositionSingularValueDecomposition5种分解方式。

每一个接口定义对应着一个***DecompositionImpl的实现类。

先来看看LU分解。

LU分解是矩阵分解的一种,可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积(有时是它们和一个置换矩阵的乘积)。LU分解主要应用在数值分析中,用来解线性方程、求逆矩阵或计算行列式。

Math库中的LU分解主要是LUP分解,即针对n*n方阵A,找出三个n*n的矩阵LUP,满足PA=LU。其中L是一个单位下三角矩阵,U是一个上三角矩阵,P是一个置换矩阵。非奇异的矩阵(可逆)都有这样一种分解(可证明)。LUP分解的计算方法就是高斯消元。具体的算法见算导第28章或者参看我列出的参考资料。

Commons Math中的实现LUP分解是这样的,主要在LUDecompositionImpl的构造方法中。见源码:

 1public 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.02.03.0},
31                2.05.03.0},
32                1.00.08.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                10 }2-5 }31 }
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


posted on 2010-12-13 09:39 changedi 阅读(4078) 评论(0)  编辑  收藏 所属分类: 数学


只有注册用户登录后才能发表评论。


网站导航: