commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From mdigg...@apache.org
Subject cvs commit: jakarta-commons/math/src/experimental/org/apache/commons/math/linear CholeskySolverTest.java CholeskySolver.java
Date Sun, 23 Nov 2003 19:53:40 GMT
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 apache@apache.org.
   *
   * 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 {@link 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 apache@apache.org.
   *
   * 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: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message