New additions of CholeskySolver contributed by Stefan Koeberle
git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk/src/experimental@141044 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Mark R. Diggory
parent
c9250274d1
commit
82397d8788
|
|
@ -0,0 +1,257 @@
|
||||||
|
/* ====================================================================
|
||||||
|
* 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
|
||||||
|
|
@ -0,0 +1,329 @@
|
||||||
|
/* ====================================================================
|
||||||
|
* 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
|
||||||
Loading…
Reference in New Issue