mdiggory 2003/11/23 11:53:40 Added: math/src/experimental/org/apache/commons/math/linear CholeskySolverTest.java CholeskySolver.java Log: New additions of CholeskySolver contributed by Stefan Koeberle Revision Changes Path 1.1 jakarta-commons/math/src/experimental/org/apache/commons/math/linear/CholeskySolverTest.java Index: CholeskySolverTest.java =================================================================== /* ==================================================================== * The Apache Software License, Version 1.1 * * Copyright (c) 2003 The Apache Software Foundation. All rights * reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in * the documentation and/or other materials provided with the * distribution. * * 3. The end-user documentation included with the redistribution, if * any, must include the following acknowledgement: * "This product includes software developed by the * Apache Software Foundation (http://www.apache.org/)." * Alternately, this acknowledgement may appear in the software itself, * if and wherever such third-party acknowledgements normally appear. * * 4. The names "The Jakarta Project", "Commons", and "Apache Software * Foundation" must not be used to endorse or promote products derived * from this software without prior written permission. For written * permission, please contact [EMAIL PROTECTED] * * 5. Products derived from this software may not be called "Apache" * nor may "Apache" appear in their name without prior written * permission of the Apache Software Foundation. * * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * DISCLAIMED. IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. * ==================================================================== * * This software consists of voluntary contributions made by many * individuals on behalf of the Apache Software Foundation. For more * information on the Apache Software Foundation, please see * <http://www.apache.org/>. */ package org.apache.commons.math.linear; import junit.framework.Test; import junit.framework.TestCase; import junit.framework.TestSuite; import junit.textui.TestRunner; /** * Test cases for the [EMAIL PROTECTED] CholeskySolver} class. * <p> * @author Stefan Koeberle, 11/2003 */ public class CholeskySolverTest extends TestCase { private double[][] m1 = {{1}}; private double m1Det = 1.0d; private double[][] m2 = {{1, 0} , {0, 2}}; private double m2Det = 2.0d; private double[][] m3 = {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}}; private double m3Det = 6.0d; private double[][] m4 = {{1, 0, 0}, {2, 3, 0}, {4, 5, 6}}; private double m4Det = 18.0d; private double[][] m5 = {{ 1, 0, 0, 0, 0}, {-2, 3, 0, 0, 0}, { 4, -5, 6, 0, 0}, { 7, 8, -9, 10, 0}, {11, 12, 13, 14, 15}}; private double m5Det = 2700.0d; private double[][] m6 = {{1, 0, 0}, {2, 0, 0}, {4, 5, 6}}; private double[][] m7 = {{1, 2, 3}, {4, 5, 6}}; /** * Creates a new instance of CholeskySolverTest */ public CholeskySolverTest(String nameOfTest) { super(nameOfTest); }//constructor CholeskySolverTest public void setUp() throws java.lang.Exception { super.setUp(); }//setUp public void tearDown() throws java.lang.Exception { super.tearDown(); }//tearDown public static Test suite() { TestSuite suite = new TestSuite(CholeskySolverTest.class); suite.setName("CholeskySolver Tests"); return suite; }//suite /** * tests CholeskySolver.setNumericalZero() */ public void testNumericalZero() { CholeskySolver solver = new CholeskySolver(); double numericalZero = 77.77d; solver.setNumericalZero(numericalZero); assertEquals(solver.getNumericalZero(), numericalZero, 0.0d); try { solver.decompose( new RealMatrixImpl(new double[][]{{numericalZero/2, 0}, {0, numericalZero/2}})); fail("testing numericalZero"); } catch (IllegalArgumentException e) {} }//testNumericalZero /** * tests CholeskySolver.decompose(...) */ public void testDecompose() { //The following decompositions should succeed. testDecompose(m1, "Decomposing matrix m1"); testDecompose(m2, "Decomposing matrix m2"); testDecompose(m3, "Decomposing matrix m3"); testDecompose(m4, "Decomposing matrix m4"); testDecompose(m5, "Decomposing matrix m5"); //The following decompositions will fail. An IllegalArgumentException //should be thrown. try { testDecompose(m6, "Decomposing matrix m6"); fail("Decomposing matrix m6"); } catch (IllegalArgumentException e) {} try { CholeskySolver solver = new CholeskySolver(); solver.decompose(new RealMatrixImpl(m7)); fail("Decomposing matrix m7"); } catch (IllegalArgumentException e) {} }//testDecomposition /** * tests CholeskySolver.solve(...) */ public void testSolve() { //If there's no matrix, there's no linear euqitation to solve ... try { CholeskySolver solver = new CholeskySolver(); solver.solve(new double[] {1,2,3}); fail("solving a liniar equitation with a missing matrix should fail"); } catch (IllegalStateException e) {} //The following operations should succeed. testSolve(m1, "Solving matrix m1"); testSolve(m2, "Solving matrix m2"); testSolve(m3, "Solving matrix m3"); testSolve(m4, "Solving matrix m4"); testSolve(m5, "Solving matrix m5"); //The following operations will fail. An IllegalArgumentException //should be thrown. try { testSolve(m6, "Solving matrix m6"); fail("Solving matrix m6"); } catch (IllegalArgumentException e) {} try { CholeskySolver solver = new CholeskySolver(); solver.solve(new RealMatrixImpl(m3), new double[] {1, 2, 3, 4}); fail("Solving matrix m3[3x3], v[4]"); } catch (IllegalArgumentException e) {} }//testDecomposition /** * tests CholeskySolver.getDeterminant(...) */ public void testGetDeterminant() { //Since no matrix was decomposed, there's no determinant. try { CholeskySolver solver = new CholeskySolver(); solver.getDeterminant(); fail("Calculating determinant of missing matrix should fail"); } catch (IllegalStateException e) {} //These test will suceed. testGetDeterminant(m1, m1Det, "Calculating determinant of m1"); testGetDeterminant(m2, m2Det, "Calculating determinant of m2"); testGetDeterminant(m3, m3Det, "Calculating determinant of m3"); testGetDeterminant(m4, m4Det, "Calculating determinant of m4"); testGetDeterminant(m5, m5Det, "Calculating determinant of m5"); }//test /** * Generates the matrix * <code>m = lowerTriangularMatrix * lowerTriangularMatrix^T</code>. * If alle diagonalelements of <code>lowerTriangularMatrix</code> are * positiv, <code>m</code> will be positiv definit. * Decomposing <code>m</code> should result in * <code>lowerTriangularMatrix</code> again. So there's a simple test ... */ private void testDecompose(double[][] lowerTriangularMatrix, String message) throws IllegalArgumentException { RealMatrix triangularMatrix = new RealMatrixImpl(lowerTriangularMatrix); RealMatrix pdMatrix = triangularMatrix.multiply(triangularMatrix.transpose()); CholeskySolver solver = new CholeskySolver(); solver.decompose(pdMatrix); assertTrue(message, areEqual(triangularMatrix, solver.getDecomposition(), 1.0E-10)); }//testDecompose /** * Similar to <code> private testDecompose(...)</code>. */ private void testSolve(double[][] lowerTriangularMatrix, String message) { RealMatrix triangularMatrix = new RealMatrixImpl(lowerTriangularMatrix); RealMatrixImpl pdMatrix = (RealMatrixImpl) triangularMatrix.multiply(triangularMatrix.transpose()); CholeskySolver solver = new CholeskySolver(); double[] c = new double[lowerTriangularMatrix.length]; for (int i=0; i<c.length; i++) for (int j=0; j<lowerTriangularMatrix[0].length; j++) c[i] += lowerTriangularMatrix[i][j]; solver.decompose(pdMatrix); RealMatrix x = new RealMatrixImpl(solver.solve(c)); assertTrue(message, areEqual(pdMatrix.multiply(x), new RealMatrixImpl(c), 1.0E-10)); }//testSolve /** * Similar to <code> private testDecompose(...)</code>. */ private void testGetDeterminant(double[][] lowerTriangularMatrix, double determinant, String message) throws IllegalArgumentException { RealMatrix triangularMatrix = new RealMatrixImpl(lowerTriangularMatrix); RealMatrix pdMatrix = triangularMatrix.multiply(triangularMatrix.transpose()); double pdDeterminant = determinant * determinant; CholeskySolver solver = new CholeskySolver(); solver.decompose(pdMatrix); assertEquals(message, solver.getDeterminant(), pdDeterminant, 1.0E-10); }//testGetDeterminant /** * Are <code>m1</code> and <code>m2</code> equal? */ private static boolean areEqual(RealMatrix m1, RealMatrix m2, double delta) { double[][] mv1 = m1.getData(); double[][] mv2 = m2.getData(); if (mv1.length != mv1.length || mv1[0].length != mv2[0].length) return false; for (int i=0; i<mv1.length; i++) for (int j=0; j<mv1[0].length; j++) if (Math.abs(mv1[i][j] -mv2[i][j]) > delta) return false; return true; }//isEqual /** * Executes all tests of this class */ public static void main(String[] args) { System.out.println("Start"); TestRunner runner = new TestRunner(); runner.doRun(CholeskySolverTest.suite()); System.out.println("End"); }//main }//class CholeskySolverTest 1.1 jakarta-commons/math/src/experimental/org/apache/commons/math/linear/CholeskySolver.java Index: CholeskySolver.java =================================================================== /* ==================================================================== * The Apache Software License, Version 1.1 * * Copyright (c) 2003 The Apache Software Foundation. All rights * reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in * the documentation and/or other materials provided with the * distribution. * * 3. The end-user documentation included with the redistribution, if * any, must include the following acknowledgement: * "This product includes software developed by the * Apache Software Foundation (http://www.apache.org/)." * Alternately, this acknowledgement may appear in the software itself, * if and wherever such third-party acknowledgements normally appear. * * 4. The names "The Jakarta Project", "Commons", and "Apache Software * Foundation" must not be used to endorse or promote products derived * from this software without prior written permission. For written * permission, please contact [EMAIL PROTECTED] * * 5. Products derived from this software may not be called "Apache" * nor may "Apache" appear in their name without prior written * permission of the Apache Software Foundation. * * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * DISCLAIMED. IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. * ==================================================================== * * This software consists of voluntary contributions made by many * individuals on behalf of the Apache Software Foundation. For more * information on the Apache Software Foundation, please see * <http://www.apache.org/>. */ package org.apache.commons.math.linear; /** * Solves a linear equitation with symmetrical, positiv definit * coefficient matrix by Cholesky decomposition. * <p> * For every symmetric, positiv definit matrix <code>M</code> there is a * lower triangular matrix <code>L</code> so that <code>L*L^T=M</code>. * <code>L</code> is called the <i>Cholesky decomposition</i> of <code>M</code>. * For any constant vector <code>c</code> it can be used to solve * the linear equitation <code>M*x=L*(L^T*x)=c</code>.<br> * Compared to the LU-decompoistion the Cholesky methods requires only half * the number of operations. * <p> * @author Stefan Koeberle, 11/2003 */ public class CholeskySolver { private double numericalZero = 10E-12; /** The lower triangular matrix */ private RealMatrixImpl decompMatrix; /** * Creates a new instance of CholeskySolver */ public CholeskySolver() { }//constructor CholeskySolver /** * Every double <code>d</code> satisfying * <code>java.lang.Math.abs(d) <= numericalZero</code> * is considered equal to <code>0.0d.</code> */ public void setNumericalZero(double numericalZero) { this.numericalZero = numericalZero; }//setNumericalZero /** * See <code>setNumericalZero</code> */ public double getNumericalZero() { return numericalZero; }//getNumericalZero /** * Calculates the Cholesky-decomposition of the symmetrical, positiv definit * matrix <code>M</code>. * <p> * The decomposition matrix is internally stored. * <p> * @throws IllegalArgumentException if <code>M</code> ist not square or * not positiv definit */ public void decompose(RealMatrix m) throws IllegalArgumentException { decompMatrix = null; double[][] mval = m.getData(); int numRows = m.getRowDimension(); int numCols = m.getColumnDimension(); if (numRows != numCols) throw new IllegalArgumentException("matrix is not square"); double[][] decomp = new double[numRows][numCols]; double sum; //for all columns for (int col=0; col<numCols; col++) { //diagonal element sum = mval[col][col]; for (int k=0; k<col; k++) sum = sum - decomp[col][k]*decomp[col][k]; if (sum <= numericalZero) { throw new IllegalArgumentException( "Matrix is not positiv definit"); } decomp[col][col] += Math.sqrt(sum); //column below diagonal for (int row=col+1; row<numRows; row++) { sum = mval[row][col]; for (int k=0; k<col; k++) sum = sum - decomp[col][k]*decomp[row][k]; decomp[row][col] = sum/decomp[col][col]; }//for }//for all columns decompMatrix = new RealMatrixImpl(decomp); }//decompose /** * Returns the last calculated decomposition matrix. * <p> * Caution: Every call of this Method will return the same object. * Decomposing another matrix will generate a new one. */ public RealMatrixImpl getDecomposition() { return decompMatrix; }//getDecomposition /** * Returns the solution for a linear system with constant vector <code>c</code>. * <p> * This method solves a linear system <code>M*x=c</code> for a symmetrical, * positiv definit coefficient matrix <code>M</code>. Before using this * method the matrix <code>M</code> must have been decomposed. * <p> * @throws IllegalStateException if this methode is called before * a matrix was decomposed * @throws IllegalArgumentException if the dimension of <code>c</code> doesn't * match the row dimension of <code>M</code> */ public double[] solve(double[] c) throws IllegalStateException, IllegalArgumentException { if (decompMatrix == null) { throw new IllegalStateException("no decomposed matrix available"); }//if if (decompMatrix.getColumnDimension() != c.length) throw new IllegalArgumentException("matrix dimension mismatch"); double[][] decomp = decompMatrix.getData(); double[] x = new double[decomp.length]; double sum; //forward elimination for (int i=0; i<x.length; i++) { sum = c[i]; for (int k=0; k<i; k++) sum = sum - decomp[i][k]*x[k]; x[i] = sum / decomp[i][i]; }//forward elimination //backward elimination for (int i=x.length-1; i>=0; i--) { sum = x[i]; for (int k=i+1; k<x.length; k++) sum = sum - decomp[k][i]*x[k]; x[i] = sum / decomp[i][i]; }//backward elimination return x; }//solve /** * Returns the solution for a linear system with a symmetrical, * positiv definit coefficient matrix <code>M</code> and * constant vector <code>c</code>. * <p> * As a side effect, the Cholesky-decomposition <code>L*L^T=M</code> is * calculated and internally stored. * <p> * This is a convenience method for <code><pre> * solver.decompose(m); * solver.solve(c); * </pre></code> * @throws IllegalArgumentException if M ist not square, not positive definit * or the dimensions of <code>M</code> and * <code>c</code> don't match. */ public double[] solve(RealMatrix m, double[] c) throws IllegalArgumentException { decompose(m); return solve(c); }//solve /** * Returns the determinant of the a matrix <code>M</code>. * <p> * Before using this method the matrix <code>M</code> must * have been decomposed. * <p> * @throws IllegalStateException if this method is called before * a matrix was decomposed */ public double getDeterminant() { if (decompMatrix == null) { throw new IllegalStateException("no decomposed matrix available"); }//if double[][] data = decompMatrix.getData(); double res = 1.0d; for (int i=0; i<data.length; i++) { res *= data[i][i]; }//for res = res*res; return res; }//getDeterminant }//class CholeskySolver
--------------------------------------------------------------------- To unsubscribe, e-mail: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]