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.

A pair of matrices L,U∈RNXN comprise a LU factorisation of A∈RNXN if the following conditions are met

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 L1,...,L1. and U as a series of rows U1,....,U1.
To start the derivation, set L1 equal to the first column of A, divided through by a11 and U1 equal to the first row.
Now construct a new matrix, A2 = A - L1U1, where L1U1 denotes a vector product (i.e. a matrix where the ith row is U1 multiplied by the ith entry of L1; see the example if this is unclear!)
Now L2 is the second column of A2 divided by A22,2; and U2 is the second row of A2.
Continue in this vein to obtain the third column and row through use of A3 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 eN i.e. has 1 as the nth 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 2N3/3 in terms of computing complexity- compare to the order N2 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 L1= [1,0,-1]T/1 and U1= [1,3,2]. We now construct A2:
 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 L2= [0,4,8]T/4 = [0,1,2]T and U2= [0,4,0]. We now construct A2:
 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 L3= [0,0,3]T/3 = [0,0,1]T and U3= [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.