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