LU factorisation (thing)

A significant proportion of the world's computing capacity and time is dedicated to the solution of systems of linear equations through manipulation of matrices, and as a result many methods have emerged to tackle the essential problem: finding a vector x that satisfies the equation Ax=b. A costly approach that is manageable for student problems (but rarely real-world ones) is to find the inverse of the matrix, as then pre-multiplication on each side yields x=A^-1b. Indeed, there are many ways to go about finding the inverse, from the motivating but painfully slow example of finding adjugate matrices to row operations; but ultimately this approach is overkill. A system can be solved without finding the inverse by means of back substitution, a method that obtains the elements of x based on eliminating successive rows in a triangular matrix. Of course, to apply back substitution you'll need a triangular matrix first, but if your initial A isn't of this type, all is not lost. A common means of obtaining a suitable matrix is to apply Gaussian elimination; but the LU factorisation offers another route which may save on computation, be easier to understand, or simpler to implement in code.

Definition
A pair of matrices L,U∈R^NXN comprise a LU factorisation of A∈R^NXN if the following conditions are met

• L is lower triangular (no entries above the main diagonal)
• L_ii=1 ∀i∈1...N (entries on the main diagonal are all 1)
• U is upper triangular (no entries below the main diagonal)
• LU=A (L and U are actually factors!)

Solving a system with LU factorisation
Given a matrix system Ax=b, and an LU factorisation of A, proceed as follows:

• Use back substitution to solve Ly=b
• Use back substitution to solve Ux=y

Then x is precisely that desired, since Ax=(LU)x=L(Ux)=L(y) by the second equation = b by the first equation.

Deriving an LU factorisation
We derive L and U iteratively, one row or column at a time. For the derivation, consider L as a series of columns L₁,...,L₁. and U as a series of rows U₁,....,U₁.
To start the derivation, set L₁ equal to the first column of A, divided through by a₁₁ and U₁ equal to the first row.
Now construct a new matrix, A² = A - L₁U₁, where L₁U₁ denotes a vector product (i.e. a matrix where the i^th row is U₁ multiplied by the i^th entry of L₁; see the example if this is unclear!)
Now L₂ is the second column of A² divided by A²_2,2; and U₂ is the second row of A².
Continue in this vein to obtain the third column and row through use of A³ etc. until all N rows and columns of L and U have been found (you should get that the final column of L is simply the elementary column vector e_N i.e. has 1 as the n^th entry and zero elsewhere, due to the definition of L)

Complications and Complexity
It is worth noting that this procedure fails if any entry on the diagonal of A is zero, since then we are required to divide through by that zero to obtain a column of L, which is obviously invalid. The same scenario will also break Gaussian elimination, as it introduces a zero pivot. To avoid such problems, rows of the matrix would need to be reordered and the problem manipulated to compensate for these row operations. Testing for and resolving this issue adds to the overheads of an implementation of either method.

An efficient implementation of LU factorisation (using knowledge that many entries in L and U are zero, and the diagonal of L contains only ones, for instance) is of order 2N³/3 in terms of computing complexity- compare to the order N² calculations required for back substitution if we are fortunate enough to start with a triangular matrix. Although some savings can be made in special cases (See the Cholesky factorisation of SPD matrices), as N grows large such precise methods become unwieldy and iterative techniques need to be employed.

---

A sample LU factorisation
Let A be the 3X3 matrix

 1 3 2
 0 4 0
-1 5 1

This gives L₁= [1,0,-1]^T/1 and U₁= [1,3,2]. We now construct A²:

 1 3 2   |1| [1,3,2]   1 3 2      1  3  2     0 0 0 
 0 4 0  -|0|        =  0 4 0   -  0  0  0  =  0 4 0
-1 5 1   |-1|         -1 5 1     -1 -3 -2     0 8 3

This gives us L₂= [0,4,8]^T/4 = [0,1,2]^T and U₂= [0,4,0]. We now construct A²:

 0 0 0   |0| [0,4,0]   0 0 0      0  0  0     0 0 0 
 0 4 0  -|1|        =  0 4 0   -  0  4  0  =  0 0 0
 0 8 3   |2|           0 8 3      0  8  0     0 0 3

From which we get L₃= [0,0,3]^T/3 = [0,0,1]^T and U₃= [0,0,3]. This enables us to build our solution:

   1 0 0          1 3 2
L= 0 1 0  and U = 0 4 0
  -1 2 1          0 0 3

You can multiply these out to check that we indeed obtain A.