handoutB.dvi
ECS130 Scientific Computation Handout B January 13, 2017
1. Gaussian elimination = LU factorization
A = LU.
where L is a unit lower triangular matrix and U a upper triangular matrix.
2. Not all matrices have the LU factorization. For example,
A =
0 1 0
0 2 1
2 3 1
6= LU.
3. A permutation matrix P is an identity matrix with permuted rows.
Let P, P1, P2 be n× n permutation matrices, and X be an n× n matrix. Then
• P TP = I, i.e., P−1 = P T .
• det(P ) = ±1.
• P1P2 is also a permutation matrix.
• PX is the same as X with its rows permuted.
• XP is the same as X with its columns permuted.
• P1XP2 reorders both rows and columns of X.
4. The need of pivoting, mathematically
The LU factorization can fail on nonsingular matrices, see the above example. But by exchanging
the first and third rows, we get
PA =
2 3 1
0 2 1
0 1 0
=
1 0 0
0 1 0
0 1
2
1
2 3 1
0 2 1
0 0 −1
2
:= LU.
5. The above simple observation is the basis for LU factorization with pivoting.
Theorem. If A is nonsingular, then there exist permutations P , a unit lower triangular matrix
L, and a nonsingular upper triangular matrix U such that
PA = LU.
6. Function lutx.m
7. Solving Ax = b using the LU factorization
1. Factorize A into PA = LU
2. Permute the entries of b: b := Pb.
3. Solve L(Ux) = b for Ux by forward substitution:
Ux = L−1b.
4. Solve Ux = L−1b for x by back substitution:
x = U−1(L−1b).
1
8. Function bslashx0.m
9. The need for pivoting, numerically
Let us apply LU factorization without pivoting to
A =
[
.0001 1
1 1
]
=
[
10−4 1
1 1
]
= LU =
[
1
l21 1
] [
u11 u12
u22
]
in three decimal-digit floating point arithmetic. We obtain
L =
[
1 0
fl(1/10−4) 1
]
=
[
1 0
104 1
]
,
U =
[
10−4 1
fl(1− 104 · 1)
]
=
[
10−4 1
−104
]
,
so
LU =
[
1 0
104 1
] [
10−4 1
−104
]
=
[
10−4 1
1 0
]
6≈ A,
where the original a22 has been entirely “lost” from the computation by subtracting 10
4 from it.
In fact, we would have gotten the same LU factors whether a22 had been 1, 0,−2, or any number
such that fl(a22 − 10
4) = −104. Since the algorithm proceeds to work only with L and U , it will
get the same answer for all these different a22, which correspond to completely different A and so
completely different x = A−1b; there is no way to guarantee an accurate answer. This is called
numerical instability. L and U are not the exact factors of a matrix close to A.
Let us see what happens when we go on to solve Ax = [1, 2]T for x using this LU factorization.
The correct answer is x ≈ [1, 1]T . Instead we get the following. Solving
Ly =
[
1
2
]
⇒ y1 = fl(1/1) = 1 and y2 = fl(2− 10
4 · 1) = −104.
Note that the value 2 has been “lost” by subtracting 104 from it. Solving
Ux = y =
[
1
−104
]
⇒ x̂2 = fl((−10
4)/(−104)) = 1 and x̂1 = fl((1− 1)/10
−4) = 0,
a completely erroneous solution.
On the other hand, the LU factorization with partial pivoting would have reversed the order of
the two equations before proceeding. You can confirm that we get
PA = LU,
where
P =
[
0 1
1 0
]
, L =
[
1 0
fl(.0001/1) 1
]
=
[
1 0
.0001 1
]
,
and
U =
[
1 1
fl(1− .0001 · 1)
]
=
[
1 1
1
]
.
The computed LU approximates A very accurately. As a result, the computed solution x is also
perfect!
2