/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.linear;

import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.util.FastMath;

Calculates the LUP-decomposition of a square matrix.

The LUP-decomposition of a matrix A consists of three matrices L, U and P that satisfy: P×A = L×U. L is lower triangular (with unit diagonal terms), U is upper triangular and P is a permutation matrix. All matrices are m×m.

As shown by the presence of the P matrix, this decomposition is implemented using partial pivoting.

This class is based on the class with similar name from the JAMA library.

  • a getP method has been added,
  • the det method has been renamed as getDeterminant,
  • the getDoublePivot method has been removed (but the int based getPivot method has been kept),
  • the solve and isNonSingular methods have been replaced by a getSolver method and the equivalent methods provided by the returned DecompositionSolver.
See Also:
Since:2.0 (changed to concrete class in 3.0)
/** * Calculates the LUP-decomposition of a square matrix. * <p>The LUP-decomposition of a matrix A consists of three matrices L, U and * P that satisfy: P&times;A = L&times;U. L is lower triangular (with unit * diagonal terms), U is upper triangular and P is a permutation matrix. All * matrices are m&times;m.</p> * <p>As shown by the presence of the P matrix, this decomposition is * implemented using partial pivoting.</p> * <p>This class is based on the class with similar name from the * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p> * <ul> * <li>a {@link #getP() getP} method has been added,</li> * <li>the {@code det} method has been renamed as {@link #getDeterminant() * getDeterminant},</li> * <li>the {@code getDoublePivot} method has been removed (but the int based * {@link #getPivot() getPivot} method has been kept),</li> * <li>the {@code solve} and {@code isNonSingular} methods have been replaced * by a {@link #getSolver() getSolver} method and the equivalent methods * provided by the returned {@link DecompositionSolver}.</li> * </ul> * * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a> * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a> * @since 2.0 (changed to concrete class in 3.0) */
public class LUDecomposition {
Default bound to determine effective singularity in LU decomposition.
/** Default bound to determine effective singularity in LU decomposition. */
private static final double DEFAULT_TOO_SMALL = 1e-11;
Entries of LU decomposition.
/** Entries of LU decomposition. */
private final double[][] lu;
Pivot permutation associated with LU decomposition.
/** Pivot permutation associated with LU decomposition. */
private final int[] pivot;
Parity of the permutation associated with the LU decomposition.
/** Parity of the permutation associated with the LU decomposition. */
private boolean even;
Singularity indicator.
/** Singularity indicator. */
private boolean singular;
Cached value of L.
/** Cached value of L. */
private RealMatrix cachedL;
Cached value of U.
/** Cached value of U. */
private RealMatrix cachedU;
Cached value of P.
/** Cached value of P. */
private RealMatrix cachedP;
Calculates the LU-decomposition of the given matrix. This constructor uses 1e-11 as default value for the singularity threshold.
Params:
  • matrix – Matrix to decompose.
Throws:
/** * Calculates the LU-decomposition of the given matrix. * This constructor uses 1e-11 as default value for the singularity * threshold. * * @param matrix Matrix to decompose. * @throws NonSquareMatrixException if matrix is not square. */
public LUDecomposition(RealMatrix matrix) { this(matrix, DEFAULT_TOO_SMALL); }
Calculates the LU-decomposition of the given matrix.
Params:
  • matrix – The matrix to decompose.
  • singularityThreshold – threshold (based on partial row norm) under which a matrix is considered singular
Throws:
/** * Calculates the LU-decomposition of the given matrix. * @param matrix The matrix to decompose. * @param singularityThreshold threshold (based on partial row norm) * under which a matrix is considered singular * @throws NonSquareMatrixException if matrix is not square */
public LUDecomposition(RealMatrix matrix, double singularityThreshold) { if (!matrix.isSquare()) { throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension()); } final int m = matrix.getColumnDimension(); lu = matrix.getData(); pivot = new int[m]; cachedL = null; cachedU = null; cachedP = null; // Initialize permutation array and parity for (int row = 0; row < m; row++) { pivot[row] = row; } even = true; singular = false; // Loop over columns for (int col = 0; col < m; col++) { // upper for (int row = 0; row < col; row++) { final double[] luRow = lu[row]; double sum = luRow[col]; for (int i = 0; i < row; i++) { sum -= luRow[i] * lu[i][col]; } luRow[col] = sum; } // lower int max = col; // permutation row double largest = Double.NEGATIVE_INFINITY; for (int row = col; row < m; row++) { final double[] luRow = lu[row]; double sum = luRow[col]; for (int i = 0; i < col; i++) { sum -= luRow[i] * lu[i][col]; } luRow[col] = sum; // maintain best permutation choice if (FastMath.abs(sum) > largest) { largest = FastMath.abs(sum); max = row; } } // Singularity check if (FastMath.abs(lu[max][col]) < singularityThreshold) { singular = true; return; } // Pivot if necessary if (max != col) { double tmp = 0; final double[] luMax = lu[max]; final double[] luCol = lu[col]; for (int i = 0; i < m; i++) { tmp = luMax[i]; luMax[i] = luCol[i]; luCol[i] = tmp; } int temp = pivot[max]; pivot[max] = pivot[col]; pivot[col] = temp; even = !even; } // Divide the lower elements by the "winning" diagonal elt. final double luDiag = lu[col][col]; for (int row = col + 1; row < m; row++) { lu[row][col] /= luDiag; } } }
Returns the matrix L of the decomposition.

L is a lower-triangular matrix

Returns:the L matrix (or null if decomposed matrix is singular)
/** * Returns the matrix L of the decomposition. * <p>L is a lower-triangular matrix</p> * @return the L matrix (or null if decomposed matrix is singular) */
public RealMatrix getL() { if ((cachedL == null) && !singular) { final int m = pivot.length; cachedL = MatrixUtils.createRealMatrix(m, m); for (int i = 0; i < m; ++i) { final double[] luI = lu[i]; for (int j = 0; j < i; ++j) { cachedL.setEntry(i, j, luI[j]); } cachedL.setEntry(i, i, 1.0); } } return cachedL; }
Returns the matrix U of the decomposition.

U is an upper-triangular matrix

Returns:the U matrix (or null if decomposed matrix is singular)
/** * Returns the matrix U of the decomposition. * <p>U is an upper-triangular matrix</p> * @return the U matrix (or null if decomposed matrix is singular) */
public RealMatrix getU() { if ((cachedU == null) && !singular) { final int m = pivot.length; cachedU = MatrixUtils.createRealMatrix(m, m); for (int i = 0; i < m; ++i) { final double[] luI = lu[i]; for (int j = i; j < m; ++j) { cachedU.setEntry(i, j, luI[j]); } } } return cachedU; }
Returns the P rows permutation matrix.

P is a sparse matrix with exactly one element set to 1.0 in each row and each column, all other elements being set to 0.0.

The positions of the 1 elements are given by the pivot permutation vector.

See Also:
Returns:the P rows permutation matrix (or null if decomposed matrix is singular)
/** * Returns the P rows permutation matrix. * <p>P is a sparse matrix with exactly one element set to 1.0 in * each row and each column, all other elements being set to 0.0.</p> * <p>The positions of the 1 elements are given by the {@link #getPivot() * pivot permutation vector}.</p> * @return the P rows permutation matrix (or null if decomposed matrix is singular) * @see #getPivot() */
public RealMatrix getP() { if ((cachedP == null) && !singular) { final int m = pivot.length; cachedP = MatrixUtils.createRealMatrix(m, m); for (int i = 0; i < m; ++i) { cachedP.setEntry(i, pivot[i], 1.0); } } return cachedP; }
Returns the pivot permutation vector.
See Also:
Returns:the pivot permutation vector
/** * Returns the pivot permutation vector. * @return the pivot permutation vector * @see #getP() */
public int[] getPivot() { return pivot.clone(); }
Return the determinant of the matrix
Returns:determinant of the matrix
/** * Return the determinant of the matrix * @return determinant of the matrix */
public double getDeterminant() { if (singular) { return 0; } else { final int m = pivot.length; double determinant = even ? 1 : -1; for (int i = 0; i < m; i++) { determinant *= lu[i][i]; } return determinant; } }
Get a solver for finding the A × X = B solution in exact linear sense.
Returns:a solver
/** * Get a solver for finding the A &times; X = B solution in exact linear * sense. * @return a solver */
public DecompositionSolver getSolver() { return new Solver(lu, pivot, singular); }
Specialized solver.
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
Entries of LU decomposition.
/** Entries of LU decomposition. */
private final double[][] lu;
Pivot permutation associated with LU decomposition.
/** Pivot permutation associated with LU decomposition. */
private final int[] pivot;
Singularity indicator.
/** Singularity indicator. */
private final boolean singular;
Build a solver from decomposed matrix.
Params:
  • lu – entries of LU decomposition
  • pivot – pivot permutation associated with LU decomposition
  • singular – singularity indicator
/** * Build a solver from decomposed matrix. * @param lu entries of LU decomposition * @param pivot pivot permutation associated with LU decomposition * @param singular singularity indicator */
private Solver(final double[][] lu, final int[] pivot, final boolean singular) { this.lu = lu; this.pivot = pivot; this.singular = singular; }
{@inheritDoc}
/** {@inheritDoc} */
public boolean isNonSingular() { return !singular; }
{@inheritDoc}
/** {@inheritDoc} */
public RealVector solve(RealVector b) { final int m = pivot.length; if (b.getDimension() != m) { throw new DimensionMismatchException(b.getDimension(), m); } if (singular) { throw new SingularMatrixException(); } final double[] bp = new double[m]; // Apply permutations to b for (int row = 0; row < m; row++) { bp[row] = b.getEntry(pivot[row]); } // Solve LY = b for (int col = 0; col < m; col++) { final double bpCol = bp[col]; for (int i = col + 1; i < m; i++) { bp[i] -= bpCol * lu[i][col]; } } // Solve UX = Y for (int col = m - 1; col >= 0; col--) { bp[col] /= lu[col][col]; final double bpCol = bp[col]; for (int i = 0; i < col; i++) { bp[i] -= bpCol * lu[i][col]; } } return new ArrayRealVector(bp, false); }
{@inheritDoc}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b) { final int m = pivot.length; if (b.getRowDimension() != m) { throw new DimensionMismatchException(b.getRowDimension(), m); } if (singular) { throw new SingularMatrixException(); } final int nColB = b.getColumnDimension(); // Apply permutations to b final double[][] bp = new double[m][nColB]; for (int row = 0; row < m; row++) { final double[] bpRow = bp[row]; final int pRow = pivot[row]; for (int col = 0; col < nColB; col++) { bpRow[col] = b.getEntry(pRow, col); } } // Solve LY = b for (int col = 0; col < m; col++) { final double[] bpCol = bp[col]; for (int i = col + 1; i < m; i++) { final double[] bpI = bp[i]; final double luICol = lu[i][col]; for (int j = 0; j < nColB; j++) { bpI[j] -= bpCol[j] * luICol; } } } // Solve UX = Y for (int col = m - 1; col >= 0; col--) { final double[] bpCol = bp[col]; final double luDiag = lu[col][col]; for (int j = 0; j < nColB; j++) { bpCol[j] /= luDiag; } for (int i = 0; i < col; i++) { final double[] bpI = bp[i]; final double luICol = lu[i][col]; for (int j = 0; j < nColB; j++) { bpI[j] -= bpCol[j] * luICol; } } } return new Array2DRowRealMatrix(bp, false); }
Get the inverse of the decomposed matrix.
Throws:
Returns:the inverse matrix.
/** * Get the inverse of the decomposed matrix. * * @return the inverse matrix. * @throws SingularMatrixException if the decomposed matrix is singular. */
public RealMatrix getInverse() { return solve(MatrixUtils.createRealIdentityMatrix(pivot.length)); } } }