174 6 Linear Algebra
LU decomposition forms the backbone of other algorithms in linear algebra, such as the
solution of linear equations given by
a 11 x 1 +a 12 x 2 +a 13 x 3 +a 14 x 4 =w 1
a 21 x 1 +a 22 x 2 +a 23 x 3 +a 24 x 4 =w 2
a 31 x 1 +a 32 x 2 +a 33 x 3 +a 34 x 4 =w 3
a 41 x 1 +a 42 x 2 +a 43 x 3 +a 44 x 4 =w 4.The above set of equations is conveniently solved by using LUdecomposition as an interme-
diate step, see the next subsection for more details on how tosolve linear equations with an
LU decomposed matrix.
The matrixA∈Rn×nhas an LU factorization if the determinant is different fromzero. If
the LU factorization exists andAis non-singular, then the LU factorization is unique and the
determinant is given by
det{A}=u 11 u 22 ...unn.
For a proof of this statement, see chapter 3.2 of Ref. [28].
The algorithm for obtainingLandUis actually quite simple. We start always with the first
column. In our simple ( 4 × 4 ) case we obtain then the following equations for the first column
a 11 = u 11
a 21 =l 21 u 11
a 31 =l 31 u 11
a 41 =l 41 u 11 ,which determine the elementsu 11 ,l 21 ,l 31 andl 41 inLandU. Writing out the equations for the
second column we get
a 12 = u 12
a 22 = l 21 u 12 +u 22
a 32 =l 31 u 12 +l 32 u 22
a 42 =l 41 u 12 +l 42 u 22.
Here the unknowns areu 12 ,u 22 ,l 32 andl 42 which can all be evaluated by means of the
results from the first column and the elements ofA. Note an important feature. When going
from the first to the second column we do not need any further information from the matrix
elementsai 1. This is a general property throughout the whole algorithm.Thus the memory
locations for the matrixAcan be used to store the calculated matrix elements ofLandU.
This saves memory.
We can generalize this procedure into three equations
i<j: li 1 u 1 j+li 2 u 2 j+···+liiui j=ai j
i=j: li 1 u 1 j+li 2 u 2 j+···+liiuj j=ai j
i>j: li 1 u 1 j+li 2 u 2 j+···+li juj j=ai jwhich gives the following algorithm:
Calculate the elements inLandUcolumnwise starting with column one. For each column
(j):
- Compute the first elementu 1 jby
u 1 j=a 1 j. - Next, we calculate all elementsui j,i= 2 ,...,j− 1