    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    }