alvinalexander.com | career | drupal | java | mac | mysql | perl | scala | uml | unix

# Commons Math example source code file (CholeskyDecompositionImpl.java)

This example Commons Math source code file (CholeskyDecompositionImpl.java) is included in the DevDaily.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

## Java - Commons Math tags/keywords

choleskydecompositionimpl, default_absolute_positivity_threshold, default_relative_symmetry_threshold, illegalargumentexception, illegalargumentexception, invalidmatrixexception, invalidmatrixexception, nonsquarematrixexception, notpositivedefinitematrixexception, notsymmetricmatrixexception, realmatrix, realmatrix, solver, solver

## The Commons Math CholeskyDecompositionImpl.java source code

```/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements.  See the NOTICE file distributed with
* 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
*
*
* Unless required by applicable law or agreed to in writing, software
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
*/

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>TQ = I). In a sense, this is the square root of A.
*
* @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld
* @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia
* @version \$Revision: 811685 \$ \$Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) \$
* @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 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_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 * 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];

// check diagonal element
if (ltI[i] < absolutePositivityThreshold) {
throw new NotPositiveDefiniteMatrixException();
}

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 {

/** Row-oriented storage for L<sup>T matrix data. */
private final double[][] lTData;

/**
* Build a solver from decomposed matrix.
* @param lTData row-oriented storage for L<sup>T 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}",
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(
"vector length mismatch: got {0} but expected {1}",
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
* @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(
"dimensions mismatch: got {0}x{1} but expected {2}x{3}",
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));
}

}

}
```

## Other Commons Math examples (source code examples)

Here is a short list of links related to this Commons Math CholeskyDecompositionImpl.java source code file:

 ... this post is sponsored by my books ... #1 New Release! FP Best Seller