001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math3.linear;
019    
020    import org.apache.commons.math3.exception.DimensionMismatchException;
021    import org.apache.commons.math3.util.FastMath;
022    
023    
024    /**
025     * Calculates the Cholesky decomposition of a matrix.
026     * <p>The Cholesky decomposition of a real symmetric positive-definite
027     * matrix A consists of a lower triangular matrix L with same size such
028     * that: A = LL<sup>T</sup>. In a sense, this is the square root of A.</p>
029     * <p>This class is based on the class with similar name from the
030     * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
031     * following changes:</p>
032     * <ul>
033     *   <li>a {@link #getLT() getLT} method has been added,</li>
034     *   <li>the {@code isspd} method has been removed, since the constructor of
035     *   this class throws a {@link NonPositiveDefiniteMatrixException} when a
036     *   matrix cannot be decomposed,</li>
037     *   <li>a {@link #getDeterminant() getDeterminant} method has been added,</li>
038     *   <li>the {@code solve} method has been replaced by a {@link #getSolver()
039     *   getSolver} method and the equivalent method provided by the returned
040     *   {@link DecompositionSolver}.</li>
041     * </ul>
042     *
043     * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
044     * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
045     * @version $Id: CholeskyDecomposition.java 1416643 2012-12-03 19:37:14Z tn $
046     * @since 2.0 (changed to concrete class in 3.0)
047     */
048    public class CholeskyDecomposition {
049        /**
050         * Default threshold above which off-diagonal elements are considered too different
051         * and matrix not symmetric.
052         */
053        public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
054        /**
055         * Default threshold below which diagonal elements are considered null
056         * and matrix not positive definite.
057         */
058        public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
059        /** Row-oriented storage for L<sup>T</sup> matrix data. */
060        private double[][] lTData;
061        /** Cached value of L. */
062        private RealMatrix cachedL;
063        /** Cached value of LT. */
064        private RealMatrix cachedLT;
065    
066        /**
067         * Calculates the Cholesky decomposition of the given matrix.
068         * <p>
069         * Calling this constructor is equivalent to call {@link
070         * #CholeskyDecomposition(RealMatrix, double, double)} with the
071         * thresholds set to the default values {@link
072         * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
073         * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
074         * </p>
075         * @param matrix the matrix to decompose
076         * @throws NonSquareMatrixException if the matrix is not square.
077         * @throws NonSymmetricMatrixException if the matrix is not symmetric.
078         * @throws NonPositiveDefiniteMatrixException if the matrix is not
079         * strictly positive definite.
080         * @see #CholeskyDecomposition(RealMatrix, double, double)
081         * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
082         * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
083         */
084        public CholeskyDecomposition(final RealMatrix matrix) {
085            this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
086                 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
087        }
088    
089        /**
090         * Calculates the Cholesky decomposition of the given matrix.
091         * @param matrix the matrix to decompose
092         * @param relativeSymmetryThreshold threshold above which off-diagonal
093         * elements are considered too different and matrix not symmetric
094         * @param absolutePositivityThreshold threshold below which diagonal
095         * elements are considered null and matrix not positive definite
096         * @throws NonSquareMatrixException if the matrix is not square.
097         * @throws NonSymmetricMatrixException if the matrix is not symmetric.
098         * @throws NonPositiveDefiniteMatrixException if the matrix is not
099         * strictly positive definite.
100         * @see #CholeskyDecomposition(RealMatrix)
101         * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
102         * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
103         */
104        public CholeskyDecomposition(final RealMatrix matrix,
105                                         final double relativeSymmetryThreshold,
106                                         final double absolutePositivityThreshold) {
107            if (!matrix.isSquare()) {
108                throw new NonSquareMatrixException(matrix.getRowDimension(),
109                                                   matrix.getColumnDimension());
110            }
111    
112            final int order = matrix.getRowDimension();
113            lTData   = matrix.getData();
114            cachedL  = null;
115            cachedLT = null;
116    
117            // check the matrix before transformation
118            for (int i = 0; i < order; ++i) {
119                final double[] lI = lTData[i];
120    
121                // check off-diagonal elements (and reset them to 0)
122                for (int j = i + 1; j < order; ++j) {
123                    final double[] lJ = lTData[j];
124                    final double lIJ = lI[j];
125                    final double lJI = lJ[i];
126                    final double maxDelta =
127                        relativeSymmetryThreshold * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI));
128                    if (FastMath.abs(lIJ - lJI) > maxDelta) {
129                        throw new NonSymmetricMatrixException(i, j, relativeSymmetryThreshold);
130                    }
131                    lJ[i] = 0;
132               }
133            }
134    
135            // transform the matrix
136            for (int i = 0; i < order; ++i) {
137    
138                final double[] ltI = lTData[i];
139    
140                // check diagonal element
141                if (ltI[i] <= absolutePositivityThreshold) {
142                    throw new NonPositiveDefiniteMatrixException(ltI[i], i, absolutePositivityThreshold);
143                }
144    
145                ltI[i] = FastMath.sqrt(ltI[i]);
146                final double inverse = 1.0 / ltI[i];
147    
148                for (int q = order - 1; q > i; --q) {
149                    ltI[q] *= inverse;
150                    final double[] ltQ = lTData[q];
151                    for (int p = q; p < order; ++p) {
152                        ltQ[p] -= ltI[q] * ltI[p];
153                    }
154                }
155            }
156        }
157    
158        /**
159         * Returns the matrix L of the decomposition.
160         * <p>L is an lower-triangular matrix</p>
161         * @return the L matrix
162         */
163        public RealMatrix getL() {
164            if (cachedL == null) {
165                cachedL = getLT().transpose();
166            }
167            return cachedL;
168        }
169    
170        /**
171         * Returns the transpose of the matrix L of the decomposition.
172         * <p>L<sup>T</sup> is an upper-triangular matrix</p>
173         * @return the transpose of the matrix L of the decomposition
174         */
175        public RealMatrix getLT() {
176    
177            if (cachedLT == null) {
178                cachedLT = MatrixUtils.createRealMatrix(lTData);
179            }
180    
181            // return the cached matrix
182            return cachedLT;
183        }
184    
185        /**
186         * Return the determinant of the matrix
187         * @return determinant of the matrix
188         */
189        public double getDeterminant() {
190            double determinant = 1.0;
191            for (int i = 0; i < lTData.length; ++i) {
192                double lTii = lTData[i][i];
193                determinant *= lTii * lTii;
194            }
195            return determinant;
196        }
197    
198        /**
199         * Get a solver for finding the A &times; X = B solution in least square sense.
200         * @return a solver
201         */
202        public DecompositionSolver getSolver() {
203            return new Solver(lTData);
204        }
205    
206        /** Specialized solver. */
207        private static class Solver implements DecompositionSolver {
208            /** Row-oriented storage for L<sup>T</sup> matrix data. */
209            private final double[][] lTData;
210    
211            /**
212             * Build a solver from decomposed matrix.
213             * @param lTData row-oriented storage for L<sup>T</sup> matrix data
214             */
215            private Solver(final double[][] lTData) {
216                this.lTData = lTData;
217            }
218    
219            /** {@inheritDoc} */
220            public boolean isNonSingular() {
221                // if we get this far, the matrix was positive definite, hence non-singular
222                return true;
223            }
224    
225            /** {@inheritDoc} */
226            public RealVector solve(final RealVector b) {
227                final int m = lTData.length;
228                if (b.getDimension() != m) {
229                    throw new DimensionMismatchException(b.getDimension(), m);
230                }
231    
232                final double[] x = b.toArray();
233    
234                // Solve LY = b
235                for (int j = 0; j < m; j++) {
236                    final double[] lJ = lTData[j];
237                    x[j] /= lJ[j];
238                    final double xJ = x[j];
239                    for (int i = j + 1; i < m; i++) {
240                        x[i] -= xJ * lJ[i];
241                    }
242                }
243    
244                // Solve LTX = Y
245                for (int j = m - 1; j >= 0; j--) {
246                    x[j] /= lTData[j][j];
247                    final double xJ = x[j];
248                    for (int i = 0; i < j; i++) {
249                        x[i] -= xJ * lTData[i][j];
250                    }
251                }
252    
253                return new ArrayRealVector(x, false);
254            }
255    
256            /** {@inheritDoc} */
257            public RealMatrix solve(RealMatrix b) {
258                final int m = lTData.length;
259                if (b.getRowDimension() != m) {
260                    throw new DimensionMismatchException(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            /** {@inheritDoc} */
303            public RealMatrix getInverse() {
304                return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
305            }
306        }
307    }