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 }