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