From 098b003e187a972710f332e04a9ed5e50de2d857 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sun, 15 Feb 2009 17:53:32 +0000 Subject: [PATCH] added Cholesky decomposition git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@744708 13f79535-47bb-0310-9956-ffa450edef68 --- .../commons/math/MessagesResources_fr.java | 6 +- .../math/linear/CholeskyDecomposition.java | 72 ++++ .../linear/CholeskyDecompositionImpl.java | 359 ++++++++++++++++++ .../linear/NotSymmetricMatrixException.java | 42 ++ src/site/xdoc/changes.xml | 3 +- src/site/xdoc/userguide/linear.xml | 3 +- .../linear/CholeskyDecompositionImplTest.java | 152 ++++++++ .../math/linear/CholeskySolverTest.java | 133 +++++++ 8 files changed, 767 insertions(+), 3 deletions(-) create mode 100644 src/java/org/apache/commons/math/linear/CholeskyDecomposition.java create mode 100644 src/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java create mode 100644 src/java/org/apache/commons/math/linear/NotSymmetricMatrixException.java create mode 100644 src/test/org/apache/commons/math/linear/CholeskyDecompositionImplTest.java create mode 100644 src/test/org/apache/commons/math/linear/CholeskySolverTest.java diff --git a/src/java/org/apache/commons/math/MessagesResources_fr.java b/src/java/org/apache/commons/math/MessagesResources_fr.java index 3ce5d70d5..b78d109df 100644 --- a/src/java/org/apache/commons/math/MessagesResources_fr.java +++ b/src/java/org/apache/commons/math/MessagesResources_fr.java @@ -68,10 +68,14 @@ public class MessagesResources_fr { "dimension mismatch {0} != {1}", "dimensions incompatibles {0} != {1}" }, - // org.apache.commons.math.random.NotPositiveDefiniteMatrixException + // org.apache.commons.math.linear.NotPositiveDefiniteMatrixException { "not positive definite matrix", "matrice non d\u00e9finie positive" }, + // org.apache.commons.math.linear.NotSymmetricMatrixException + { "not symmetric matrix", + "matrice non symm\u00e9trique" }, + // org.apache.commons.math.fraction.FractionConversionException { "Unable to convert {0} to fraction after {1} iterations", "Impossible de convertir {0} en fraction apr\u00e8s {1} it\u00e9rations" }, diff --git a/src/java/org/apache/commons/math/linear/CholeskyDecomposition.java b/src/java/org/apache/commons/math/linear/CholeskyDecomposition.java new file mode 100644 index 000000000..d788517ba --- /dev/null +++ b/src/java/org/apache/commons/math/linear/CholeskyDecomposition.java @@ -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. + *

This interface is based on the class with similar name from the now defunct + * JAMA library, with the + * following changes:

+ * + * + * @see MathWorld + * @see Wikipedia + * @version $Revision$ $Date$ + * @since 2.0 + */ +public interface CholeskyDecomposition extends Serializable { + + /** + * Returns the matrix L of the decomposition. + *

L is an lower-triangular matrix

+ * @return the L matrix + */ + RealMatrix getL(); + + /** + * Returns the transpose of the matrix L of the decomposition. + *

LT is an upper-triangular matrix

+ * @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 × X = B solution in least square sense. + * @return a solver + */ + DecompositionSolver getSolver(); + +} diff --git a/src/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java b/src/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java new file mode 100644 index 000000000..164f9b7c0 --- /dev/null +++ b/src/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java @@ -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. + *

The Cholesky decomposition of a real symmetric positive-definite + * matrix A consists of a lower triangular matrix L with same size that + * satisfy: A = LLTQ = I). In a sense, this is the square root of A.

+ * + * @see MathWorld + * @see Wikipedia + * @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 LT 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. + *

+ * 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} + *

+ * @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 LT matrix data. */ + private final double[][] lTData; + + /** + * Build a solver from decomposed matrix. + * @param lData row-oriented storage for LT 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 × X = B. + *

The A matrix is implicit here. It is

+ * @param b right-hand side of the equation A × X = B + * @return a vector X such that A × 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)); + } + + } + +} diff --git a/src/java/org/apache/commons/math/linear/NotSymmetricMatrixException.java b/src/java/org/apache/commons/math/linear/NotSymmetricMatrixException.java new file mode 100644 index 000000000..5dbe1f5af --- /dev/null +++ b/src/java/org/apache/commons/math/linear/NotSymmetricMatrixException.java @@ -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); + } + +} diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 6802795c5..f67f5d9eb 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -170,7 +170,8 @@ The type attribute can be add,update,fix,remove. decompose a matrix as a product of several other matrices with predefined properties and shapes depending on the algorithm. These algorithms allow to 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). Added removeData methods for the SimpleRegression class. This allows diff --git a/src/site/xdoc/userguide/linear.xml b/src/site/xdoc/userguide/linear.xml index 441cd6744..ee95902b0 100644 --- a/src/site/xdoc/userguide/linear.xml +++ b/src/site/xdoc/userguide/linear.xml @@ -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 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 - the one obtained from LUDecomposition can only find the solution + the one obtained from LUDecomposition or + CholeskyDecompositioncode> can only find the solution for square matrices and when the solution is an exact linear solution. Other solvers like the one obtained from QRDecomposition are more versatile and can also find solutions with non-square matrix A or when diff --git a/src/test/org/apache/commons/math/linear/CholeskyDecompositionImplTest.java b/src/test/org/apache/commons/math/linear/CholeskyDecompositionImplTest.java new file mode 100644 index 000000000..f051e41ce --- /dev/null +++ b/src/test/org/apache/commons/math/linear/CholeskyDecompositionImplTest.java @@ -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()); + + } + +} diff --git a/src/test/org/apache/commons/math/linear/CholeskySolverTest.java b/src/test/org/apache/commons/math/linear/CholeskySolverTest.java new file mode 100644 index 000000000..dc4be6015 --- /dev/null +++ b/src/test/org/apache/commons/math/linear/CholeskySolverTest.java @@ -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(); + } + +}