1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math4.legacy.linear; 19 20 import org.apache.commons.math4.core.jdkmath.JdkMath; 21 22 /** 23 * Calculates the rectangular Cholesky decomposition of a matrix. 24 * <p>The rectangular Cholesky decomposition of a real symmetric positive 25 * semidefinite matrix A consists of a rectangular matrix B with the same 26 * number of rows such that: A is almost equal to BB<sup>T</sup>, depending 27 * on a user-defined tolerance. In a sense, this is the square root of A.</p> 28 * <p>The difference with respect to the regular {@link CholeskyDecomposition} 29 * is that rows/columns may be permuted (hence the rectangular shape instead 30 * of the traditional triangular shape) and there is a threshold to ignore 31 * small diagonal elements. This is used for example to generate {@link 32 * org.apache.commons.math4.legacy.random.CorrelatedVectorFactory correlated 33 * random n-dimensions vectors} in a p-dimension subspace (p < n). 34 * In other words, it allows generating random vectors from a covariance 35 * matrix that is only positive semidefinite, and not positive definite.</p> 36 * <p>Rectangular Cholesky decomposition is <em>not</em> suited for solving 37 * linear systems, so it does not provide any {@link DecompositionSolver 38 * decomposition solver}.</p> 39 * 40 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a> 41 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a> 42 * @since 2.0 (changed to concrete class in 3.0) 43 */ 44 public class RectangularCholeskyDecomposition { 45 46 /** Permutated Cholesky root of the symmetric positive semidefinite matrix. */ 47 private final RealMatrix root; 48 49 /** Rank of the symmetric positive semidefinite matrix. */ 50 private int rank; 51 52 /** 53 * Decompose a symmetric positive semidefinite matrix. 54 * <p> 55 * <b>Note:</b> this constructor follows the linpack method to detect dependent 56 * columns by proceeding with the Cholesky algorithm until a nonpositive diagonal 57 * element is encountered. 58 * 59 * @see <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf"> 60 * Analysis of the Cholesky Decomposition of a Semi-definite Matrix</a> 61 * 62 * @param matrix Symmetric positive semidefinite matrix. 63 * @exception NonPositiveDefiniteMatrixException if the matrix is not 64 * positive semidefinite. 65 * @since 3.1 66 */ 67 public RectangularCholeskyDecomposition(RealMatrix matrix) 68 throws NonPositiveDefiniteMatrixException { 69 this(matrix, 0); 70 } 71 72 /** 73 * Decompose a symmetric positive semidefinite matrix. 74 * 75 * @param matrix Symmetric positive semidefinite matrix. 76 * @param small Diagonal elements threshold under which columns are 77 * considered to be dependent on previous ones and are discarded. 78 * @exception NonPositiveDefiniteMatrixException if the matrix is not 79 * positive semidefinite. 80 */ 81 public RectangularCholeskyDecomposition(RealMatrix matrix, double small) 82 throws NonPositiveDefiniteMatrixException { 83 84 final int order = matrix.getRowDimension(); 85 final double[][] c = matrix.getData(); 86 final double[][] b = new double[order][order]; 87 88 int[] index = new int[order]; 89 for (int i = 0; i < order; ++i) { 90 index[i] = i; 91 } 92 93 int r = 0; 94 for (boolean loop = true; loop;) { 95 96 // find maximal diagonal element 97 int swapR = r; 98 for (int i = r + 1; i < order; ++i) { 99 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 } else { 138 139 // transform the matrix 140 final double sqrt = JdkMath.sqrt(c[ir][ir]); 141 b[r][r] = sqrt; 142 final double inverse = 1 / sqrt; 143 final double inverse2 = 1 / c[ir][ir]; 144 for (int i = r + 1; i < order; ++i) { 145 final int ii = index[i]; 146 final double e = inverse * c[ii][ir]; 147 b[i][r] = e; 148 c[ii][ii] -= c[ii][ir] * c[ii][ir] * inverse2; 149 for (int j = r + 1; j < i; ++j) { 150 final int ij = index[j]; 151 final double f = c[ii][ij] - e * b[j][r]; 152 c[ii][ij] = f; 153 c[ij][ii] = f; 154 } 155 } 156 157 // prepare next iteration 158 loop = ++r < order; 159 } 160 } 161 162 // build the root matrix 163 rank = r; 164 root = MatrixUtils.createRealMatrix(order, r); 165 for (int i = 0; i < order; ++i) { 166 for (int j = 0; j < r; ++j) { 167 root.setEntry(index[i], j, b[i][j]); 168 } 169 } 170 } 171 172 /** Get the root of the covariance matrix. 173 * The root is the rectangular matrix <code>B</code> such that 174 * the covariance matrix is equal to <code>B.B<sup>T</sup></code> 175 * @return root of the square matrix 176 * @see #getRank() 177 */ 178 public RealMatrix getRootMatrix() { 179 return root; 180 } 181 182 /** Get the rank of the symmetric positive semidefinite matrix. 183 * The r is the number of independent rows in the symmetric positive semidefinite 184 * matrix, it is also the number of columns of the rectangular 185 * matrix of the decomposition. 186 * @return r of the square matrix. 187 * @see #getRootMatrix() 188 */ 189 public int getRank() { 190 return rank; 191 } 192 }