LU decomposition, animated

Matrix A

Pivot matrix P

Instructions

L matrix is in lower left, U matrix is top right, and PA (original matrix multiplied by pivot) is lower right. This layout shows the interaction of each element in PA with the rows and columns of each factor L and U.

The "Slower" button toggles between three animation modes: fast, slow and ultra-slow. There is also a step-by-step mode: press "Step-by-step" button to enter this mode, and keep clicking it to advance the calculation. (Some steps don't change the screen at all; keep clicking!)

When the LUP decomposition is not unique, some elements in in L are filled with 0.666, making them easy to spot. (Such elements could be filled with any other value and the decomposition would still be valid.)

To "pivotize" is to exchange some rows of original matrix, in a way that diagonal elements are the biggest number in its own column.

Press the "Pivotize" button to find a permutation matrix, or press "Do not pivotize" to make P=identity and not make any changes. The default status is non-pivotized.

Important: if you don't press "Pivotize" and P is the identity matrix, you are actually calculating LU decomposition, not LUP. If the given matrix has no LU decomposition, some L or U elements will show "Error", and this means that LU has failed.

In the other hand, LUP decomposition cannot fail; just press "Pivotize" to find P for the current matrix, or press "Automatic" to keep P always updated for LUP.

Theory

Generally speaking, a matrix decomposition algorithm takes a matrix and tries to find a set of two or more matrixes, whose recombination is equal to the original. In the case of LU decomposition:

A = L.U

There are many types of decomposition, but the general idea is to decompose into factors that are "simpler" under some criteria. In the case of LU decomposition, L and U are triangular matrixes (L is lower-triangular and U is upper-triangular).

LU decomposition was invented by Alan Turing. It is useful when it is desirable to deal only with triangular matrixes. For example, it is very easy to calculate the determinant of a triangular matrix: just multiply the diagonal elements. Decomposing a matrix and then calculating det(L).det(U) is way more efficient than trying to calculate det(A) directly.

Not all matrixes can be factored in two triangular factors. Some divisions by zero may appear in the process. On the other hand, if we are allowed to rearrange rows first, the decomposition is always possible. The LU augmented by rearranging rows, or "pivoting", is the LUP decomposition.

P.A=L.U
A = P^-1.L.U

where P is a permutation matrix (an identity matrix with rows permuted), that can do row permutations under ordinary matrix multiplication. Every square matrix has a LUP decomposition. Invertible matrixes have unique LU factors. The LUP technique also improves numerical stability.

Thanks to the desirable features, the LUP decomposition is employed more often than pure LU. In most instances, references to "LU decomposition" actually mean "LUP decomposition".

Singular matrixes (e.g. all-zeros) can be LUP-decomposed but the result is not unique. Some implementations break under singular matrixes because of 0/0 divisons. Our version detects this condition and adopts the value 0.666 as an (arbitrary) result of the indetermination.

Algorithm

The Javascript of this page is not obfuscated, but it has been "asynchronized" in order to support animation, and it is not beautiful nor very readable. The straightforward form of the algorithm is this:

// LU decomposition algorithm

for (var j = 1; j <= n; ++j) {
	var i, k;
	for (var i = 1; i < (j + 1); ++i) {
		var s1 = 0;
		for (k = 1; k < i; ++k) {
			s1 += get("U", k, j) * get("L", i, k);
		}
		set("U", i, j, get("PA", i, j) - s1);
	}

	for (i = j + 1; i <= n; ++i) {
		var s2 = 0;
		for (k = 1; k < j; ++k) {
			s2 += get("U", k, j) * get("L", i, k);
		}

		var dividend = get("PA", i, j) - s2;
		var divisor = get("U", j, j);
		var res = dividend / divisor;

		if (res !== res) {
			res = 0.666;
		}

		set("L", i, j, res);
	}
}

The LU decompositon expects that the PA matrix is properly pivoted; otherwise, some division by zero may occur. The algorithm below finds the best permutation P for a given matrix A:

// Pivotize algoritm

var exchanges = 0;

for (var j = 1; j <= n; ++j) {
	var row = j;
	var val = 0;
	for (var i = j; i <= n; ++i) {
		var cand = Math.abs(get("A", i, j));
			console.log(cand);
		if (val < cand) {
			val = cand;
			row = i;
		}
	}
	if (j !== row) {
		++exchanges;
		for (var i = 1; i <= n; ++i) {
			var tmp = get("P", j, i);
			set("P", j, i, get("P", row, i));
			set("P", row, i, tmp);
		}
	}
}

Both code snippets are loosely based on Python implementation found at RosettaCode.