2.1 LUP Decomposition

The basic idea of LUP decomposition, is that given a non-singular (invertible) square matrix M , we want to find matrices L , U and P such that,

Where L is unit lower triangular, U is upper triangular and P is a permutation matrix. A lower triangular matrix has all entries above the main diagonal 0, and an upper triangular matrix has all entries below the main diagonal 0. A unit lower or upper triangular matrix is lower or upper triangular with all entries on the main diagonal equal to 1. For example, a unit lower triangular matrix,

and an upper triangular matrix,

A permutation matrix has one entry per row that is equal to 1, and the rest equal to 0. When pre-multiplying a matrix M with a permutation matrix P , the rows of M are swapped. For example,

We want to find L and U to allow us to solve the system by using backwards, then forwards substitution. Given a system of equations in upper triangular form, we can solve using backwards substitution as follows,

And forwards substitution with a lower triangular works the same way, only we start from the first row of the matrix, not the last. With a unit lower triangular matrix, we know the entries main diagonal are 1, so it's even simpler,

Once L , U and P are found (and they always exist for non-singular square matrices, though they may not be unique), it becomes very simple to solve MV = C .

Let Y = UV and solve LY = PC for Y using forward substitution. Then solve UV = Y for V using backward substitution.

So now we need to find L , U and P . This is done using Crout's method .

First we'll look at the simpler case of LU decomposition (no permutation matrix).

Let M be an invertible n×n matrix.

Where

and

is called the Schur complement of M with respect to m₁₁. If M is invertible, then so is the Schur complement (see [CLRS01] for proof). This then gives us a recursive way to find the LU decomposition. Let L ' U ' equal the Schur complement, where L' is unit lower triangular and U' is upper triangular. Then we have,

To implement this algorithm, the tail recursion of decomposing the Schur complement into L ' U ' can be converted into iteration.

The LU decomposition will fail if there exists i in {1,..., n } such that m_ii = 0. Even if no such i exists, the LU decomposition will still have numerical stability problems if any m_ii is a very small value.

LUP decomposition permutes the matrix so that at every step we are dividing by the largest possible value. This additional step of finding the best (largest) value and swapping rows appropriately, is known as pivoting. If no non-zero value is found in the first column of any Schur complement, then the matrix M is not invertible.

Listing 3 implements LUP decomposition, and using the decomposition to solve systems of linear equations. The permutation matrix P is stored as a permutation array, since this is a more compact representation for exactly the same information.

2.2 Calculating the New Shadow Caster Volume Planes