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.