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 */
17
18 package org.hipparchus.linear;
19
20 import java.util.Arrays;
21
22 import org.hipparchus.exception.LocalizedCoreFormats;
23 import org.hipparchus.exception.MathIllegalArgumentException;
24 import org.hipparchus.util.FastMath;
25
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 {
40
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;
47
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 }
58
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 }
71
72 final int order = matrix.getRowDimension();
73 final double[][] lTData = matrix.getData();
74 cachedL = MatrixUtils.createRealMatrix(lTData);
75 int def = 1;
76
77 final double[] zeroArray = new double[order];
78 Arrays.fill(zeroArray, 0.);
79
80 for (int i = 0; i < order; ++i) {
81 cachedL.setColumn(i, zeroArray);
82 }
83
84 cachedLT = cachedL.transpose();
85
86 for (int i = 0; i < order; ++i) {
87 for (int j = 0; j < i + 1; j++) {
88 double s = lTData[i][j];
89 for (int k = 0; k < j; ++k) {
90 s = s - cachedL.getEntry(i, k) * cachedL.getEntry(j, k);
91 }
92 if (j < i) {
93 if (cachedL.getEntry(j, j) > FastMath.ulp(1.0)) {
94 cachedL.setEntry(i, j, s / cachedL.getEntry(j, j));
95 } else {
96 cachedL.setEntry(i, j, 0.);
97 }
98 } else {
99 if (s < -positivityThreshold) {
100 s = 0;
101 def = -1;
102 } else if (s < positivityThreshold) {
103 s = 0;
104 def = FastMath.min(0, def);
105 }
106 cachedL.setEntry(j, j, FastMath.sqrt(s));
107 }
108
109 }
110 }
111
112 if (def < 0) {
113 throw new MathIllegalArgumentException(LocalizedCoreFormats.NEGATIVE_DEFINITE_MATRIX);
114 }
115
116 }
117
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 }
127
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 cachedLT = getL().transpose();
135 // return the cached matrix
136 return cachedLT;
137 }
138
139 }