You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by md...@apache.org on 2003/11/23 20:53:40 UTC
cvs commit: jakarta-commons/math/src/experimental/org/apache/commons/math/linear CholeskySolverTest.java CholeskySolver.java
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