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.math.linear; 019 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.exception.util.LocalizedFormats; 022 import org.apache.commons.math.util.FastMath; 023 024 025 /** 026 * Calculates the Cholesky decomposition of a matrix. 027 * <p>The Cholesky decomposition of a real symmetric positive-definite 028 * matrix A consists of a lower triangular matrix L with same size that 029 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p> 030 * 031 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a> 032 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a> 033 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $ 034 * @since 2.0 035 */ 036 public class CholeskyDecompositionImpl implements CholeskyDecomposition { 037 038 /** Default threshold above which off-diagonal elements are considered too different 039 * and matrix not symmetric. */ 040 public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15; 041 042 /** Default threshold below which diagonal elements are considered null 043 * and matrix not positive definite. */ 044 public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10; 045 046 /** Row-oriented storage for L<sup>T</sup> matrix data. */ 047 private double[][] lTData; 048 049 /** Cached value of L. */ 050 private RealMatrix cachedL; 051 052 /** Cached value of LT. */ 053 private RealMatrix cachedLT; 054 055 /** 056 * Calculates the Cholesky decomposition of the given matrix. 057 * <p> 058 * Calling this constructor is equivalent to call {@link 059 * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the 060 * thresholds set to the default values {@link 061 * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link 062 * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD} 063 * </p> 064 * @param matrix the matrix to decompose 065 * @exception NonSquareMatrixException if matrix is not square 066 * @exception NotSymmetricMatrixException if matrix is not symmetric 067 * @exception NotPositiveDefiniteMatrixException if the matrix is not 068 * strictly positive definite 069 * @see #CholeskyDecompositionImpl(RealMatrix, double, double) 070 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD 071 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD 072 */ 073 public CholeskyDecompositionImpl(final RealMatrix matrix) 074 throws NonSquareMatrixException, 075 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException { 076 this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD, 077 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD); 078 } 079 080 /** 081 * Calculates the Cholesky decomposition of the given matrix. 082 * @param matrix the matrix to decompose 083 * @param relativeSymmetryThreshold threshold above which off-diagonal 084 * elements are considered too different and matrix not symmetric 085 * @param absolutePositivityThreshold threshold below which diagonal 086 * elements are considered null and matrix not positive definite 087 * @exception NonSquareMatrixException if matrix is not square 088 * @exception NotSymmetricMatrixException if matrix is not symmetric 089 * @exception NotPositiveDefiniteMatrixException if the matrix is not 090 * strictly positive definite 091 * @see #CholeskyDecompositionImpl(RealMatrix) 092 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD 093 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD 094 */ 095 public CholeskyDecompositionImpl(final RealMatrix matrix, 096 final double relativeSymmetryThreshold, 097 final double absolutePositivityThreshold) 098 throws NonSquareMatrixException, 099 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException { 100 101 if (!matrix.isSquare()) { 102 throw new NonSquareMatrixException(matrix.getRowDimension(), 103 matrix.getColumnDimension()); 104 } 105 106 final int order = matrix.getRowDimension(); 107 lTData = matrix.getData(); 108 cachedL = null; 109 cachedLT = null; 110 111 // check the matrix before transformation 112 for (int i = 0; i < order; ++i) { 113 114 final double[] lI = lTData[i]; 115 116 // check off-diagonal elements (and reset them to 0) 117 for (int j = i + 1; j < order; ++j) { 118 final double[] lJ = lTData[j]; 119 final double lIJ = lI[j]; 120 final double lJI = lJ[i]; 121 final double maxDelta = 122 relativeSymmetryThreshold * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI)); 123 if (FastMath.abs(lIJ - lJI) > maxDelta) { 124 throw new NotSymmetricMatrixException(); 125 } 126 lJ[i] = 0; 127 } 128 } 129 130 // transform the matrix 131 for (int i = 0; i < order; ++i) { 132 133 final double[] ltI = lTData[i]; 134 135 // check diagonal element 136 if (ltI[i] < absolutePositivityThreshold) { 137 throw new NotPositiveDefiniteMatrixException(); 138 } 139 140 ltI[i] = FastMath.sqrt(ltI[i]); 141 final double inverse = 1.0 / ltI[i]; 142 143 for (int q = order - 1; q > i; --q) { 144 ltI[q] *= inverse; 145 final double[] ltQ = lTData[q]; 146 for (int p = q; p < order; ++p) { 147 ltQ[p] -= ltI[q] * ltI[p]; 148 } 149 } 150 151 } 152 153 } 154 155 /** {@inheritDoc} */ 156 public RealMatrix getL() { 157 if (cachedL == null) { 158 cachedL = getLT().transpose(); 159 } 160 return cachedL; 161 } 162 163 /** {@inheritDoc} */ 164 public RealMatrix getLT() { 165 166 if (cachedLT == null) { 167 cachedLT = MatrixUtils.createRealMatrix(lTData); 168 } 169 170 // return the cached matrix 171 return cachedLT; 172 173 } 174 175 /** {@inheritDoc} */ 176 public double getDeterminant() { 177 double determinant = 1.0; 178 for (int i = 0; i < lTData.length; ++i) { 179 double lTii = lTData[i][i]; 180 determinant *= lTii * lTii; 181 } 182 return determinant; 183 } 184 185 /** {@inheritDoc} */ 186 public DecompositionSolver getSolver() { 187 return new Solver(lTData); 188 } 189 190 /** Specialized solver. */ 191 private static class Solver implements DecompositionSolver { 192 193 /** Row-oriented storage for L<sup>T</sup> matrix data. */ 194 private final double[][] lTData; 195 196 /** 197 * Build a solver from decomposed matrix. 198 * @param lTData row-oriented storage for L<sup>T</sup> matrix data 199 */ 200 private Solver(final double[][] lTData) { 201 this.lTData = lTData; 202 } 203 204 /** {@inheritDoc} */ 205 public boolean isNonSingular() { 206 // if we get this far, the matrix was positive definite, hence non-singular 207 return true; 208 } 209 210 /** {@inheritDoc} */ 211 public double[] solve(double[] b) 212 throws IllegalArgumentException, InvalidMatrixException { 213 214 final int m = lTData.length; 215 if (b.length != m) { 216 throw MathRuntimeException.createIllegalArgumentException( 217 LocalizedFormats.VECTOR_LENGTH_MISMATCH, 218 b.length, m); 219 } 220 221 final double[] x = b.clone(); 222 223 // Solve LY = b 224 for (int j = 0; j < m; j++) { 225 final double[] lJ = lTData[j]; 226 x[j] /= lJ[j]; 227 final double xJ = x[j]; 228 for (int i = j + 1; i < m; i++) { 229 x[i] -= xJ * lJ[i]; 230 } 231 } 232 233 // Solve LTX = Y 234 for (int j = m - 1; j >= 0; j--) { 235 x[j] /= lTData[j][j]; 236 final double xJ = x[j]; 237 for (int i = 0; i < j; i++) { 238 x[i] -= xJ * lTData[i][j]; 239 } 240 } 241 242 return x; 243 244 } 245 246 /** {@inheritDoc} */ 247 public RealVector solve(RealVector b) 248 throws IllegalArgumentException, InvalidMatrixException { 249 try { 250 return solve((ArrayRealVector) b); 251 } catch (ClassCastException cce) { 252 253 final int m = lTData.length; 254 if (b.getDimension() != m) { 255 throw MathRuntimeException.createIllegalArgumentException( 256 LocalizedFormats.VECTOR_LENGTH_MISMATCH, 257 b.getDimension(), m); 258 } 259 260 final double[] x = b.getData(); 261 262 // Solve LY = b 263 for (int j = 0; j < m; j++) { 264 final double[] lJ = lTData[j]; 265 x[j] /= lJ[j]; 266 final double xJ = x[j]; 267 for (int i = j + 1; i < m; i++) { 268 x[i] -= xJ * lJ[i]; 269 } 270 } 271 272 // Solve LTX = Y 273 for (int j = m - 1; j >= 0; j--) { 274 x[j] /= lTData[j][j]; 275 final double xJ = x[j]; 276 for (int i = 0; i < j; i++) { 277 x[i] -= xJ * lTData[i][j]; 278 } 279 } 280 281 return new ArrayRealVector(x, false); 282 283 } 284 } 285 286 /** Solve the linear equation A × X = B. 287 * <p>The A matrix is implicit here. It is </p> 288 * @param b right-hand side of the equation A × X = B 289 * @return a vector X such that A × X = B 290 * @exception IllegalArgumentException if matrices dimensions don't match 291 * @exception InvalidMatrixException if decomposed matrix is singular 292 */ 293 public ArrayRealVector solve(ArrayRealVector b) 294 throws IllegalArgumentException, InvalidMatrixException { 295 return new ArrayRealVector(solve(b.getDataRef()), false); 296 } 297 298 /** {@inheritDoc} */ 299 public RealMatrix solve(RealMatrix b) 300 throws IllegalArgumentException, InvalidMatrixException { 301 302 final int m = lTData.length; 303 if (b.getRowDimension() != m) { 304 throw MathRuntimeException.createIllegalArgumentException( 305 LocalizedFormats.DIMENSIONS_MISMATCH_2x2, 306 b.getRowDimension(), b.getColumnDimension(), m, "n"); 307 } 308 309 final int nColB = b.getColumnDimension(); 310 double[][] x = b.getData(); 311 312 // Solve LY = b 313 for (int j = 0; j < m; j++) { 314 final double[] lJ = lTData[j]; 315 final double lJJ = lJ[j]; 316 final double[] xJ = x[j]; 317 for (int k = 0; k < nColB; ++k) { 318 xJ[k] /= lJJ; 319 } 320 for (int i = j + 1; i < m; i++) { 321 final double[] xI = x[i]; 322 final double lJI = lJ[i]; 323 for (int k = 0; k < nColB; ++k) { 324 xI[k] -= xJ[k] * lJI; 325 } 326 } 327 } 328 329 // Solve LTX = Y 330 for (int j = m - 1; j >= 0; j--) { 331 final double lJJ = lTData[j][j]; 332 final double[] xJ = x[j]; 333 for (int k = 0; k < nColB; ++k) { 334 xJ[k] /= lJJ; 335 } 336 for (int i = 0; i < j; i++) { 337 final double[] xI = x[i]; 338 final double lIJ = lTData[i][j]; 339 for (int k = 0; k < nColB; ++k) { 340 xI[k] -= xJ[k] * lIJ; 341 } 342 } 343 } 344 345 return new Array2DRowRealMatrix(x, false); 346 347 } 348 349 /** {@inheritDoc} */ 350 public RealMatrix getInverse() throws InvalidMatrixException { 351 return solve(MatrixUtils.createRealIdentityMatrix(lTData.length)); 352 } 353 354 } 355 356 }