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 }