mirror of
https://github.com/apache/commons-math.git
synced 2025-02-06 10:09:26 +00:00
Merged CholeskyDecomposition and CholeskyDecompositionImpl (see MATH-662).
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1173481 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
3c316a16d8
commit
60e328c213
@ -19,7 +19,7 @@ package org.apache.commons.math.filter;
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.linear.ArrayRealVector;
|
||||
import org.apache.commons.math.linear.CholeskyDecompositionImpl;
|
||||
import org.apache.commons.math.linear.CholeskyDecomposition;
|
||||
import org.apache.commons.math.linear.DecompositionSolver;
|
||||
import org.apache.commons.math.linear.MatrixDimensionMismatchException;
|
||||
import org.apache.commons.math.linear.MatrixUtils;
|
||||
@ -355,7 +355,7 @@ public class KalmanFilter {
|
||||
// invert S
|
||||
// as the error covariance matrix is a symmetric positive
|
||||
// semi-definite matrix, we can use the cholesky decomposition
|
||||
DecompositionSolver solver = new CholeskyDecompositionImpl(s).getSolver();
|
||||
DecompositionSolver solver = new CholeskyDecomposition(s).getSolver();
|
||||
RealMatrix invertedS = solver.getInverse();
|
||||
|
||||
// Inn = z(k) - H * xHat(k)-
|
||||
|
||||
@ -17,25 +17,27 @@
|
||||
|
||||
package org.apache.commons.math.linear;
|
||||
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate the
|
||||
* Cholesky decomposition of a real symmetric positive-definite matrix.
|
||||
* 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>
|
||||
* <p>This interface is based on the class with similar name from the
|
||||
* <p>This class is based on the class with similar name from the
|
||||
* <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 NonPositiveDefiniteMatrixException}
|
||||
* when a matrix cannot be decomposed,</li>
|
||||
* <li>the {@code isspd} method has been removed, since the constructor of
|
||||
* this class throws a {@link NonPositiveDefiniteMatrixException} 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>
|
||||
* <li>the {@code solve} 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>
|
||||
@ -43,30 +45,263 @@ package org.apache.commons.math.linear;
|
||||
* @version $Id$
|
||||
* @since 2.0
|
||||
*/
|
||||
public interface CholeskyDecomposition {
|
||||
public class 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 CholeskyDecomposition(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 CholeskyDecomposition(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(ltI[i], 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];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix L of the decomposition.
|
||||
* <p>L is an lower-triangular matrix</p>
|
||||
* @return the L matrix
|
||||
*/
|
||||
RealMatrix getL();
|
||||
public RealMatrix getL() {
|
||||
if (cachedL == null) {
|
||||
cachedL = getLT().transpose();
|
||||
}
|
||||
return cachedL;
|
||||
}
|
||||
|
||||
/**
|
||||
* 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();
|
||||
public RealMatrix getLT() {
|
||||
|
||||
if (cachedLT == null) {
|
||||
cachedLT = MatrixUtils.createRealMatrix(lTData);
|
||||
}
|
||||
|
||||
// return the cached matrix
|
||||
return cachedLT;
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the determinant of the matrix
|
||||
* @return determinant of the matrix
|
||||
*/
|
||||
double getDeterminant();
|
||||
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;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get a solver for finding the A × X = B solution in least square sense.
|
||||
* @return a solver
|
||||
*/
|
||||
DecompositionSolver getSolver();
|
||||
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 RealVector solve(final RealVector b) {
|
||||
final int m = lTData.length;
|
||||
if (b.getDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getDimension(), m);
|
||||
}
|
||||
|
||||
final double[] x = b.toArray();
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix solve(RealMatrix b) {
|
||||
final int m = lTData.length;
|
||||
if (b.getRowDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getRowDimension(), m);
|
||||
}
|
||||
|
||||
final int nColB = b.getColumnDimension();
|
||||
final 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);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getInverse() {
|
||||
return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,280 +0,0 @@
|
||||
/*
|
||||
* 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$
|
||||
* @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(ltI[i], 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 RealVector solve(final RealVector b) {
|
||||
final int m = lTData.length;
|
||||
if (b.getDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getDimension(), m);
|
||||
}
|
||||
|
||||
final double[] x = b.toArray();
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix solve(RealMatrix b) {
|
||||
final int m = lTData.length;
|
||||
if (b.getRowDimension() != m) {
|
||||
throw new DimensionMismatchException(b.getRowDimension(), m);
|
||||
}
|
||||
|
||||
final int nColB = b.getColumnDimension();
|
||||
final 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);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getInverse() {
|
||||
return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -20,7 +20,7 @@ package org.apache.commons.math.linear;
|
||||
import org.junit.Test;
|
||||
import org.junit.Assert;
|
||||
|
||||
public class CholeskyDecompositionImplTest {
|
||||
public class CholeskyDecompositionTest {
|
||||
|
||||
private double[][] testData = new double[][] {
|
||||
{ 1, 2, 4, 7, 11 },
|
||||
@ -33,8 +33,8 @@ public class CholeskyDecompositionImplTest {
|
||||
/** test dimensions */
|
||||
@Test
|
||||
public void testDimensions() {
|
||||
CholeskyDecompositionImpl llt =
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData));
|
||||
CholeskyDecomposition llt =
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(testData));
|
||||
Assert.assertEquals(testData.length, llt.getL().getRowDimension());
|
||||
Assert.assertEquals(testData.length, llt.getL().getColumnDimension());
|
||||
Assert.assertEquals(testData.length, llt.getLT().getRowDimension());
|
||||
@ -44,7 +44,7 @@ public class CholeskyDecompositionImplTest {
|
||||
/** test non-square matrix */
|
||||
@Test(expected = NonSquareMatrixException.class)
|
||||
public void testNonSquare() {
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(new double[3][2]));
|
||||
}
|
||||
|
||||
/** test non-symmetric matrix */
|
||||
@ -52,13 +52,13 @@ public class CholeskyDecompositionImplTest {
|
||||
public void testNotSymmetricMatrixException() {
|
||||
double[][] changed = testData.clone();
|
||||
changed[0][changed[0].length - 1] += 1.0e-5;
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(changed));
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(changed));
|
||||
}
|
||||
|
||||
/** test non positive definite matrix */
|
||||
@Test(expected = NonPositiveDefiniteMatrixException.class)
|
||||
public void testNotPositiveDefinite() {
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[][] {
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 14, 11, 13, 15, 24 },
|
||||
{ 11, 34, 13, 8, 25 },
|
||||
{ 13, 13, 14, 15, 21 },
|
||||
@ -69,7 +69,7 @@ public class CholeskyDecompositionImplTest {
|
||||
|
||||
@Test(expected = NonPositiveDefiniteMatrixException.class)
|
||||
public void testMath274() {
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[][] {
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 0.40434286, -0.09376327, 0.30328980, 0.04909388 },
|
||||
{-0.09376327, 0.10400408, 0.07137959, 0.04762857 },
|
||||
{ 0.30328980, 0.07137959, 0.30458776, 0.04882449 },
|
||||
@ -82,7 +82,7 @@ public class CholeskyDecompositionImplTest {
|
||||
@Test
|
||||
public void testAEqualLLT() {
|
||||
RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
|
||||
CholeskyDecompositionImpl llt = new CholeskyDecompositionImpl(matrix);
|
||||
CholeskyDecomposition llt = new CholeskyDecomposition(matrix);
|
||||
RealMatrix l = llt.getL();
|
||||
RealMatrix lt = llt.getLT();
|
||||
double norm = l.multiply(lt).subtract(matrix).getNorm();
|
||||
@ -93,7 +93,7 @@ public class CholeskyDecompositionImplTest {
|
||||
@Test
|
||||
public void testLLowerTriangular() {
|
||||
RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
|
||||
RealMatrix l = new CholeskyDecompositionImpl(matrix).getL();
|
||||
RealMatrix l = new CholeskyDecomposition(matrix).getL();
|
||||
for (int i = 0; i < l.getRowDimension(); i++) {
|
||||
for (int j = i + 1; j < l.getColumnDimension(); j++) {
|
||||
Assert.assertEquals(0.0, l.getEntry(i, j), 0.0);
|
||||
@ -105,7 +105,7 @@ public class CholeskyDecompositionImplTest {
|
||||
@Test
|
||||
public void testLTTransposed() {
|
||||
RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
|
||||
CholeskyDecompositionImpl llt = new CholeskyDecompositionImpl(matrix);
|
||||
CholeskyDecomposition llt = new CholeskyDecomposition(matrix);
|
||||
RealMatrix l = llt.getL();
|
||||
RealMatrix lt = llt.getLT();
|
||||
double norm = l.subtract(lt.transpose()).getNorm();
|
||||
@ -122,8 +122,8 @@ public class CholeskyDecompositionImplTest {
|
||||
{ 7, 8, 9, 10, 0 },
|
||||
{ 11, 12, 13, 14, 15 }
|
||||
});
|
||||
CholeskyDecompositionImpl llt =
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData));
|
||||
CholeskyDecomposition llt =
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(testData));
|
||||
|
||||
// check values against known references
|
||||
RealMatrix l = llt.getL();
|
||||
@ -36,7 +36,7 @@ public class CholeskySolverTest {
|
||||
@Test
|
||||
public void testSolveDimensionErrors() {
|
||||
DecompositionSolver solver =
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(testData)).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
|
||||
try {
|
||||
solver.solve(b);
|
||||
@ -62,7 +62,7 @@ public class CholeskySolverTest {
|
||||
@Test
|
||||
public void testSolve() {
|
||||
DecompositionSolver solver =
|
||||
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
|
||||
new CholeskyDecomposition(MatrixUtils.createRealMatrix(testData)).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 78, -13, 1 },
|
||||
{ 414, -62, -1 },
|
||||
@ -106,7 +106,7 @@ public class CholeskySolverTest {
|
||||
}
|
||||
|
||||
private double getDeterminant(RealMatrix m) {
|
||||
return new CholeskyDecompositionImpl(m).getDeterminant();
|
||||
return new CholeskyDecomposition(m).getDeterminant();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user