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.math.linear;
019    
020    import org.apache.commons.math.MathRuntimeException;
021    
022    
023    /**
024     * Calculates the Cholesky decomposition of a matrix.
025     * <p>The Cholesky decomposition of a real symmetric positive-definite
026     * matrix A consists of a lower triangular matrix L with same size that
027     * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
028     *
029     * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
030     * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
031     * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
032     * @since 2.0
033     */
034    public class CholeskyDecompositionImpl implements CholeskyDecomposition {
035    
036        /** Default threshold above which off-diagonal elements are considered too different
037         * and matrix not symmetric. */
038        public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
039    
040        /** Default threshold below which diagonal elements are considered null
041         * and matrix not positive definite. */
042        public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
043    
044        /** Row-oriented storage for L<sup>T</sup> matrix data. */
045        private double[][] lTData;
046    
047        /** Cached value of L. */
048        private RealMatrix cachedL;
049    
050        /** Cached value of LT. */
051        private RealMatrix cachedLT;
052    
053        /**
054         * Calculates the Cholesky decomposition of the given matrix.
055         * <p>
056         * Calling this constructor is equivalent to call {@link
057         * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
058         * thresholds set to the default values {@link
059         * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
060         * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
061         * </p>
062         * @param matrix the matrix to decompose
063         * @exception NonSquareMatrixException if matrix is not square
064         * @exception NotSymmetricMatrixException if matrix is not symmetric
065         * @exception NotPositiveDefiniteMatrixException if the matrix is not
066         * strictly positive definite
067         * @see #CholeskyDecompositionImpl(RealMatrix, double, double)
068         * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
069         * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
070         */
071        public CholeskyDecompositionImpl(final RealMatrix matrix)
072            throws NonSquareMatrixException,
073                   NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
074            this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
075                 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
076        }
077    
078        /**
079         * Calculates the Cholesky decomposition of the given matrix.
080         * @param matrix the matrix to decompose
081         * @param relativeSymmetryThreshold threshold above which off-diagonal
082         * elements are considered too different and matrix not symmetric
083         * @param absolutePositivityThreshold threshold below which diagonal
084         * elements are considered null and matrix not positive definite
085         * @exception NonSquareMatrixException if matrix is not square
086         * @exception NotSymmetricMatrixException if matrix is not symmetric
087         * @exception NotPositiveDefiniteMatrixException if the matrix is not
088         * strictly positive definite
089         * @see #CholeskyDecompositionImpl(RealMatrix)
090         * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
091         * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
092         */
093        public CholeskyDecompositionImpl(final RealMatrix matrix,
094                                         final double relativeSymmetryThreshold,
095                                         final double absolutePositivityThreshold)
096            throws NonSquareMatrixException,
097                   NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
098    
099            if (!matrix.isSquare()) {
100                throw new NonSquareMatrixException(matrix.getRowDimension(),
101                                                   matrix.getColumnDimension());
102            }
103    
104            final int order = matrix.getRowDimension();
105            lTData   = matrix.getData();
106            cachedL  = null;
107            cachedLT = null;
108    
109            // check the matrix before transformation
110            for (int i = 0; i < order; ++i) {
111    
112                final double[] lI = lTData[i];
113    
114                // check off-diagonal elements (and reset them to 0)
115                for (int j = i + 1; j < order; ++j) {
116                    final double[] lJ = lTData[j];
117                    final double lIJ = lI[j];
118                    final double lJI = lJ[i];
119                    final double maxDelta =
120                        relativeSymmetryThreshold * Math.max(Math.abs(lIJ), Math.abs(lJI));
121                    if (Math.abs(lIJ - lJI) > maxDelta) {
122                        throw new NotSymmetricMatrixException();
123                    }
124                    lJ[i] = 0;
125               }
126            }
127    
128            // transform the matrix
129            for (int i = 0; i < order; ++i) {
130    
131                final double[] ltI = lTData[i];
132    
133                // check diagonal element
134                if (ltI[i] < absolutePositivityThreshold) {
135                    throw new NotPositiveDefiniteMatrixException();
136                }
137    
138                ltI[i] = Math.sqrt(ltI[i]);
139                final double inverse = 1.0 / ltI[i];
140    
141                for (int q = order - 1; q > i; --q) {
142                    ltI[q] *= inverse;
143                    final double[] ltQ = lTData[q];
144                    for (int p = q; p < order; ++p) {
145                        ltQ[p] -= ltI[q] * ltI[p];
146                    }
147                }
148    
149            }
150    
151        }
152    
153        /** {@inheritDoc} */
154        public RealMatrix getL() {
155            if (cachedL == null) {
156                cachedL = getLT().transpose();
157            }
158            return cachedL;
159        }
160    
161        /** {@inheritDoc} */
162        public RealMatrix getLT() {
163    
164            if (cachedLT == null) {
165                cachedLT = MatrixUtils.createRealMatrix(lTData);
166            }
167    
168            // return the cached matrix
169            return cachedLT;
170    
171        }
172    
173        /** {@inheritDoc} */
174        public double getDeterminant() {
175            double determinant = 1.0;
176            for (int i = 0; i < lTData.length; ++i) {
177                double lTii = lTData[i][i];
178                determinant *= lTii * lTii;
179            }
180            return determinant;
181        }
182    
183        /** {@inheritDoc} */
184        public DecompositionSolver getSolver() {
185            return new Solver(lTData);
186        }
187    
188        /** Specialized solver. */
189        private static class Solver implements DecompositionSolver {
190    
191            /** Row-oriented storage for L<sup>T</sup> matrix data. */
192            private final double[][] lTData;
193    
194            /**
195             * Build a solver from decomposed matrix.
196             * @param lTData row-oriented storage for L<sup>T</sup> matrix data
197             */
198            private Solver(final double[][] lTData) {
199                this.lTData = lTData;
200            }
201    
202            /** {@inheritDoc} */
203            public boolean isNonSingular() {
204                // if we get this far, the matrix was positive definite, hence non-singular
205                return true;
206            }
207    
208            /** {@inheritDoc} */
209            public double[] solve(double[] b)
210                throws IllegalArgumentException, InvalidMatrixException {
211    
212                final int m = lTData.length;
213                if (b.length != m) {
214                    throw MathRuntimeException.createIllegalArgumentException(
215                            "vector length mismatch: got {0} but expected {1}",
216                            b.length, m);
217                }
218    
219                final double[] x = b.clone();
220    
221                // Solve LY = b
222                for (int j = 0; j < m; j++) {
223                    final double[] lJ = lTData[j];
224                    x[j] /= lJ[j];
225                    final double xJ = x[j];
226                    for (int i = j + 1; i < m; i++) {
227                        x[i] -= xJ * lJ[i];
228                    }
229                }
230    
231                // Solve LTX = Y
232                for (int j = m - 1; j >= 0; j--) {
233                    x[j] /= lTData[j][j];
234                    final double xJ = x[j];
235                    for (int i = 0; i < j; i++) {
236                        x[i] -= xJ * lTData[i][j];
237                    }
238                }
239    
240                return x;
241    
242            }
243    
244            /** {@inheritDoc} */
245            public RealVector solve(RealVector b)
246                throws IllegalArgumentException, InvalidMatrixException {
247                try {
248                    return solve((ArrayRealVector) b);
249                } catch (ClassCastException cce) {
250    
251                    final int m = lTData.length;
252                    if (b.getDimension() != m) {
253                        throw MathRuntimeException.createIllegalArgumentException(
254                                "vector length mismatch: got {0} but expected {1}",
255                                b.getDimension(), m);
256                    }
257    
258                    final double[] x = b.getData();
259    
260                    // Solve LY = b
261                    for (int j = 0; j < m; j++) {
262                        final double[] lJ = lTData[j];
263                        x[j] /= lJ[j];
264                        final double xJ = x[j];
265                        for (int i = j + 1; i < m; i++) {
266                            x[i] -= xJ * lJ[i];
267                        }
268                    }
269    
270                    // Solve LTX = Y
271                    for (int j = m - 1; j >= 0; j--) {
272                        x[j] /= lTData[j][j];
273                        final double xJ = x[j];
274                        for (int i = 0; i < j; i++) {
275                            x[i] -= xJ * lTData[i][j];
276                        }
277                    }
278    
279                    return new ArrayRealVector(x, false);
280    
281                }
282            }
283    
284            /** Solve the linear equation A &times; X = B.
285             * <p>The A matrix is implicit here. It is </p>
286             * @param b right-hand side of the equation A &times; X = B
287             * @return a vector X such that A &times; X = B
288             * @exception IllegalArgumentException if matrices dimensions don't match
289             * @exception InvalidMatrixException if decomposed matrix is singular
290             */
291            public ArrayRealVector solve(ArrayRealVector b)
292                throws IllegalArgumentException, InvalidMatrixException {
293                return new ArrayRealVector(solve(b.getDataRef()), false);
294            }
295    
296            /** {@inheritDoc} */
297            public RealMatrix solve(RealMatrix b)
298                throws IllegalArgumentException, InvalidMatrixException {
299    
300                final int m = lTData.length;
301                if (b.getRowDimension() != m) {
302                    throw MathRuntimeException.createIllegalArgumentException(
303                            "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
304                            b.getRowDimension(), b.getColumnDimension(), m, "n");
305                }
306    
307                final int nColB = b.getColumnDimension();
308                double[][] x = b.getData();
309    
310                // Solve LY = b
311                for (int j = 0; j < m; j++) {
312                    final double[] lJ = lTData[j];
313                    final double lJJ = lJ[j];
314                    final double[] xJ = x[j];
315                    for (int k = 0; k < nColB; ++k) {
316                        xJ[k] /= lJJ;
317                    }
318                    for (int i = j + 1; i < m; i++) {
319                        final double[] xI = x[i];
320                        final double lJI = lJ[i];
321                        for (int k = 0; k < nColB; ++k) {
322                            xI[k] -= xJ[k] * lJI;
323                        }
324                    }
325                }
326    
327                // Solve LTX = Y
328                for (int j = m - 1; j >= 0; j--) {
329                    final double lJJ = lTData[j][j];
330                    final double[] xJ = x[j];
331                    for (int k = 0; k < nColB; ++k) {
332                        xJ[k] /= lJJ;
333                    }
334                    for (int i = 0; i < j; i++) {
335                        final double[] xI = x[i];
336                        final double lIJ = lTData[i][j];
337                        for (int k = 0; k < nColB; ++k) {
338                            xI[k] -= xJ[k] * lIJ;
339                        }
340                    }
341                }
342    
343                return new Array2DRowRealMatrix(x, false);
344    
345            }
346    
347            /** {@inheritDoc} */
348            public RealMatrix getInverse() throws InvalidMatrixException {
349                return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
350            }
351    
352        }
353    
354    }