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.exception.DimensionMismatchException; 021 import org.apache.commons.math3.util.FastMath; 022 023 /** 024 * Calculates the LUP-decomposition of a square matrix. 025 * <p>The LUP-decomposition of a matrix A consists of three matrices L, U and 026 * P that satisfy: P×A = L×U. L is lower triangular (with unit 027 * diagonal terms), U is upper triangular and P is a permutation matrix. All 028 * matrices are m×m.</p> 029 * <p>As shown by the presence of the P matrix, this decomposition is 030 * implemented using partial pivoting.</p> 031 * <p>This class is based on the class with similar name from the 032 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p> 033 * <ul> 034 * <li>a {@link #getP() getP} method has been added,</li> 035 * <li>the {@code det} method has been renamed as {@link #getDeterminant() 036 * getDeterminant},</li> 037 * <li>the {@code getDoublePivot} method has been removed (but the int based 038 * {@link #getPivot() getPivot} method has been kept),</li> 039 * <li>the {@code solve} and {@code isNonSingular} methods have been replaced 040 * by a {@link #getSolver() getSolver} method and the equivalent methods 041 * provided by the returned {@link DecompositionSolver}.</li> 042 * </ul> 043 * 044 * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a> 045 * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a> 046 * @version $Id: LUDecomposition.java 1416643 2012-12-03 19:37:14Z tn $ 047 * @since 2.0 (changed to concrete class in 3.0) 048 */ 049 public class LUDecomposition { 050 /** Default bound to determine effective singularity in LU decomposition. */ 051 private static final double DEFAULT_TOO_SMALL = 1e-11; 052 /** Entries of LU decomposition. */ 053 private final double[][] lu; 054 /** Pivot permutation associated with LU decomposition. */ 055 private final int[] pivot; 056 /** Parity of the permutation associated with the LU decomposition. */ 057 private boolean even; 058 /** Singularity indicator. */ 059 private boolean singular; 060 /** Cached value of L. */ 061 private RealMatrix cachedL; 062 /** Cached value of U. */ 063 private RealMatrix cachedU; 064 /** Cached value of P. */ 065 private RealMatrix cachedP; 066 067 /** 068 * Calculates the LU-decomposition of the given matrix. 069 * This constructor uses 1e-11 as default value for the singularity 070 * threshold. 071 * 072 * @param matrix Matrix to decompose. 073 * @throws NonSquareMatrixException if matrix is not square. 074 */ 075 public LUDecomposition(RealMatrix matrix) { 076 this(matrix, DEFAULT_TOO_SMALL); 077 } 078 079 /** 080 * Calculates the LU-decomposition of the given matrix. 081 * @param matrix The matrix to decompose. 082 * @param singularityThreshold threshold (based on partial row norm) 083 * under which a matrix is considered singular 084 * @throws NonSquareMatrixException if matrix is not square 085 */ 086 public LUDecomposition(RealMatrix matrix, double singularityThreshold) { 087 if (!matrix.isSquare()) { 088 throw new NonSquareMatrixException(matrix.getRowDimension(), 089 matrix.getColumnDimension()); 090 } 091 092 final int m = matrix.getColumnDimension(); 093 lu = matrix.getData(); 094 pivot = new int[m]; 095 cachedL = null; 096 cachedU = null; 097 cachedP = null; 098 099 // Initialize permutation array and parity 100 for (int row = 0; row < m; row++) { 101 pivot[row] = row; 102 } 103 even = true; 104 singular = false; 105 106 // Loop over columns 107 for (int col = 0; col < m; col++) { 108 109 // upper 110 for (int row = 0; row < col; row++) { 111 final double[] luRow = lu[row]; 112 double sum = luRow[col]; 113 for (int i = 0; i < row; i++) { 114 sum -= luRow[i] * lu[i][col]; 115 } 116 luRow[col] = sum; 117 } 118 119 // lower 120 int max = col; // permutation row 121 double largest = Double.NEGATIVE_INFINITY; 122 for (int row = col; row < m; row++) { 123 final double[] luRow = lu[row]; 124 double sum = luRow[col]; 125 for (int i = 0; i < col; i++) { 126 sum -= luRow[i] * lu[i][col]; 127 } 128 luRow[col] = sum; 129 130 // maintain best permutation choice 131 if (FastMath.abs(sum) > largest) { 132 largest = FastMath.abs(sum); 133 max = row; 134 } 135 } 136 137 // Singularity check 138 if (FastMath.abs(lu[max][col]) < singularityThreshold) { 139 singular = true; 140 return; 141 } 142 143 // Pivot if necessary 144 if (max != col) { 145 double tmp = 0; 146 final double[] luMax = lu[max]; 147 final double[] luCol = lu[col]; 148 for (int i = 0; i < m; i++) { 149 tmp = luMax[i]; 150 luMax[i] = luCol[i]; 151 luCol[i] = tmp; 152 } 153 int temp = pivot[max]; 154 pivot[max] = pivot[col]; 155 pivot[col] = temp; 156 even = !even; 157 } 158 159 // Divide the lower elements by the "winning" diagonal elt. 160 final double luDiag = lu[col][col]; 161 for (int row = col + 1; row < m; row++) { 162 lu[row][col] /= luDiag; 163 } 164 } 165 } 166 167 /** 168 * Returns the matrix L of the decomposition. 169 * <p>L is a lower-triangular matrix</p> 170 * @return the L matrix (or null if decomposed matrix is singular) 171 */ 172 public RealMatrix getL() { 173 if ((cachedL == null) && !singular) { 174 final int m = pivot.length; 175 cachedL = MatrixUtils.createRealMatrix(m, m); 176 for (int i = 0; i < m; ++i) { 177 final double[] luI = lu[i]; 178 for (int j = 0; j < i; ++j) { 179 cachedL.setEntry(i, j, luI[j]); 180 } 181 cachedL.setEntry(i, i, 1.0); 182 } 183 } 184 return cachedL; 185 } 186 187 /** 188 * Returns the matrix U of the decomposition. 189 * <p>U is an upper-triangular matrix</p> 190 * @return the U matrix (or null if decomposed matrix is singular) 191 */ 192 public RealMatrix getU() { 193 if ((cachedU == null) && !singular) { 194 final int m = pivot.length; 195 cachedU = MatrixUtils.createRealMatrix(m, m); 196 for (int i = 0; i < m; ++i) { 197 final double[] luI = lu[i]; 198 for (int j = i; j < m; ++j) { 199 cachedU.setEntry(i, j, luI[j]); 200 } 201 } 202 } 203 return cachedU; 204 } 205 206 /** 207 * Returns the P rows permutation matrix. 208 * <p>P is a sparse matrix with exactly one element set to 1.0 in 209 * each row and each column, all other elements being set to 0.0.</p> 210 * <p>The positions of the 1 elements are given by the {@link #getPivot() 211 * pivot permutation vector}.</p> 212 * @return the P rows permutation matrix (or null if decomposed matrix is singular) 213 * @see #getPivot() 214 */ 215 public RealMatrix getP() { 216 if ((cachedP == null) && !singular) { 217 final int m = pivot.length; 218 cachedP = MatrixUtils.createRealMatrix(m, m); 219 for (int i = 0; i < m; ++i) { 220 cachedP.setEntry(i, pivot[i], 1.0); 221 } 222 } 223 return cachedP; 224 } 225 226 /** 227 * Returns the pivot permutation vector. 228 * @return the pivot permutation vector 229 * @see #getP() 230 */ 231 public int[] getPivot() { 232 return pivot.clone(); 233 } 234 235 /** 236 * Return the determinant of the matrix 237 * @return determinant of the matrix 238 */ 239 public double getDeterminant() { 240 if (singular) { 241 return 0; 242 } else { 243 final int m = pivot.length; 244 double determinant = even ? 1 : -1; 245 for (int i = 0; i < m; i++) { 246 determinant *= lu[i][i]; 247 } 248 return determinant; 249 } 250 } 251 252 /** 253 * Get a solver for finding the A × X = B solution in exact linear 254 * sense. 255 * @return a solver 256 */ 257 public DecompositionSolver getSolver() { 258 return new Solver(lu, pivot, singular); 259 } 260 261 /** Specialized solver. */ 262 private static class Solver implements DecompositionSolver { 263 264 /** Entries of LU decomposition. */ 265 private final double[][] lu; 266 267 /** Pivot permutation associated with LU decomposition. */ 268 private final int[] pivot; 269 270 /** Singularity indicator. */ 271 private final boolean singular; 272 273 /** 274 * Build a solver from decomposed matrix. 275 * @param lu entries of LU decomposition 276 * @param pivot pivot permutation associated with LU decomposition 277 * @param singular singularity indicator 278 */ 279 private Solver(final double[][] lu, final int[] pivot, final boolean singular) { 280 this.lu = lu; 281 this.pivot = pivot; 282 this.singular = singular; 283 } 284 285 /** {@inheritDoc} */ 286 public boolean isNonSingular() { 287 return !singular; 288 } 289 290 /** {@inheritDoc} */ 291 public RealVector solve(RealVector b) { 292 final int m = pivot.length; 293 if (b.getDimension() != m) { 294 throw new DimensionMismatchException(b.getDimension(), m); 295 } 296 if (singular) { 297 throw new SingularMatrixException(); 298 } 299 300 final double[] bp = new double[m]; 301 302 // Apply permutations to b 303 for (int row = 0; row < m; row++) { 304 bp[row] = b.getEntry(pivot[row]); 305 } 306 307 // Solve LY = b 308 for (int col = 0; col < m; col++) { 309 final double bpCol = bp[col]; 310 for (int i = col + 1; i < m; i++) { 311 bp[i] -= bpCol * lu[i][col]; 312 } 313 } 314 315 // Solve UX = Y 316 for (int col = m - 1; col >= 0; col--) { 317 bp[col] /= lu[col][col]; 318 final double bpCol = bp[col]; 319 for (int i = 0; i < col; i++) { 320 bp[i] -= bpCol * lu[i][col]; 321 } 322 } 323 324 return new ArrayRealVector(bp, false); 325 } 326 327 /** {@inheritDoc} */ 328 public RealMatrix solve(RealMatrix b) { 329 330 final int m = pivot.length; 331 if (b.getRowDimension() != m) { 332 throw new DimensionMismatchException(b.getRowDimension(), m); 333 } 334 if (singular) { 335 throw new SingularMatrixException(); 336 } 337 338 final int nColB = b.getColumnDimension(); 339 340 // Apply permutations to b 341 final double[][] bp = new double[m][nColB]; 342 for (int row = 0; row < m; row++) { 343 final double[] bpRow = bp[row]; 344 final int pRow = pivot[row]; 345 for (int col = 0; col < nColB; col++) { 346 bpRow[col] = b.getEntry(pRow, col); 347 } 348 } 349 350 // Solve LY = b 351 for (int col = 0; col < m; col++) { 352 final double[] bpCol = bp[col]; 353 for (int i = col + 1; i < m; i++) { 354 final double[] bpI = bp[i]; 355 final double luICol = lu[i][col]; 356 for (int j = 0; j < nColB; j++) { 357 bpI[j] -= bpCol[j] * luICol; 358 } 359 } 360 } 361 362 // Solve UX = Y 363 for (int col = m - 1; col >= 0; col--) { 364 final double[] bpCol = bp[col]; 365 final double luDiag = lu[col][col]; 366 for (int j = 0; j < nColB; j++) { 367 bpCol[j] /= luDiag; 368 } 369 for (int i = 0; i < col; i++) { 370 final double[] bpI = bp[i]; 371 final double luICol = lu[i][col]; 372 for (int j = 0; j < nColB; j++) { 373 bpI[j] -= bpCol[j] * luICol; 374 } 375 } 376 } 377 378 return new Array2DRowRealMatrix(bp, false); 379 } 380 381 /** {@inheritDoc} */ 382 public RealMatrix getInverse() { 383 return solve(MatrixUtils.createRealIdentityMatrix(pivot.length)); 384 } 385 } 386 }