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 *      http://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
18package org.apache.commons.math.linear;
19
20import org.apache.commons.math.MathRuntimeException;
21import org.apache.commons.math.exception.util.LocalizedFormats;
22import org.apache.commons.math.util.FastMath;
23
24
25/**
26 * Calculates the Cholesky decomposition of a matrix.
27 * <p>The Cholesky decomposition of a real symmetric positive-definite
28 * matrix A consists of a lower triangular matrix L with same size that
29 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
30 *
31 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
32 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
34 * @since 2.0
35 */
36public class CholeskyDecompositionImpl implements CholeskyDecomposition {
37
38    /** Default threshold above which off-diagonal elements are considered too different
39     * and matrix not symmetric. */
40    public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
41
42    /** Default threshold below which diagonal elements are considered null
43     * and matrix not positive definite. */
44    public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
45
46    /** Row-oriented storage for L<sup>T</sup> matrix data. */
47    private double[][] lTData;
48
49    /** Cached value of L. */
50    private RealMatrix cachedL;
51
52    /** Cached value of LT. */
53    private RealMatrix cachedLT;
54
55    /**
56     * Calculates the Cholesky decomposition of the given matrix.
57     * <p>
58     * Calling this constructor is equivalent to call {@link
59     * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
60     * thresholds set to the default values {@link
61     * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
62     * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
63     * </p>
64     * @param matrix the matrix to decompose
65     * @exception NonSquareMatrixException if matrix is not square
66     * @exception NotSymmetricMatrixException if matrix is not symmetric
67     * @exception NotPositiveDefiniteMatrixException if the matrix is not
68     * strictly positive definite
69     * @see #CholeskyDecompositionImpl(RealMatrix, double, double)
70     * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
71     * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
72     */
73    public CholeskyDecompositionImpl(final RealMatrix matrix)
74        throws NonSquareMatrixException,
75               NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
76        this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
77             DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
78    }
79
80    /**
81     * Calculates the Cholesky decomposition of the given matrix.
82     * @param matrix the matrix to decompose
83     * @param relativeSymmetryThreshold threshold above which off-diagonal
84     * elements are considered too different and matrix not symmetric
85     * @param absolutePositivityThreshold threshold below which diagonal
86     * elements are considered null and matrix not positive definite
87     * @exception NonSquareMatrixException if matrix is not square
88     * @exception NotSymmetricMatrixException if matrix is not symmetric
89     * @exception NotPositiveDefiniteMatrixException if the matrix is not
90     * strictly positive definite
91     * @see #CholeskyDecompositionImpl(RealMatrix)
92     * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
93     * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
94     */
95    public CholeskyDecompositionImpl(final RealMatrix matrix,
96                                     final double relativeSymmetryThreshold,
97                                     final double absolutePositivityThreshold)
98        throws NonSquareMatrixException,
99               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}
357