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
-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∈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 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 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 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
L2= [0,4,8]
T/4 = [0,1,2]
T and
U2= [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
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.