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