8Numerical linear algebra
IB Numerical Analysis
8.2 LU factorization
In general, we don’t always have triangular matrices. The idea is, for every
matrix A, find some lower and triangular matrices L, U such that A = LU.
=
If we can do this, then we can solve
Ax
=
b
in two steps — we first find a
y
such that Ly = b. Then we find an x such that Ux = y. Then
Ax = LUx = Ly = b.
So what we want is such a factorization. To guarantee uniqueness of the
factorization, we require that
L
is unit, i.e. the diagonals are all 1. Otherwise,
given any such factorization, we can divide
U
by some (non-zero) constant and
multiply L by the same constant to get another factorization.
Definition
(LU factorization)
. A
=
LU
is an LU factorization if
U
is upper
triangular and L is unit lower triangular (i.e. the diagonals of L are all 1).
Note that since
L
has to be unit, it must be non-singular. However, we still
allow A and U to be singular. Note that
det(A) = det(L) det(U ) = det(U ).
So A is singular if and only if U is.
Unfortunately, even if
A
is non-singular, it may not have an LU factorization.
We take a simple example of
0 1
1 1
This is clearly non-singular, with determinant
−
1. However, we can manually
check that there is no LU factorization of this.
On the other hand, while we don’t really like singular matrices, singular
matrices can still have LU factorizations. For example,
0 0
0 1
=
1 0
0 1
0 0
0 1
is trivially an LU factorization of a singular matrix.
Fortunately, if
A
is non-singular and has an LU factorization, then this
factorization is unique (this is not necessarily true if
A
is singular. For example,
we know
0 0
0 1
=
1 0
a 1
0 0
0 1
for any real number a).
To understand when LU factorizations exist, we try to construct one, and
see when it could fail.
We write L in terms of its columns, and U in terms of its rows:
L =
l
1
l
2
··· l
n
, U =
u
T
1
u
T
2
.
.
.
u
T
n
Clearly, these rows and columns cannot be arbitrary, since
L
and
U
are triangular.
In particular, l
i
, u
i
must be zero in the first i −1 entries.
Suppose this is an LU factorization of A. Then we can write
A = L · U = l
1
u
T
1
+ l
2
u
T
2
+ ··· + l
n
u
T
n
.
What do these matrices look like? For each
i
, we know
l
i
and
u
i
have the first
i −
1 entries zero. So the first
i −
1 rows and columns of
l
i
u
T
i
are zero. In
particular, the first row and columns only have contributions from
l
1
u
T
1
, the
second row/column only has contributions from l
2
u
T
2
and l
1
u
T
1
etc.
The plan is as follows:
(i)
Obtain
l
1
and
u
1
from the first row and column of
A
. Since the first entry
of
l
1
is 1,
u
T
1
is exactly the first row of
A
. We can then obtain
l
2
by taking
the first column of A and dividing by U
11
= A
11
.
(ii) Obtain l
2
and u
2
form the second row and column of A − l
1
u
T
1
similarly.
(iii) ···
(iv) Obtain l
n
and u
T
n
from the nth row and column of A −
P
n−1
i=1
l
i
u
T
i
.
We can turn this into an algorithm. We define the intermediate matrices, starting
with
A
(0)
= A.
For k = 1, ··· , n, we let
U
kj
= A
(k−1)
kj
j = k, ··· , n
L
ik
=
A
ik
A
(k−1)
kk
i = k, ··· , n
A
(k)
ij
= A
(k−1)
ij
− L
ik
U
kj
i, j ≥ k
When
k
=
n
, we end up with a zero matrix, and then
U
and
L
are completely
filled.
We can now see when this will break down. A sufficient condition for
A
=
LU
to exist is that
A
(k−1)
kk
6
= 0 for all
k
. Since
A
(k−1)
kk
=
U
kk
, this sufficient condition
ensures
U
, and hence
A
is non-singular. Conversely, if
A
is non-singular and
an LU factorization exists, then this would always work, since we must have
A
(k−1)
kk
=
U
kk
6
= 0. Moreover, the LU factorization must be given by this
algorithm. So we get uniqueness.
The problem with this sufficient condition is that most of these coefficients
do not appear in the matrix
A
. They are constructed during the algorithm. We
don’t know easily what they are in terms of the coefficients of
A
. We will later
come up with an equivalent condition on our original A that is easier to check.
Note that as long as this method does not break down, we need
O
(
n
3
)
operations to perform this factorization. Recall we only needed
O
(
n
2
) operations
to solve the equation after factorization. So the bulk of the work in solving
Ax = b is in doing the LU factorization.
As before, this allows us to find the inverse of
A
if it is non-singular. In
particular, solving
Ax
j
=
e
j
gives the
j
th column of
A
−1
. Note that we are
solving the system for the same
A
for each
j
. So we only have to perform the
LU factorization once, and then solve
n
different equations. So in total we need
O(n
3
) operations.
However, we still have the problem that factorization is not always possible.
Requiring that we must factorize
A
as
LU
is too restrictive. The idea is to factor
something closely related, but not exactly A. Instead, we want a factorization
P A = LU,
where
P
is a permutation matrix. Recall a permutation matrix acting on a
column vector just permutes the elements of the vector, and
P
acting on
A
would just permute the rows of
A
. So we want to factor
A
up to a permutation
of rows.
Our goal now is to extend the previous algorithm to allow permutations of
rows, and then we shall show that we will be able to perform this factorization
all the time.
Suppose our breakdown occurs at
k
= 1, i.e.
A
(0)
11
=
A
11
= 0. We find a
permutation matrix
P
1
and let it act via
P
1
A
(0)
. The idea is to look down the
first column of
A
, and find a row starting with a non-zero element, say
p
. Then
we use
P
1
to interchange rows 1 and
p
such that
P
1
A
(0)
has a non-zero top-most
entry. For simplicity, we assume we always need a
P
1
, and if
A
(0)
11
is non-zero in
the first place, we just take P
1
to be the identity.
After that, we can carry on. We construct
l
1
and
u
1
from
P
1
A
(0)
as before,
and set A
(1)
= P
1
A
(0)
− l
1
u
T
1
.
But what happens if the first column of
A
is completely zero? Then no
interchange will make the (1
,
1) entry non-zero. However, in this case, we don’t
actually have to do anything. We can immediately find our
l
1
and
u
1
, namely
set
l
1
=
e
1
(or anything) and let
u
T
1
be the first row of
A
(0)
. Then this already
works. Note however that this corresponds to
A
(and hence
U
) being singular,
and we are not too interested with these.
The later steps are exactly analogous. Suppose we have
A
(k−1)
kk
= 0. Again
we find a
P
k
such that
P
k
A
(k−1)
has a non-zero (
k, k
) entry. We then construct
l
k
and u
k
from P
k
A
(k−1)
and set
A
(k)
= P
k
A
(n−1)
− l
k
u
t
k
.
Again, if the
k
th column of
A
(k−1)
is completely zero, we set
l
k
=
e
k
and
u
T
k
to
be the kth row of A
(k−1)
. But again this implies A and U will be singular.
However, as we do this, the permutation matrices appear all over the place
inside the algorithm. It is not immediately clear that we do get a factorization
of the form
P A
=
LU
. Fortunately, keeping track of the interchanges, we do
have an LU factorization
P A =
˜
LU,
where U is what we got from the algorithm,
P = P
n−1
···P
2
P
1
,
while
˜
L is given by
˜
L =
˜
l
1
···
˜
l
n
,
˜
l
k
= P
n−1
···P
k−1
l
k
.
Note that in particular, we have
˜
l
n−1
= l
n−1
,
˜
l
n
= l
n
.
One problem we have not considered is the problem of inexact arithmetic. While
these formula are correct mathematically, when we actually implement things,
we do them on computers with finite precision. As we go through the algorithm,
errors will accumulate, and the error might be amplified to a significant amount
when we get the reach the end. We want an algorithm that is insensitive to
errors. In order to work safely in inexact arithmetic, we will put the element of
largest modulus in the (
k, k
)th position, not just an arbitrary non-zero one, as
this minimizes the error when dividing.