CholeskyDecomposition.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math4.legacy.linear;

  20. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  21. import org.apache.commons.math4.core.jdkmath.JdkMath;


  24. /**
  25.  * Calculates the Cholesky decomposition of a matrix.
  26.  * <p>The Cholesky decomposition of a real symmetric positive-definite
  27.  * matrix A consists of a lower triangular matrix L with same size such
  28.  * that: A = LL<sup>T</sup>. In a sense, this is the square root of A.</p>
  29.  * <p>This class is based on the class with similar name from the
  30.  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
  31.  * following changes:</p>
  32.  * <ul>
  33.  *   <li>a {@link #getLT() getLT} method has been added,</li>
  34.  *   <li>the {@code isspd} method has been removed, since the constructor of
  35.  *   this class throws a {@link NonPositiveDefiniteMatrixException} when a
  36.  *   matrix cannot be decomposed,</li>
  37.  *   <li>a {@link #getDeterminant() getDeterminant} method has been added,</li>
  38.  *   <li>the {@code solve} method has been replaced by a {@link #getSolver()
  39.  *   getSolver} method and the equivalent method provided by the returned
  40.  *   {@link DecompositionSolver}.</li>
  41.  * </ul>
  42.  *
  43.  * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
  44.  * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
  45.  * @since 2.0 (changed to concrete class in 3.0)
  46.  */
  47. public class CholeskyDecomposition {
  48.     /**
  49.      * Default threshold above which off-diagonal elements are considered too different
  50.      * and matrix not symmetric.
  51.      */
  52.     public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
  53.     /**
  54.      * Default threshold below which diagonal elements are considered null
  55.      * and matrix not positive definite.
  56.      */
  57.     public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
  58.     /** Row-oriented storage for L<sup>T</sup> matrix data. */
  59.     private final double[][] lTData;
  60.     /** Cached value of L. */
  61.     private RealMatrix cachedL;
  62.     /** Cached value of LT. */
  63.     private RealMatrix cachedLT;

  65.     /**
  66.      * Calculates the Cholesky decomposition of the given matrix.
  67.      * <p>
  68.      * Calling this constructor is equivalent to call {@link
  69.      * #CholeskyDecomposition(RealMatrix, double, double)} with the
  70.      * thresholds set to the default values {@link
  71.      * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
  72.      * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
  73.      * </p>
  74.      * @param matrix the matrix to decompose
  75.      * @throws NonSquareMatrixException if the matrix is not square.
  76.      * @throws NonSymmetricMatrixException if the matrix is not symmetric.
  77.      * @throws NonPositiveDefiniteMatrixException if the matrix is not
  78.      * strictly positive definite.
  79.      * @see #CholeskyDecomposition(RealMatrix, double, double)
  80.      * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
  81.      * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
  82.      */
  83.     public CholeskyDecomposition(final RealMatrix matrix) {
  84.         this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
  85.              DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
  86.     }

  88.     /**
  89.      * Calculates the Cholesky decomposition of the given matrix.
  90.      * @param matrix the matrix to decompose
  91.      * @param relativeSymmetryThreshold threshold above which off-diagonal
  92.      * elements are considered too different and matrix not symmetric
  93.      * @param absolutePositivityThreshold threshold below which diagonal
  94.      * elements are considered null and matrix not positive definite
  95.      * @throws NonSquareMatrixException if the matrix is not square.
  96.      * @throws NonSymmetricMatrixException if the matrix is not symmetric.
  97.      * @throws NonPositiveDefiniteMatrixException if the matrix is not
  98.      * strictly positive definite.
  99.      * @see #CholeskyDecomposition(RealMatrix)
  100.      * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
  101.      * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
  102.      */
  103.     public CholeskyDecomposition(final RealMatrix matrix,
  104.                                      final double relativeSymmetryThreshold,
  105.                                      final double absolutePositivityThreshold) {
  106.         if (!matrix.isSquare()) {
  107.             throw new NonSquareMatrixException(matrix.getRowDimension(),
  108.                                                matrix.getColumnDimension());
  109.         }

  111.         final int order = matrix.getRowDimension();
  112.         lTData   = matrix.getData();
  113.         cachedL  = null;
  114.         cachedLT = null;

  116.         // check the matrix before transformation
  117.         for (int i = 0; i < order; ++i) {
  118.             final double[] lI = lTData[i];

  120.             // check off-diagonal elements (and reset them to 0)
  121.             for (int j = i + 1; j < order; ++j) {
  122.                 final double[] lJ = lTData[j];
  123.                 final double lIJ = lI[j];
  124.                 final double lJI = lJ[i];
  125.                 final double maxDelta =
  126.                     relativeSymmetryThreshold * JdkMath.max(JdkMath.abs(lIJ), JdkMath.abs(lJI));
  127.                 if (JdkMath.abs(lIJ - lJI) > maxDelta) {
  128.                     throw new NonSymmetricMatrixException(i, j, relativeSymmetryThreshold);
  129.                 }
  130.                 lJ[i] = 0;
  131.            }
  132.         }

  134.         // transform the matrix
  135.         for (int i = 0; i < order; ++i) {

  137.             final double[] ltI = lTData[i];

  139.             // check diagonal element
  140.             if (ltI[i] <= absolutePositivityThreshold) {
  141.                 throw new NonPositiveDefiniteMatrixException(ltI[i], i, absolutePositivityThreshold);
  142.             }

  144.             ltI[i] = JdkMath.sqrt(ltI[i]);
  145.             final double inverse = 1.0 / ltI[i];

  147.             for (int q = order - 1; q > i; --q) {
  148.                 ltI[q] *= inverse;
  149.                 final double[] ltQ = lTData[q];
  150.                 for (int p = q; p < order; ++p) {
  151.                     ltQ[p] -= ltI[q] * ltI[p];
  152.                 }
  153.             }
  154.         }
  155.     }

  157.     /**
  158.      * Returns the matrix L of the decomposition.
  159.      * <p>L is an lower-triangular matrix</p>
  160.      * @return the L matrix
  161.      */
  162.     public RealMatrix getL() {
  163.         if (cachedL == null) {
  164.             cachedL = getLT().transpose();
  165.         }
  166.         return cachedL;
  167.     }

  169.     /**
  170.      * Returns the transpose of the matrix L of the decomposition.
  171.      * <p>L<sup>T</sup> is an upper-triangular matrix</p>
  172.      * @return the transpose of the matrix L of the decomposition
  173.      */
  174.     public RealMatrix getLT() {

  176.         if (cachedLT == null) {
  177.             cachedLT = MatrixUtils.createRealMatrix(lTData);
  178.         }

  180.         // return the cached matrix
  181.         return cachedLT;
  182.     }

  184.     /**
  185.      * Return the determinant of the matrix.
  186.      * @return determinant of the matrix
  187.      */
  188.     public double getDeterminant() {
  189.         double determinant = 1.0;
  190.         for (int i = 0; i < lTData.length; ++i) {
  191.             double lTii = lTData[i][i];
  192.             determinant *= lTii * lTii;
  193.         }
  194.         return determinant;
  195.     }

  197.     /**
  198.      * Get a solver for finding the A &times; X = B solution in least square sense.
  199.      * @return a solver
  200.      */
  201.     public DecompositionSolver getSolver() {
  202.         return new Solver(lTData);
  203.     }

  205.     /** Specialized solver. */
  206.     private static final class Solver implements DecompositionSolver {
  207.         /** Row-oriented storage for L<sup>T</sup> matrix data. */
  208.         private final double[][] lTData;

  210.         /**
  211.          * Build a solver from decomposed matrix.
  212.          * @param lTData row-oriented storage for L<sup>T</sup> matrix data
  213.          */
  214.         private Solver(final double[][] lTData) {
  215.             this.lTData = lTData;
  216.         }

  218.         /** {@inheritDoc} */
  219.         @Override
  220.         public boolean isNonSingular() {
  221.             // if we get this far, the matrix was positive definite, hence non-singular
  222.             return true;
  223.         }

  225.         /** {@inheritDoc} */
  226.         @Override
  227.         public RealVector solve(final RealVector b) {
  228.             final int m = lTData.length;
  229.             if (b.getDimension() != m) {
  230.                 throw new DimensionMismatchException(b.getDimension(), m);
  231.             }

  233.             final double[] x = b.toArray();

  235.             // Solve LY = b
  236.             for (int j = 0; j < m; j++) {
  237.                 final double[] lJ = lTData[j];
  238.                 x[j] /= lJ[j];
  239.                 final double xJ = x[j];
  240.                 for (int i = j + 1; i < m; i++) {
  241.                     x[i] -= xJ * lJ[i];
  242.                 }
  243.             }

  245.             // Solve LTX = Y
  246.             for (int j = m - 1; j >= 0; j--) {
  247.                 x[j] /= lTData[j][j];
  248.                 final double xJ = x[j];
  249.                 for (int i = 0; i < j; i++) {
  250.                     x[i] -= xJ * lTData[i][j];
  251.                 }
  252.             }

  254.             return new ArrayRealVector(x, false);
  255.         }

  257.         /** {@inheritDoc} */
  258.         @Override
  259.         public RealMatrix solve(RealMatrix b) {
  260.             final int m = lTData.length;
  261.             if (b.getRowDimension() != m) {
  262.                 throw new DimensionMismatchException(b.getRowDimension(), m);
  263.             }

  265.             final int nColB = b.getColumnDimension();
  266.             final double[][] x = b.getData();

  268.             // Solve LY = b
  269.             for (int j = 0; j < m; j++) {
  270.                 final double[] lJ = lTData[j];
  271.                 final double lJJ = lJ[j];
  272.                 final double[] xJ = x[j];
  273.                 for (int k = 0; k < nColB; ++k) {
  274.                     xJ[k] /= lJJ;
  275.                 }
  276.                 for (int i = j + 1; i < m; i++) {
  277.                     final double[] xI = x[i];
  278.                     final double lJI = lJ[i];
  279.                     for (int k = 0; k < nColB; ++k) {
  280.                         xI[k] -= xJ[k] * lJI;
  281.                     }
  282.                 }
  283.             }

  285.             // Solve LTX = Y
  286.             for (int j = m - 1; j >= 0; j--) {
  287.                 final double lJJ = lTData[j][j];
  288.                 final double[] xJ = x[j];
  289.                 for (int k = 0; k < nColB; ++k) {
  290.                     xJ[k] /= lJJ;
  291.                 }
  292.                 for (int i = 0; i < j; i++) {
  293.                     final double[] xI = x[i];
  294.                     final double lIJ = lTData[i][j];
  295.                     for (int k = 0; k < nColB; ++k) {
  296.                         xI[k] -= xJ[k] * lIJ;
  297.                     }
  298.                 }
  299.             }

  301.             return new Array2DRowRealMatrix(x);
  302.         }

  304.         /**
  305.          * Get the inverse of the decomposed matrix.
  306.          *
  307.          * @return the inverse matrix.
  308.          */
  309.         @Override
  310.         public RealMatrix getInverse() {
  311.             return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
  312.         }
  313.     }
  314. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT