/*
* 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.exception.DimensionMismatchException;
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 such
* that: A = LL<sup>T</sup>. 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 $Id: CholeskyDecompositionImpl.java 1131229 2011-06-03 20:49:25Z luc $
* @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
* @throws NonSquareMatrixException if the matrix is not square.
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
* @throws NonPositiveDefiniteMatrixException 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) {
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
* @throws NonSquareMatrixException if the matrix is not square.
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
* @throws NonPositiveDefiniteMatrixException 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) {
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 NonSymmetricMatrixException(i, j, relativeSymmetryThreshold);
}
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 NonPositiveDefiniteMatrixException(i, absolutePositivityThreshold);
}
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) {
final int m = lTData.length;
if (b.length != m) {
throw new DimensionMismatchException(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) {
try {
return solve((ArrayRealVector) b);
} catch (ClassCastException cce) {
final int m = lTData.length;
if (b.getDimension() != m) {
throw new DimensionMismatchException(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
* @throws DimensionMismatchException if the matrices dimensions do not match.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public ArrayRealVector solve(ArrayRealVector b) {
return new ArrayRealVector(solve(b.getDataRef()), false);
}
/** Solve the linear equation A × X = B for matrices A.
* <p>The A matrix is implicit, it is provided by the underlying
* decomposition algorithm.</p>
* @param b right-hand side of the equation A × X = B
* @param reuseB if true, the b array will be reused and returned,
* instead of being copied
* @return a matrix X that minimizes the two norm of A × X - B
* @throws org.apache.commons.math.exception.DimensionMismatchException
* if the matrices dimensions do not match.
* @throws SingularMatrixException
* if the decomposed matrix is singular.
*/
private double[][] solve(double[][] b, boolean reuseB) {
final int m = lTData.length;
if (b.length != m) {
throw new DimensionMismatchException(b.length, m);
}
final int nColB = b[0].length;
final double[][] x;
if (reuseB) {
x = b;
} else {
x = new double[b.length][nColB];
for (int i = 0; i < b.length; ++i) {
System.arraycopy(b[i], 0, x[i], 0, nColB);
}
}
// 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 x;
}
/** {@inheritDoc} */
public double[][] solve(double[][] b) {
return solve(b, false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b) {
return new Array2DRowRealMatrix(solve(b.getData(), true), false);
}
/** {@inheritDoc} */
public RealMatrix getInverse() {
return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
}
}
}