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
018package org.apache.commons.math3.linear;
019
020import 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 * @since 2.0 (changed to concrete class in 3.0)
043 */
044public class RectangularCholeskyDecomposition {
045
046    /** Permutated Cholesky root of the symmetric positive semidefinite matrix. */
047    private final RealMatrix root;
048
049    /** Rank of the symmetric positive semidefinite matrix. */
050    private int rank;
051
052    /**
053     * Decompose a symmetric positive semidefinite matrix.
054     * <p>
055     * <b>Note:</b> this constructor follows the linpack method to detect dependent
056     * columns by proceeding with the Cholesky algorithm until a nonpositive diagonal
057     * element is encountered.
058     *
059     * @see <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf">
060     * Analysis of the Cholesky Decomposition of a Semi-definite Matrix</a>
061     *
062     * @param matrix Symmetric positive semidefinite matrix.
063     * @exception NonPositiveDefiniteMatrixException if the matrix is not
064     * positive semidefinite.
065     * @since 3.1
066     */
067    public RectangularCholeskyDecomposition(RealMatrix matrix)
068        throws NonPositiveDefiniteMatrixException {
069        this(matrix, 0);
070    }
071
072    /**
073     * Decompose a symmetric positive semidefinite matrix.
074     *
075     * @param matrix Symmetric positive semidefinite matrix.
076     * @param small Diagonal elements threshold under which columns are
077     * considered to be dependent on previous ones and are discarded.
078     * @exception NonPositiveDefiniteMatrixException if the matrix is not
079     * positive semidefinite.
080     */
081    public RectangularCholeskyDecomposition(RealMatrix matrix, double small)
082        throws NonPositiveDefiniteMatrixException {
083
084        final int order = matrix.getRowDimension();
085        final double[][] c = matrix.getData();
086        final double[][] b = new double[order][order];
087
088        int[] index = new int[order];
089        for (int i = 0; i < order; ++i) {
090            index[i] = i;
091        }
092
093        int r = 0;
094        for (boolean loop = true; loop;) {
095
096            // find maximal diagonal element
097            int swapR = r;
098            for (int i = r + 1; i < order; ++i) {
099                int ii  = index[i];
100                int isr = index[swapR];
101                if (c[ii][ii] > c[isr][isr]) {
102                    swapR = i;
103                }
104            }
105
106
107            // swap elements
108            if (swapR != r) {
109                final int tmpIndex    = index[r];
110                index[r]              = index[swapR];
111                index[swapR]          = tmpIndex;
112                final double[] tmpRow = b[r];
113                b[r]                  = b[swapR];
114                b[swapR]              = tmpRow;
115            }
116
117            // check diagonal element
118            int ir = index[r];
119            if (c[ir][ir] <= small) {
120
121                if (r == 0) {
122                    throw new NonPositiveDefiniteMatrixException(c[ir][ir], ir, small);
123                }
124
125                // check remaining diagonal elements
126                for (int i = r; i < order; ++i) {
127                    if (c[index[i]][index[i]] < -small) {
128                        // there is at least one sufficiently negative diagonal element,
129                        // the symmetric positive semidefinite matrix is wrong
130                        throw new NonPositiveDefiniteMatrixException(c[index[i]][index[i]], i, small);
131                    }
132                }
133
134                // all remaining diagonal elements are close to zero, we consider we have
135                // found the rank of the symmetric positive semidefinite matrix
136                loop = false;
137
138            } else {
139
140                // transform the matrix
141                final double sqrt = FastMath.sqrt(c[ir][ir]);
142                b[r][r] = sqrt;
143                final double inverse  = 1 / sqrt;
144                final double inverse2 = 1 / c[ir][ir];
145                for (int i = r + 1; i < order; ++i) {
146                    final int ii = index[i];
147                    final double e = inverse * c[ii][ir];
148                    b[i][r] = e;
149                    c[ii][ii] -= c[ii][ir] * c[ii][ir] * inverse2;
150                    for (int j = r + 1; j < i; ++j) {
151                        final int ij = index[j];
152                        final double f = c[ii][ij] - e * b[j][r];
153                        c[ii][ij] = f;
154                        c[ij][ii] = f;
155                    }
156                }
157
158                // prepare next iteration
159                loop = ++r < order;
160            }
161        }
162
163        // build the root matrix
164        rank = r;
165        root = MatrixUtils.createRealMatrix(order, r);
166        for (int i = 0; i < order; ++i) {
167            for (int j = 0; j < r; ++j) {
168                root.setEntry(index[i], j, b[i][j]);
169            }
170        }
171
172    }
173
174    /** Get the root of the covariance matrix.
175     * The root is the rectangular matrix <code>B</code> such that
176     * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
177     * @return root of the square matrix
178     * @see #getRank()
179     */
180    public RealMatrix getRootMatrix() {
181        return root;
182    }
183
184    /** Get the rank of the symmetric positive semidefinite matrix.
185     * The r is the number of independent rows in the symmetric positive semidefinite
186     * matrix, it is also the number of columns of the rectangular
187     * matrix of the decomposition.
188     * @return r of the square matrix.
189     * @see #getRootMatrix()
190     */
191    public int getRank() {
192        return rank;
193    }
194
195}