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

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

**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**. and U as a series of rows

_{1}**U**,....,

_{1}**U**.

_{1}To start the derivation, set

**L**equal to the first column of A, divided through by a

_{1}_{11}and

**U**equal to the first row.

_{1}Now construct a new matrix, A

^{2}= A -

**L**, where

_{1}U_{1}**L**denotes a vector product (i.e. a matrix where the i

_{1}U_{1}^{th}row is

**U**multiplied by the i

_{1}^{th}entry of

**L**; see the example if this is unclear!)

_{1}Now

**L**is the second column of A

_{2}^{2}divided by A

^{2}

_{2,2}; and

**U**is the second row of A

_{2}^{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**i.e. has 1 as the n

_{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 1This gives

**L**= [1,0,-1]

_{1}^{T}/1 and

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

_{1}^{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 3This gives us

**L**= [0,4,8]

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

^{T}and

**U**= [0,4,0]. We now construct A

_{2}^{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 3From which we get

**L**= [0,0,3]

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

^{T}and

**U**= [0,0,3]. This enables us to build our solution:

_{3}1 0 0 1 3 2 L= 0 1 0 and U = 0 4 0 -1 2 1 0 0 3You can multiply these out to check that we indeed obtain A.