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.math.linear; 19 20 import org.apache.commons.math.MathRuntimeException; 21 import org.apache.commons.math.exception.util.LocalizedFormats; 22 import org.apache.commons.math.util.FastMath; 23 24 25 /** 26 * Calculates the Cholesky decomposition of a matrix. 27 * <p>The Cholesky decomposition of a real symmetric positive-definite 28 * matrix A consists of a lower triangular matrix L with same size that 29 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p> 30 * 31 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a> 32 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a> 33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $ 34 * @since 2.0 35 */ 36 public class CholeskyDecompositionImpl implements CholeskyDecomposition { 37 38 /** Default threshold above which off-diagonal elements are considered too different 39 * and matrix not symmetric. */ 40 public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15; 41 42 /** Default threshold below which diagonal elements are considered null 43 * and matrix not positive definite. */ 44 public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10; 45 46 /** Row-oriented storage for L<sup>T</sup> matrix data. */ 47 private double[][] lTData; 48 49 /** Cached value of L. */ 50 private RealMatrix cachedL; 51 52 /** Cached value of LT. */ 53 private RealMatrix cachedLT; 54 55 /** 56 * Calculates the Cholesky decomposition of the given matrix. 57 * <p> 58 * Calling this constructor is equivalent to call {@link 59 * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the 60 * thresholds set to the default values {@link 61 * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link 62 * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD} 63 * </p> 64 * @param matrix the matrix to decompose 65 * @exception NonSquareMatrixException if matrix is not square 66 * @exception NotSymmetricMatrixException if matrix is not symmetric 67 * @exception NotPositiveDefiniteMatrixException if the matrix is not 68 * strictly positive definite 69 * @see #CholeskyDecompositionImpl(RealMatrix, double, double) 70 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD 71 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD 72 */ 73 public CholeskyDecompositionImpl(final RealMatrix matrix) 74 throws NonSquareMatrixException, 75 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException { 76 this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD, 77 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD); 78 } 79 80 /** 81 * Calculates the Cholesky decomposition of the given matrix. 82 * @param matrix the matrix to decompose 83 * @param relativeSymmetryThreshold threshold above which off-diagonal 84 * elements are considered too different and matrix not symmetric 85 * @param absolutePositivityThreshold threshold below which diagonal 86 * elements are considered null and matrix not positive definite 87 * @exception NonSquareMatrixException if matrix is not square 88 * @exception NotSymmetricMatrixException if matrix is not symmetric 89 * @exception NotPositiveDefiniteMatrixException if the matrix is not 90 * strictly positive definite 91 * @see #CholeskyDecompositionImpl(RealMatrix) 92 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD 93 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD 94 */ 95 public CholeskyDecompositionImpl(final RealMatrix matrix, 96 final double relativeSymmetryThreshold, 97 final double absolutePositivityThreshold) 98 throws NonSquareMatrixException, 99 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 } 357
