    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.util.FastMath;
    021
    022    /**
    023     * Calculates the rectangular Cholesky decomposition of a matrix.
    024     * <p>The rectangular Cholesky decomposition of a real symmetric positive
    025     * semidefinite matrix A consists of a rectangular matrix B with the same
    026     * number of rows such that: A is almost equal to BB<sup>T</sup>, depending
    027     * on a user-defined tolerance. In a sense, this is the square root of A.</p>
    028     * <p>The difference with respect to the regular {@link CholeskyDecomposition}
    029     * is that rows/columns may be permuted (hence the rectangular shape instead
    030     * of the traditional triangular shape) and there is a threshold to ignore
    031     * small diagonal elements. This is used for example to generate {@link
    032     * org.apache.commons.math3.random.CorrelatedRandomVectorGenerator correlated
    033     * random n-dimensions vectors} in a p-dimension subspace (p < n).
    034     * In other words, it allows generating random vectors from a covariance
    035     * matrix that is only positive semidefinite, and not positive definite.</p>
    036     * <p>Rectangular Cholesky decomposition is <em>not</em> suited for solving
    037     * linear systems, so it does not provide any {@link DecompositionSolver
    038     * decomposition solver}.</p>
    039     *
    040     * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
    041     * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
    042     * @version $Id: RectangularCholeskyDecomposition.java 1422313 2012-12-15 18:53:41Z psteitz $
    043     * @since 2.0 (changed to concrete class in 3.0)
    044     */
    045    public class RectangularCholeskyDecomposition {
    046
    047        /** Permutated Cholesky root of the symmetric positive semidefinite matrix. */
    048        private final RealMatrix root;
    049
    050        /** Rank of the symmetric positive semidefinite matrix. */
    051        private int rank;
    052
    053        /**
    054         * Decompose a symmetric positive semidefinite matrix.
    055         * <p>
    056         * <b>Note:</b> this constructor follows the linpack method to detect dependent
    057         * columns by proceeding with the Cholesky algorithm until a nonpositive diagonal
    058         * element is encountered.
    059         *
    060         * @see <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf">
    061         * Analysis of the Cholesky Decomposition of a Semi-definite Matrix</a>
    062         *
    063         * @param matrix Symmetric positive semidefinite matrix.
    064         * @exception NonPositiveDefiniteMatrixException if the matrix is not
    065         * positive semidefinite.
    066         * @since 3.1
    067         */
    068        public RectangularCholeskyDecomposition(RealMatrix matrix)
    069            throws NonPositiveDefiniteMatrixException {
    070            this(matrix, 0);
    071        }
    072
    073        /**
    074         * Decompose a symmetric positive semidefinite matrix.
    075         *
    076         * @param matrix Symmetric positive semidefinite matrix.
    077         * @param small Diagonal elements threshold under which columns are
    078         * considered to be dependent on previous ones and are discarded.
    079         * @exception NonPositiveDefiniteMatrixException if the matrix is not
    080         * positive semidefinite.
    081         */
    082        public RectangularCholeskyDecomposition(RealMatrix matrix, double small)
    083            throws NonPositiveDefiniteMatrixException {
    084
    085            final int order = matrix.getRowDimension();
    086            final double[][] c = matrix.getData();
    087            final double[][] b = new double[order][order];
    088
    089            int[] index = new int[order];
    090            for (int i = 0; i < order; ++i) {
    091                index[i] = i;
    092            }
    093
    094            int r = 0;
    095            for (boolean loop = true; loop;) {
    096
    097                // find maximal diagonal element
    098                int swapR = r;
    099                for (int i = r + 1; i < order; ++i) {
    100                    int ii  = index[i];
    101                    int isr = index[swapR];
    102                    if (c[ii][ii] > c[isr][isr]) {
    103                        swapR = i;
    104                    }
    105                }
    106
    107
    108                // swap elements
    109                if (swapR != r) {
    110                    final int tmpIndex    = index[r];
    111                    index[r]              = index[swapR];
    112                    index[swapR]          = tmpIndex;
    113                    final double[] tmpRow = b[r];
    114                    b[r]                  = b[swapR];
    115                    b[swapR]              = tmpRow;
    116                }
    117
    118                // check diagonal element
    119                int ir = index[r];
    120                if (c[ir][ir] <= small) {
    121
    122                    if (r == 0) {
    123                        throw new NonPositiveDefiniteMatrixException(c[ir][ir], ir, small);
    124                    }
    125
    126                    // check remaining diagonal elements
    127                    for (int i = r; i < order; ++i) {
    128                        if (c[index[i]][index[i]] < -small) {
    129                            // there is at least one sufficiently negative diagonal element,
    130                            // the symmetric positive semidefinite matrix is wrong
    131                            throw new NonPositiveDefiniteMatrixException(c[index[i]][index[i]], i, small);
    132                        }
    133                    }
    134
    135                    // all remaining diagonal elements are close to zero, we consider we have
    136                    // found the rank of the symmetric positive semidefinite matrix
    137                    loop = false;
    138
    139                } else {
    140
    141                    // transform the matrix
    142                    final double sqrt = FastMath.sqrt(c[ir][ir]);
    143                    b[r][r] = sqrt;
    144                    final double inverse  = 1 / sqrt;
    145                    final double inverse2 = 1 / c[ir][ir];
    146                    for (int i = r + 1; i < order; ++i) {
    147                        final int ii = index[i];
    148                        final double e = inverse * c[ii][ir];
    149                        b[i][r] = e;
    150                        c[ii][ii] -= c[ii][ir] * c[ii][ir] * inverse2;
    151                        for (int j = r + 1; j < i; ++j) {
    152                            final int ij = index[j];
    153                            final double f = c[ii][ij] - e * b[j][r];
    154                            c[ii][ij] = f;
    155                            c[ij][ii] = f;
    156                        }
    157                    }
    158
    159                    // prepare next iteration
    160                    loop = ++r < order;
    161                }
    162            }
    163
    164            // build the root matrix
    165            rank = r;
    166            root = MatrixUtils.createRealMatrix(order, r);
    167            for (int i = 0; i < order; ++i) {
    168                for (int j = 0; j < r; ++j) {
    169                    root.setEntry(index[i], j, b[i][j]);
    170                }
    171            }
    172
    173        }
    174
    175        /** Get the root of the covariance matrix.
    176         * The root is the rectangular matrix <code>B</code> such that
    177         * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
    178         * @return root of the square matrix
    179         * @see #getRank()
    180         */
    181        public RealMatrix getRootMatrix() {
    182            return root;
    183        }
    184
    185        /** Get the rank of the symmetric positive semidefinite matrix.
    186         * The r is the number of independent rows in the symmetric positive semidefinite
    187         * matrix, it is also the number of columns of the rectangular
    188         * matrix of the decomposition.
    189         * @return r of the square matrix.
    190         * @see #getRootMatrix()
    191         */
    192        public int getRank() {
    193            return rank;
    194        }
    195
    196    }