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 }