diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java | 356 |
1 files changed, 356 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java new file mode 100644 index 0000000..bde380e --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java @@ -0,0 +1,356 @@ +/* + * 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; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.util.FastMath; + + +/** + * 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: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + * @since 2.0 + */ +public class CholeskyDecompositionImpl implements CholeskyDecomposition { + + /** 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 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 * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI)); + if (FastMath.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]; + + // check diagonal element + if (ltI[i] < absolutePositivityThreshold) { + throw new NotPositiveDefiniteMatrixException(); + } + + ltI[i] = FastMath.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 { + + /** Row-oriented storage for L<sup>T</sup> matrix data. */ + private final double[][] lTData; + + /** + * Build a solver from decomposed matrix. + * @param lTData 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( + LocalizedFormats.VECTOR_LENGTH_MISMATCH, + 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((ArrayRealVector) b); + } catch (ClassCastException cce) { + + final int m = lTData.length; + if (b.getDimension() != m) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.VECTOR_LENGTH_MISMATCH, + 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 ArrayRealVector(x, false); + + } + } + + /** Solve the linear equation A × X = B. + * <p>The A matrix is implicit here. It is </p> + * @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 ArrayRealVector solve(ArrayRealVector b) + throws IllegalArgumentException, InvalidMatrixException { + return new ArrayRealVector(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( + LocalizedFormats.DIMENSIONS_MISMATCH_2x2, + 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 Array2DRowRealMatrix(x, false); + + } + + /** {@inheritDoc} */ + public RealMatrix getInverse() throws InvalidMatrixException { + return solve(MatrixUtils.createRealIdentityMatrix(lTData.length)); + } + + } + +} |