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 A

**x**=

**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

^{-1}**b**. 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

**Solving a system with LU factorisation**

Given a matrix system A**x**=**b**, and an LU factorisation of A, proceed as follows:

- Use back substitution to solve L
**y**=**b**
- Use back substitution to solve U
**x**=**y**

Then

**x** is precisely that desired, since A

**x**=(LU)

**x**=L(U

**x**)=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**_{1},...,**L**_{1}. and U as a series of rows **U**_{1},....,**U**_{1}.

To start the derivation, set **L**_{1} equal to the first column of A, divided through by a_{11} and **U**_{1} equal to the first row.

Now construct a new matrix, A^{2} = A - **L**_{1}U_{1}, where **L**_{1}U_{1} denotes a vector product (i.e. a matrix where the i^{th} row is **U**_{1} multiplied by the i^{th} entry of **L**_{1}; see the example if this is unclear!)

Now **L**_{2} is the second column of A^{2} divided by A^{2}_{2,2}; and **U**_{2} is the second row of A^{2}.

Continue in this vein to obtain the third column and row through use of A^{3} 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}/3 in terms of computing complexity- compare to the order N^{2} 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}= [1,0,-1]

^{T}/1 and

**U**_{1}= [1,3,2]. We now construct A

^{2}:

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**_{2}= [0,4,8]

^{T}/4 = [0,1,2]

^{T} and

**U**_{2}= [0,4,0]. We now construct A

^{2}:

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**_{3}= [0,0,3]

^{T}/3 = [0,0,1]

^{T} and

**U**_{3}= [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.