Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added Cholesky decomposition 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@744708 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-02-15 17:53:32 +00:00
parent 23a8fc2592
commit 098b003e18
8 changed files with 767 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/java/org/apache/commons/math/MessagesResources_fr.java
View File

@ -68,10 +68,14 @@ public class MessagesResources_fr
{ "dimension mismatch {0} != {1}", { "dimension mismatch {0} != {1}",
"dimensions incompatibles {0} != {1}" }, "dimensions incompatibles {0} != {1}" },
// org.apache.commons.math.random.NotPositiveDefiniteMatrixException // org.apache.commons.math.linear.NotPositiveDefiniteMatrixException
{ "not positive definite matrix", { "not positive definite matrix",
"matrice non d\u00e9finie positive" }, "matrice non d\u00e9finie positive" },
// org.apache.commons.math.linear.NotSymmetricMatrixException
{ "not symmetric matrix",
"matrice non symm\u00e9trique" },
// org.apache.commons.math.fraction.FractionConversionException // org.apache.commons.math.fraction.FractionConversionException
{ "Unable to convert {0} to fraction after {1} iterations", { "Unable to convert {0} to fraction after {1} iterations",
"Impossible de convertir {0} en fraction apr\u00e8s {1} it\u00e9rations" }, "Impossible de convertir {0} en fraction apr\u00e8s {1} it\u00e9rations" },

72
src/java/org/apache/commons/math/linear/CholeskyDecomposition.java Normal file
View File

@ -0,0 +1,72 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import java.io.Serializable;
/**
* An interface to classes that implement an algorithm to calculate the
* Cholesky decomposition of a real symmetric positive-definite matrix.
* <p>This interface is based on the class with similar name from the now defunct
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
* following changes:</p>
* <ul>
* <li>a {@link #getLT() getLT} method has been added,</li>
* <li>the <code>isspd</code> method has been removed, the constructors of
* implementation classes being expected to throw {@link
* NotPositiveDefiniteMatrixException} when a matrix cannot be decomposed,</li>
* <li>a {@link #getDeterminant() getDeterminant} method has been added,</li>
* <li>the <code>solve</code> method has been replaced by a {@link
* #getSolver() getSolver} method and the equivalent method provided by
* the returned {@link DecompositionSolver}.</li>
* </ul>
*
* @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
* @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
* @version $Revision$ $Date$
* @since 2.0
*/
public interface CholeskyDecomposition extends Serializable {
/**
* Returns the matrix L of the decomposition.
* <p>L is an lower-triangular matrix</p>
* @return the L matrix
*/
RealMatrix getL();
/**
* Returns the transpose of the matrix L of the decomposition.
* <p>L<sup>T</sup> is an upper-triangular matrix</p>
* @return the transpose of the matrix L of the decomposition
*/
RealMatrix getLT();
/**
* Return the determinant of the matrix
* @return determinant of the matrix
*/
double getDeterminant();
/**
* Get a solver for finding the A &times; X = B solution in least square sense.
* @return a solver
*/
DecompositionSolver getSolver();
}

359
src/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java Normal file
View File

@ -0,0 +1,359 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import org.apache.commons.math.MathRuntimeException;
/**
* Calculates the Cholesky decomposition of a matrix.
* <p>The Cholesky decomposition of a real symmetric positive-definite
* matrix A consists of a lower triangular matrix L with same size that
* satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
*
* @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
* @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
* @version $Revision$ $Date$
* @since 2.0
*/
public class CholeskyDecompositionImpl implements CholeskyDecomposition {
/** Serializable version identifier. */
private static final long serialVersionUID = -2036131698031167221L;
/** Default threshold above which off-diagonal elements are considered too different
* and matrix not symmetric. */
public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
/** Default threshold below which diagonal elements are considered null
* and matrix not positive definite. */
public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
/** Row-oriented storage for L<sup>T</sup> matrix data. */
private double[][] lTData;
/** Cached value of L. */
private RealMatrix cachedL;
/** Cached value of LT. */
private RealMatrix cachedLT;
/**
* Calculates the Cholesky decomposition of the given matrix.
* <p>
* Calling this constructor is equivalent to call {@link
* #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
* thresholds set to the default values {@link
* #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
* #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
* </p>
* @param matrix the matrix to decompose
* @exception NonSquareMatrixException if matrix is not square
* @exception NotSymmetricMatrixException if matrix is not symmetric
* @exception NotPositiveDefiniteMatrixException if the matrix is not
* strictly positive definite
* @see #CholeskyDecompositionImpl(RealMatrix, double, double)
* @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
* @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
*/
public CholeskyDecompositionImpl(final RealMatrix matrix)
throws NonSquareMatrixException,
NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
}
/**
* Calculates the Cholesky decomposition of the given matrix.
* @param matrix the matrix to decompose
* @param relativeSymmetryThreshold threshold above which off-diagonal
* elements are considered too different and matrix not symmetric
* @param absolutePositivityThreshold threshold below which diagonal
* elements are considered null and matrix not positive definite
* @exception NonSquareMatrixException if matrix is not square
* @exception NotSymmetricMatrixException if matrix is not symmetric
* @exception NotPositiveDefiniteMatrixException if the matrix is not
* strictly positive definite
* @see #CholeskyDecompositionImpl(RealMatrix)
* @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
* @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
*/
public CholeskyDecompositionImpl(final RealMatrix matrix,
final double relativeSymmetryThreshold,
final double absolutePositivityThreshold)
throws NonSquareMatrixException,
NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
if (!matrix.isSquare()) {
throw new NonSquareMatrixException(matrix.getRowDimension(),
matrix.getColumnDimension());
}
final int order = matrix.getRowDimension();
lTData = matrix.getData();
cachedL = null;
cachedLT = null;
// check the matrix before transformation
for (int i = 0; i < order; ++i) {
final double[] lI = lTData[i];
// check diagonal element
if (lTData[i][i] < absolutePositivityThreshold) {
throw new NotPositiveDefiniteMatrixException();
}
// check off-diagonal elements (and reset them to 0)
for (int j = i + 1; j < order; ++j) {
final double[] lJ = lTData[j];
final double lIJ = lI[j];
final double lJI = lJ[i];
final double maxDelta =
relativeSymmetryThreshold * Math.max(Math.abs(lIJ), Math.abs(lJI));
if (Math.abs(lIJ - lJI) > maxDelta) {
throw new NotSymmetricMatrixException();
}
lJ[i] = 0;
}
}
// transform the matrix
for (int i = 0; i < order; ++i) {
final double[] ltI = lTData[i];
ltI[i] = Math.sqrt(ltI[i]);
final double inverse = 1.0 / ltI[i];
for (int q = order - 1; q > i; --q) {
ltI[q] *= inverse;
final double[] ltQ = lTData[q];
for (int p = q; p < order; ++p) {
ltQ[p] -= ltI[q] * ltI[p];
}
}
}
}
/** {@inheritDoc} */
public RealMatrix getL() {
if (cachedL == null) {
cachedL = getLT().transpose();
}
return cachedL;
}
/** {@inheritDoc} */
public RealMatrix getLT() {
if (cachedLT == null) {
cachedLT = MatrixUtils.createRealMatrix(lTData);
}
// return the cached matrix
return cachedLT;
}
/** {@inheritDoc} */
public double getDeterminant() {
double determinant = 1.0;
for (int i = 0; i < lTData.length; ++i) {
double lTii = lTData[i][i];
determinant *= lTii * lTii;
}
return determinant;
}
/** {@inheritDoc} */
public DecompositionSolver getSolver() {
return new Solver(lTData);
}
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = -7288829864732555901L;
/** Row-oriented storage for L<sup>T</sup> matrix data. */
private final double[][] lTData;
/**
* Build a solver from decomposed matrix.
* @param lData row-oriented storage for L<sup>T</sup> matrix data
*/
private Solver(final double[][] lTData) {
this.lTData = lTData;
}
/** {@inheritDoc} */
public boolean isNonSingular() {
// if we get this far, the matrix was positive definite, hence non-singular
return true;
}
/** {@inheritDoc} */
public double[] solve(double[] b)
throws IllegalArgumentException, InvalidMatrixException {
final int m = lTData.length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
"vector length mismatch: got {0} but expected {1}",
new Object[] { b.length, m });
}
final double[] x = b.clone();
// Solve LY = b
for (int j = 0; j < m; j++) {
final double[] lJ = lTData[j];
x[j] /= lJ[j];
final double xJ = x[j];
for (int i = j + 1; i < m; i++) {
x[i] -= xJ * lJ[i];
}
}
// Solve LTX = Y
for (int j = m - 1; j >= 0; j--) {
x[j] /= lTData[j][j];
final double xJ = x[j];
for (int i = 0; i < j; i++) {
x[i] -= xJ * lTData[i][j];
}
}
return x;
}
/** {@inheritDoc} */
public RealVector solve(RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
try {
return solve((RealVectorImpl) b);
} catch (ClassCastException cce) {
final int m = lTData.length;
if (b.getDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
"vector length mismatch: got {0} but expected {1}",
new Object[] { b.getDimension(), m });
}
final double[] x = b.getData();
// Solve LY = b
for (int j = 0; j < m; j++) {
final double[] lJ = lTData[j];
x[j] /= lJ[j];
final double xJ = x[j];
for (int i = j + 1; i < m; i++) {
x[i] -= xJ * lJ[i];
}
}
// Solve LTX = Y
for (int j = m - 1; j >= 0; j--) {
x[j] /= lTData[j][j];
final double xJ = x[j];
for (int i = 0; i < j; i++) {
x[i] -= xJ * lTData[i][j];
}
}
return new RealVectorImpl(x, false);
}
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It is </p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X such that A &times; X = B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public RealVectorImpl solve(RealVectorImpl b)
throws IllegalArgumentException, InvalidMatrixException {
return new RealVectorImpl(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
final int m = lTData.length;
if (b.getRowDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
"dimensions mismatch: got {0}x{1} but expected {2}x{3}",
new Object[] { b.getRowDimension(), b.getColumnDimension(), m, "n"});
}
final int nColB = b.getColumnDimension();
double[][] x = b.getData();
// Solve LY = b
for (int j = 0; j < m; j++) {
final double[] lJ = lTData[j];
final double lJJ = lJ[j];
final double[] xJ = x[j];
for (int k = 0; k < nColB; ++k) {
xJ[k] /= lJJ;
}
for (int i = j + 1; i < m; i++) {
final double[] xI = x[i];
final double lJI = lJ[i];
for (int k = 0; k < nColB; ++k) {
xI[k] -= xJ[k] * lJI;
}
}
}
// Solve LTX = Y
for (int j = m - 1; j >= 0; j--) {
final double lJJ = lTData[j][j];
final double[] xJ = x[j];
for (int k = 0; k < nColB; ++k) {
xJ[k] /= lJJ;
}
for (int i = 0; i < j; i++) {
final double[] xI = x[i];
final double lIJ = lTData[i][j];
for (int k = 0; k < nColB; ++k) {
xI[k] -= xJ[k] * lIJ;
}
}
}
return new RealMatrixImpl(x, false);
}
/** {@inheritDoc} */
public RealMatrix getInverse() throws InvalidMatrixException {
return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
}
}
}

42
src/java/org/apache/commons/math/linear/NotSymmetricMatrixException.java Normal file
View File

@ -0,0 +1,42 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import org.apache.commons.math.MathException;
/**
* This class represents exceptions thrown when a matrix expected to
* be symmetric is not
*
* @since 2.0
* @version $Revision$ $Date$
*/
public class NotSymmetricMatrixException extends MathException {
/** Serializable version identifier */
private static final long serialVersionUID = -7012803946709786097L;
/** Simple constructor.
* build an exception with a default message.
*/
public NotSymmetricMatrixException() {
super("not symmetric matrix", null);
}
}

3
src/site/xdoc/changes.xml
View File

@ -170,7 +170,8 @@ The <action> type attribute can be add,update,fix,remove.
decompose a matrix as a product of several other matrices with predefined decompose a matrix as a product of several other matrices with predefined
properties and shapes depending on the algorithm. These algorithms allow to properties and shapes depending on the algorithm. These algorithms allow to
solve the equation A * X = B, either for an exact linear solution solve the equation A * X = B, either for an exact linear solution
(LU-decomposition) or an exact or least-squares solution (QR-decomposition). (LU-decomposition, Cholesky decomposition) or an exact or least-squares
solution (QR-decomposition).
</action> </action>
<action dev="luc" type="add" issue="MATH-219" due-to="Andrew Berry"> <action dev="luc" type="add" issue="MATH-219" due-to="Andrew Berry">
Added removeData methods for the SimpleRegression class. This allows Added removeData methods for the SimpleRegression class. This allows

3
src/site/xdoc/userguide/linear.xml
View File

@ -132,7 +132,8 @@ RealVector solution = solver.solve(constants);
for X is such that the residual AX-B has minimal norm. If an exact solution for X is such that the residual AX-B has minimal norm. If an exact solution
exist (i.e. if for some X the residual AX-B is exactly 0), then this exact exist (i.e. if for some X the residual AX-B is exactly 0), then this exact
solution is also the solution in least square sense. Some solvers like solution is also the solution in least square sense. Some solvers like
the one obtained from <code>LUDecomposition</code> can only find the solution the one obtained from <code>LUDecomposition</code> or
<code>CholeskyDecomposition</code>code> can only find the solution
for square matrices and when the solution is an exact linear solution. Other for square matrices and when the solution is an exact linear solution. Other
solvers like the one obtained from <code>QRDecomposition</code> solvers like the one obtained from <code>QRDecomposition</code>
are more versatile and can also find solutions with non-square matrix A or when are more versatile and can also find solutions with non-square matrix A or when

152
src/test/org/apache/commons/math/linear/CholeskyDecompositionImplTest.java Normal file
View File

@ -0,0 +1,152 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import org.apache.commons.math.MathException;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class CholeskyDecompositionImplTest extends TestCase {
private double[][] testData = new double[][] {
{ 1, 2, 4, 7, 11 },
{ 2, 13, 23, 38, 58 },
{ 4, 23, 77, 122, 182 },
{ 7, 38, 122, 294, 430 },
{ 11, 58, 182, 430, 855 }
};
public CholeskyDecompositionImplTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(CholeskyDecompositionImplTest.class);
suite.setName("CholeskyDecompositionImpl Tests");
return suite;
}
/** test dimensions */
public void testDimensions() throws MathException {
CholeskyDecomposition llt =
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData));
assertEquals(testData.length, llt.getL().getRowDimension());
assertEquals(testData.length, llt.getL().getColumnDimension());
assertEquals(testData.length, llt.getLT().getRowDimension());
assertEquals(testData.length, llt.getLT().getColumnDimension());
}
/** test non-square matrix */
public void testNonSquare() {
try {
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
} catch (NonSquareMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test non-symmetric matrix */
public void testNotSymmetricMatrixException() {
try {
double[][] changed = testData.clone();
changed[0][changed[0].length - 1] += 1.0e-5;
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(changed));
} catch (NotSymmetricMatrixException e) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test non positive definite matrix */
public void testNotPositiveDefinite() {
try {
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[][] {
{ 14, 11, 13, 15, 24 },
{ 11, 34, 13, 8, 25 },
{ 13, 13, 14, 15, 21 },
{ 15, 8, 15, 18, 23 },
{ 24, 25, 21, 23, 45 }
}));
} catch (NotPositiveDefiniteMatrixException e) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test A = LLT */
public void testAEqualLLT() throws MathException {
RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
CholeskyDecomposition llt = new CholeskyDecompositionImpl(matrix);
RealMatrix l = llt.getL();
RealMatrix lt = llt.getLT();
double norm = l.multiply(lt).subtract(matrix).getNorm();
assertEquals(0, norm, 1.0e-15);
}
/** test that L is lower triangular */
public void testLLowerTriangular() throws MathException {
RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
RealMatrix l = new CholeskyDecompositionImpl(matrix).getL();
for (int i = 0; i < l.getRowDimension(); i++) {
for (int j = i + 1; j < l.getColumnDimension(); j++) {
assertEquals(0.0, l.getEntry(i, j));
}
}
}
/** test that LT is transpose of L */
public void testLTTransposed() throws MathException {
RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
CholeskyDecomposition llt = new CholeskyDecompositionImpl(matrix);
RealMatrix l = llt.getL();
RealMatrix lt = llt.getLT();
double norm = l.subtract(lt.transpose()).getNorm();
assertEquals(0, norm, 1.0e-15);
}
/** test matrices values */
public void testMatricesValues() throws MathException {
RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
{ 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 }
});
CholeskyDecomposition llt =
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData));
// check values against known references
RealMatrix l = llt.getL();
assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
RealMatrix lt = llt.getLT();
assertEquals(0, lt.subtract(lRef.transpose()).getNorm(), 1.0e-13);
// check the same cached instance is returned the second time
assertTrue(l == llt.getL());
assertTrue(lt == llt.getLT());
}
}

133
src/test/org/apache/commons/math/linear/CholeskySolverTest.java Normal file
View File

@ -0,0 +1,133 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
import org.apache.commons.math.MathException;
public class CholeskySolverTest extends TestCase {
private double[][] testData = new double[][] {
{ 1, 2, 4, 7, 11 },
{ 2, 13, 23, 38, 58 },
{ 4, 23, 77, 122, 182 },
{ 7, 38, 122, 294, 430 },
{ 11, 58, 182, 430, 855 }
};
public CholeskySolverTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(CholeskySolverTest.class);
suite.setName("LUSolver Tests");
return suite;
}
/** test solve dimension errors */
public void testSolveDimensionErrors() throws MathException {
DecompositionSolver solver =
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() throws MathException {
DecompositionSolver solver =
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
{ 78, -13, 1 },
{ 414, -62, -1 },
{ 1312, -202, -37 },
{ 2989, -542, 145 },
{ 5510, -1465, 201 }
});
RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
{ 1, 0, 1 },
{ 0, 1, 1 },
{ 2, 1, -4 },
{ 2, 2, 2 },
{ 5, -3, 0 }
});
// using RealMatrix
assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVectorImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVector with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test determinant */
public void testDeterminant() throws MathException {
assertEquals(7290000.0, getDeterminant(MatrixUtils.createRealMatrix(testData)), 1.0e-15);
}
private double getDeterminant(RealMatrix m) throws MathException {
return new CholeskyDecompositionImpl(m).getDeterminant();
}
}