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    }