SemiDefinitePositiveCholeskyDecomposition.java

  1. /*
  2.  * Licensed to the Hipparchus project 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 Hipparchus project 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.  *      https://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.hipparchus.linear;

  20. import java.util.Arrays;

  22. import org.hipparchus.exception.LocalizedCoreFormats;
  23. import org.hipparchus.exception.MathIllegalArgumentException;
  24. import org.hipparchus.util.FastMath;

  26. /**
  27.  * Calculates the Cholesky decomposition of a positive semidefinite matrix.
  28.  * <p>
  29.  * The classic Cholesky decomposition ({@link CholeskyDecomposition}) applies to real
  30.  * symmetric positive-definite matrix. This class extends the Cholesky decomposition to
  31.  * positive semidefinite matrix. The main application is for estimation based on the
  32.  * Unscented Kalman Filter.
  33.  *
  34.  * @see "J. Hartikainen, A. Solin, and S. Särkkä. Optimal filtering with Kalman filters and smoothers,
  35.  *       Dept. of Biomedica Engineering and Computational Sciences, Aalto University School of Science, Aug. 2011."
  36.  *
  37.  * @since 2.2
  38.  */
  39. public class SemiDefinitePositiveCholeskyDecomposition {

  41.     /** Default threshold below which elements are not considered positive. */
  42.     public static final double POSITIVITY_THRESHOLD = 1.0e-15;
  43.     /** Cached value of L. */
  44.     private RealMatrix cachedL;
  45.     /** Cached value of LT. */
  46.     private RealMatrix cachedLT;

  48.     /**
  49.      * Calculates the Cholesky decomposition of the given matrix.
  50.      * @param matrix the matrix to decompose
  51.      * @throws MathIllegalArgumentException if the matrix is not square.
  52.      * @see #SemiDefinitePositiveCholeskyDecomposition(RealMatrix, double)
  53.      * @see #POSITIVITY_THRESHOLD
  54.      */
  55.     public SemiDefinitePositiveCholeskyDecomposition(final RealMatrix matrix) {
  56.         this(matrix, POSITIVITY_THRESHOLD);
  57.     }

  59.     /**
  60.      * Calculates the Cholesky decomposition of the given matrix.
  61.      * @param matrix the matrix to decompose
  62.      * @param positivityThreshold threshold below which elements are not considered positive
  63.      * @throws MathIllegalArgumentException if the matrix is not square.
  64.      */
  65.     public SemiDefinitePositiveCholeskyDecomposition(final RealMatrix matrix,
  66.                                                      final double positivityThreshold) {
  67.         if (!matrix.isSquare()) {
  68.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  69.                                                    matrix.getRowDimension(), matrix.getColumnDimension());
  70.         }

  72.         final int order = matrix.getRowDimension();
  73.         final double[][] lTData = matrix.getData();
  74.         cachedL  = MatrixUtils.createRealMatrix(lTData);
  75.         int def  = 1;

  77.         final double[] zeroArray = new double[order];
  78.         Arrays.fill(zeroArray, 0.);

  80.         for (int i = 0; i < order; ++i) {
  81.             cachedL.setColumn(i, zeroArray);
  82.         }

  84.         for (int i = 0; i < order; ++i) {
  85.             for (int j = 0; j < i + 1; j++) {
  86.                 double s = lTData[i][j];
  87.                 for (int k = 0; k < j; ++k) {
  88.                     s = s - cachedL.getEntry(i, k) * cachedL.getEntry(j, k);
  89.                 }
  90.                 if (j < i) {
  91.                     if (cachedL.getEntry(j, j) > FastMath.ulp(1.0)) {
  92.                         cachedL.setEntry(i, j, s / cachedL.getEntry(j, j));
  93.                     } else {
  94.                         cachedL.setEntry(i, j, 0.);
  95.                     }
  96.                 } else {
  97.                     if (s < -positivityThreshold) {
  98.                         s = 0;
  99.                         def = -1;
  100.                     } else if (s < positivityThreshold) {
  101.                         s = 0;
  102.                         def = FastMath.min(0, def);
  103.                     }
  104.                     cachedL.setEntry(j, j, FastMath.sqrt(s));
  105.                 }

  107.             }
  108.         }

  110.         cachedLT = cachedL.transpose();

  112.         if (def < 0) {
  113.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NEGATIVE_DEFINITE_MATRIX);
  114.         }

  116.     }

  118.     /**
  119.      * Returns the matrix L of the decomposition.
  120.      * <p>L is an lower-triangular matrix</p>
  121.      * @return the L matrix
  122.      */
  123.     public RealMatrix getL() {
  124.         // return the cached matrix
  125.         return cachedL;
  126.     }

  128.     /**
  129.      * Returns the transpose of the matrix L of the decomposition.
  130.      * <p>L<sup>T</sup> is an upper-triangular matrix</p>
  131.      * @return the transpose of the matrix L of the decomposition
  132.      */
  133.     public RealMatrix getLT() {
  134.          // return the cached matrix
  135.         return cachedLT;
  136.     }

  138. }
Created with JaCoCo 0.8.12.202403310830