View Javadoc
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    *      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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.linear;
24  
25  import org.hipparchus.exception.LocalizedCoreFormats;
26  import org.hipparchus.exception.MathIllegalArgumentException;
27  import org.hipparchus.util.FastMath;
28  
29  
30  /**
31   * Calculates the Cholesky decomposition of a matrix.
32   * <p>The Cholesky decomposition of a real symmetric positive-definite
33   * matrix A consists of a lower triangular matrix L with same size such
34   * that: A = LL<sup>T</sup>. In a sense, this is the square root of A.</p>
35   * <p>This class is based on the class with similar name from the
36   * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
37   * following changes:</p>
38   * <ul>
39   *   <li>a {@link #getLT() getLT} method has been added,</li>
40   *   <li>the {@code isspd} method has been removed, since the constructor of
41   *   this class throws a {@link MathIllegalArgumentException} when a
42   *   matrix cannot be decomposed,</li>
43   *   <li>a {@link #getDeterminant() getDeterminant} method has been added,</li>
44   *   <li>the {@code solve} method has been replaced by a {@link #getSolver()
45   *   getSolver} method and the equivalent method provided by the returned
46   *   {@link DecompositionSolver}.</li>
47   * </ul>
48   *
49   * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
50   * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
51   */
52  public class CholeskyDecomposition {
53      /**
54       * Default threshold above which off-diagonal elements are considered too different
55       * and matrix not symmetric.
56       */
57      public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
58      /**
59       * Default threshold below which diagonal elements are considered null
60       * and matrix not positive definite.
61       */
62      public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
63      /** Row-oriented storage for L<sup>T</sup> matrix data. */
64      private final double[][] lTData;
65      /** Cached value of L. */
66      private RealMatrix cachedL;
67      /** Cached value of LT. */
68      private RealMatrix cachedLT;
69  
70      /**
71       * Calculates the Cholesky decomposition of the given matrix.
72       * <p>
73       * Calling this constructor is equivalent to call {@link
74       * #CholeskyDecomposition(RealMatrix, double, double)} with the
75       * thresholds set to the default values {@link
76       * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
77       * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
78       * </p>
79       * @param matrix the matrix to decompose
80       * @throws MathIllegalArgumentException if the matrix is not square.
81       * @throws MathIllegalArgumentException if the matrix is not symmetric.
82       * @throws MathIllegalArgumentException if the matrix is not
83       * strictly positive definite.
84       * @see #CholeskyDecomposition(RealMatrix, double, double)
85       * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
86       * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
87       */
88      public CholeskyDecomposition(final RealMatrix matrix) {
89          this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
90               DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
91      }
92  
93      /**
94       * Calculates the Cholesky decomposition of the given matrix.
95       * @param matrix the matrix to decompose
96       * @param relativeSymmetryThreshold threshold above which off-diagonal
97       * elements are considered too different and matrix not symmetric
98       * @param absolutePositivityThreshold threshold below which diagonal
99       * elements are considered null and matrix not positive definite
100      * @throws MathIllegalArgumentException if the matrix is not square.
101      * @throws MathIllegalArgumentException if the matrix is not symmetric.
102      * @throws MathIllegalArgumentException if the matrix is not
103      * strictly positive definite.
104      * @see #CholeskyDecomposition(RealMatrix)
105      * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
106      * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
107      */
108     public CholeskyDecomposition(final RealMatrix matrix,
109                                  final double relativeSymmetryThreshold,
110                                  final double absolutePositivityThreshold) {
111         if (!matrix.isSquare()) {
112             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
113                                                    matrix.getRowDimension(), matrix.getColumnDimension());
114         }
115 
116         final int order = matrix.getRowDimension();
117         lTData   = matrix.getData();
118         cachedL  = null;
119         cachedLT = null;
120 
121         // check the matrix before transformation
122         for (int i = 0; i < order; ++i) {
123             final double[] lI = lTData[i];
124 
125             // check off-diagonal elements (and reset them to 0)
126             for (int j = i + 1; j < order; ++j) {
127                 final double[] lJ = lTData[j];
128                 final double lIJ = lI[j];
129                 final double lJI = lJ[i];
130                 final double maxDelta =
131                     relativeSymmetryThreshold * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI));
132                 if (FastMath.abs(lIJ - lJI) > maxDelta) {
133                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SYMMETRIC_MATRIX,
134                                                            i, j, relativeSymmetryThreshold);
135                 }
136                 lJ[i] = 0;
137            }
138         }
139 
140         // transform the matrix
141         for (int i = 0; i < order; ++i) {
142 
143             final double[] ltI = lTData[i];
144 
145             // check diagonal element
146             if (ltI[i] <= absolutePositivityThreshold) {
147                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_POSITIVE_DEFINITE_MATRIX);
148             }
149 
150             ltI[i] = FastMath.sqrt(ltI[i]);
151             final double inverse = 1.0 / ltI[i];
152 
153             for (int q = order - 1; q > i; --q) {
154                 ltI[q] *= inverse;
155                 final double[] ltQ = lTData[q];
156                 for (int p = q; p < order; ++p) {
157                     ltQ[p] -= ltI[q] * ltI[p];
158                 }
159             }
160         }
161     }
162 
163     /**
164      * Returns the matrix L of the decomposition.
165      * <p>L is an lower-triangular matrix</p>
166      * @return the L matrix
167      */
168     public RealMatrix getL() {
169         if (cachedL == null) {
170             cachedL = getLT().transpose();
171         }
172         return cachedL;
173     }
174 
175     /**
176      * Returns the transpose of the matrix L of the decomposition.
177      * <p>L<sup>T</sup> is an upper-triangular matrix</p>
178      * @return the transpose of the matrix L of the decomposition
179      */
180     public RealMatrix getLT() {
181 
182         if (cachedLT == null) {
183             cachedLT = MatrixUtils.createRealMatrix(lTData);
184         }
185 
186         // return the cached matrix
187         return cachedLT;
188     }
189 
190     /**
191      * Return the determinant of the matrix
192      * @return determinant of the matrix
193      */
194     public double getDeterminant() {
195         double determinant = 1.0;
196         for (int i = 0; i < lTData.length; ++i) {
197             double lTii = lTData[i][i];
198             determinant *= lTii * lTii;
199         }
200         return determinant;
201     }
202 
203     /**
204      * Get a solver for finding the A &times; X = B solution in least square sense.
205      * @return a solver
206      */
207     public DecompositionSolver getSolver() {
208         return new Solver();
209     }
210 
211     /** Specialized solver. */
212     private class Solver implements DecompositionSolver {
213 
214         /** {@inheritDoc} */
215         @Override
216         public boolean isNonSingular() {
217             // if we get this far, the matrix was positive definite, hence non-singular
218             return true;
219         }
220 
221         /** {@inheritDoc} */
222         @Override
223         public RealVector solve(final RealVector b) {
224             final int m = lTData.length;
225             if (b.getDimension() != m) {
226                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
227                                                        b.getDimension(), m);
228             }
229 
230             final double[] x = b.toArray();
231 
232             // Solve LY = b
233             for (int j = 0; j < m; j++) {
234                 final double[] lJ = lTData[j];
235                 x[j] /= lJ[j];
236                 final double xJ = x[j];
237                 for (int i = j + 1; i < m; i++) {
238                     x[i] -= xJ * lJ[i];
239                 }
240             }
241 
242             // Solve LTX = Y
243             for (int j = m - 1; j >= 0; j--) {
244                 x[j] /= lTData[j][j];
245                 final double xJ = x[j];
246                 for (int i = 0; i < j; i++) {
247                     x[i] -= xJ * lTData[i][j];
248                 }
249             }
250 
251             return new ArrayRealVector(x, false);
252         }
253 
254         /** {@inheritDoc} */
255         @Override
256         public RealMatrix solve(RealMatrix b) {
257             final int m = lTData.length;
258             if (b.getRowDimension() != m) {
259                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
260                                                        b.getRowDimension(), m);
261             }
262 
263             final int nColB = b.getColumnDimension();
264             final double[][] x = b.getData();
265 
266             // Solve LY = b
267             for (int j = 0; j < m; j++) {
268                 final double[] lJ = lTData[j];
269                 final double lJJ = lJ[j];
270                 final double[] xJ = x[j];
271                 for (int k = 0; k < nColB; ++k) {
272                     xJ[k] /= lJJ;
273                 }
274                 for (int i = j + 1; i < m; i++) {
275                     final double[] xI = x[i];
276                     final double lJI = lJ[i];
277                     for (int k = 0; k < nColB; ++k) {
278                         xI[k] -= xJ[k] * lJI;
279                     }
280                 }
281             }
282 
283             // Solve LTX = Y
284             for (int j = m - 1; j >= 0; j--) {
285                 final double lJJ = lTData[j][j];
286                 final double[] xJ = x[j];
287                 for (int k = 0; k < nColB; ++k) {
288                     xJ[k] /= lJJ;
289                 }
290                 for (int i = 0; i < j; i++) {
291                     final double[] xI = x[i];
292                     final double lIJ = lTData[i][j];
293                     for (int k = 0; k < nColB; ++k) {
294                         xI[k] -= xJ[k] * lIJ;
295                     }
296                 }
297             }
298 
299             return new Array2DRowRealMatrix(x);
300         }
301 
302         /**
303          * Get the inverse of the decomposed matrix.
304          *
305          * @return the inverse matrix.
306          */
307         @Override
308         public RealMatrix getInverse() {
309             return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
310         }
311 
312         /** {@inheritDoc} */
313         @Override
314         public int getRowDimension() {
315             return lTData.length;
316         }
317 
318         /** {@inheritDoc} */
319         @Override
320         public int getColumnDimension() {
321             return lTData[0].length;
322         }
323 
324     }
325 
326 }