LU Decomposition

This documentation was generated with the assistance of AI. Please report any inaccuracies.

LUDecomposer factors a matrix A into a lower triangular matrix L (elements only on the diagonal and below) and an upper triangular matrix U (elements only on the diagonal and above), so that (eq. 2.3.1):

Once this factorization is available, solving A * x = b reduces to two cheap triangular solves: forward substitution to solve L * y = b (eq. 2.3.6), followed by backsubstitution to solve U * x = y (eq. 2.3.7):

Because the decomposition itself, rather than the right-hand side, does the expensive work, the same LUDecomposer instance can solve for any number of right-hand sides, one at a time, without repeating the factorization; and inverting a matrix or computing its determinant reduces to cheap extra steps once the factorization is stored.

The algorithm

LUDecomposer computes the factorization with Crout’s algorithm, which solves for the entries of L and of U using the recurrences (eq. 2.3.12-2.3.13), fixing :

Applying these column by column, for , first solving for and then for , only ever needs entries already computed. This means each entry of A is read exactly once and can be overwritten in place by the corresponding or : since Crout’s method fixes the diagonal of L to 1, both triangular factors are packed into a single matrix the same size as A (eq. 2.3.14):

Numerical Recipes' Figure 2.3.1 shows the order in which this combined matrix is filled: by columns, from left to right, and within each column from top to bottom. The diagram below reproduces that figure for the 4x4 case, labeling each entry with the step at which it is computed and highlighting, for two typical entries, exactly which already-computed entries (shaded) they depend on:

Order in which Crout’s algorithm fills the combined L/U matrix one column at a time top to bottom within each column

Pivoting is essential for numerical stability: at each column, the algorithm selects the largest-magnitude candidate in that column (partial pivoting, i.e., row interchanges only, scaled as if every row had been normalized to a unit maximum coefficient — implicit pivoting) and promotes it to the diagonal before continuing. LUDecomposer therefore does not decompose A itself but a row-wise permutation of it, and keeps track of that permutation (getPivot()) so that it can be undone when solving or when reconstructing L.

LUDecomposer requires about a third of the operations of Gauss-Jordan elimination for solving one or a few right-hand sides, while matching it for matrix inversion; it is the preferred way to solve a linear system when A is square, which is why Utils.solve and Utils.inverse fall back to it whenever the input matrix is square.

Usage

import com.irurueta.algebra.LUDecomposer;
import com.irurueta.algebra.Matrix;

Matrix a = new Matrix(2, 2);
a.setElementAt(0, 0, 2.0);
a.setElementAt(0, 1, 1.0);
a.setElementAt(1, 0, 1.0);
a.setElementAt(1, 1, 3.0);

final var decomposer = new LUDecomposer(a);
decomposer.decompose();

Matrix l = decomposer.getL();
Matrix u = decomposer.getU();
Matrix pivottedL = decomposer.getPivottedL(); // L before pivot correction
int[] pivot = decomposer.getPivot();

Matrix b = Matrix.newFromArray(new double[]{5.0, 10.0});
Matrix x = decomposer.solve(b);

double determinant = decomposer.determinant();
boolean singular = decomposer.isSingular();

LUDecomposer does not expose a dedicated inverse method; use Utils.inverse, which relies on LUDecomposer internally for square matrices.

Once decompose() has run, solve can be called repeatedly with different right-hand sides without recomputing the factorization. determinant() implements eq. 2.3.15, multiplying the diagonal of U and the sign accumulated from row interchanges during pivoting:

and isSingular() checks whether any of those diagonal elements is (numerically) zero.

Error handling

decompose() throws NotReadyException if no input matrix has been set, LockedException if the instance is already decomposing, and DecomposerException if the input matrix has more columns than rows. Result accessors throw NotAvailableException if called before decompose(). solve, determinant() and isSingular() throw WrongSizeException if the input matrix is not square (solve additionally if b does not have a matching number of rows), and solve throws SingularMatrixException if the input matrix turns out to be singular.

Reference

LUDecomposer is based on the algorithm described in Numerical Recipes, 3rd Edition, section 2.3, "LU Decomposition and Its Applications".