logo
Tsung-Ming Huang
Department of Mathematics National Taiwan Normal University, Taiwan
E-mail: min@math.ntnu.edu.tw
September 12, 2015
1 / 89
logo
Outline
1 Linear systems of equations
2 Pivoting Strategies
3 Matrix factorization
4 Special types of matrices
logo
Linear systems of equations
Three operations to simplify the linear system:
1 (λEi) → (Ei): Equation Ei can be multiplied by λ 6= 0 with the resulting equation used in place of Ei.
2 (Ei+ λEj) → (Ei): Equation Ej can be multiplied by λ 6= 0 and added to equation Eiwith the resulting equation used in place of Ei.
3 (Ei) ↔ (Ej): Equation Ei and Ej can be transposed in order.
Example 1
E1 : x1 + x2 + 3x4 = 4,
E2 : 2x1 + x2 − x3 + x4 = 1, E3 : 3x1 − x2 − x3 + 2x4 = −3, E4 : −x1 + 2x2 + 3x3 − x4 = 4.
3 / 89
logo
Solution:
(E2− 2E1) → (E2), (E3− 3E1) → (E3)and (E4+ E1) → (E4):
E1: x1 + x2 + 3x4 = 4,
E2: − x2 − x3 − 5x4 = −7, E3: − 4x2 − x3 − 7x4 = −15, E4: 3x2 + 3x3 + 2x4 = 8.
(E3− 4E2) → (E3)and (E4+ 3E2) → (E4):
E1: x1 + x2 + 3x4 = 4,
E2: − x2 − x3 − 5x4 = −7,
E3: 3x3 + 13x4 = 13,
E4: − 13x4 = −13.
logo
Backward-substitution process:
1 E4 ⇒ x4= 1
2 Solve E3for x3: x3=1
3(13 − 13x4) =1
3(13 − 13) = 0.
3 E2gives
x2= −(−7 + 5x4+ x3) = −(−7 + 5 + 0) = 2.
4 E1gives
x1= 4 − 3x4− x2= 4 − 3 − 2 = −1.
5 / 89
logo
Solve linear systems of equations
a11x1+ a12x2+ · · · + a1nxn = b1
a21x1+ a22x2+ · · · + a2nxn = b2
...
an1x1+ an2x2+ · · · + annxn = bn
Rewrite in the matrix form
Ax = b, (1)
where
A =
a11 a12 · · · a1n a21 a22 · · · a2n
... ... . .. ... an1 an2 · · · ann
, b =
b1 b2
... bn
, x =
x1 x2
... xn
and [A, b] is called the augmented matrix.
logo
Gaussian elimination with backward substitution
The augmented matrix in previous example is
1 1 0 3 4
2 1 −1 1 1
3 −1 −1 2 −3
−1 2 3 −1 4
.
(E2− 2E1) → (E2), (E3− 3E1) → (E3)and (E4+ E1) → (E4):
1 1 0 3 4
0 −1 −1 −5 −7
0 −4 −1 −7 −15
0 3 3 2 8
.
(E3− 4E2) → (E3)and (E4+ 3E2) → (E4):
1 1 0 3 4
0 −1 −1 −5 −7
0 0 3 13 13
0 0 0 −13 −13
.
7 / 89
logo
The general Gaussian elimination procedure
Provided a116= 0, for each i = 2, 3, . . . , n,
Ei− ai1 a11
E1
→ (Ei).
Transform all the entries in the first col. below the diagonal are zero. Denote the new entry in the ith row and jth col. by aij. For i = 2, 3 . . . , n − 1, providedaii 6= 0,
Ej−aji
aiiEi
→ (Ej), ∀ j = i + 1, i + 2, . . . , n.
Transform all the entries in the ith column below the diagonal are zero.
Result anupper triangularmatrix:
a11 a12 · · · a1n b1
0 a22 · · · a2n b2
... . .. ... ... ... 0 · · · 0 ann bn
.
logo
The process of Gaussian elimination result in a sequence of matrices as follows:
A = A(1) → A(2)→ · · · → A(n)=upper triangular matrix The matrix A(k)has the following form:
A(k)=
a(1)11 · · · a(1)1,k−1 a(1)1k · · · a(1)1j · · · a(1)1n
... . .. ... ... ... ...
0 · · · a(k−1)k−1,k−1 a(k−1)k−1,k · · · a(k−1)k−1,j · · · a(k−1)k−1,n 0 · · · 0 a(k)kk · · · a(k)kj · · · a(k)kn
... ... ... ... ...
0 · · · 0 a(k)ik · · · a(k)ij · · · a(k)in
... ... ... ... ...
0 · · · 0 a(k)nk · · · a(k)nj · · · a(k)nn
9 / 89
logo
The entries of A(k)are produced by the formula
a(k)ij =
a(k−1)ij , for i = 1, . . . , k − 1, j = 1, . . . , n;
0, for i = k, . . . , n, j = 1, . . . , k − 1;
a(k−1)ij − a
(k−1) i,k−1
a(k−1)k−1,k−1 × a(k−1)k−1,j, for i = k, . . . , n, j = k, . . . , n.
The procedure will fail if one of the elements a(1)11, a(2)22, . . . , a(n)nn is zero.
a(i)ii is called the pivot element.
logo
Backward substitution
The new linear system is triangular:
a11x1 + a12x2 + · · · + a1nxn = b1, a22x2 + · · · + a2nxn = b2,
...
annxn = bn Solving the nth equation for xngives
xn= bn ann
.
Solving the (n − 1)th equation for xn−1and using the value for xnyields
xn−1= bn−1− an−1,nxn
an−1,n−1
. In general,
xi= bi−Pn
j=i+1aijxj aii
, ∀ i = n − 1, n − 2, . . . , 1.
11 / 89
logo
Algorithm 1 (Backward Substitution)
Suppose thatU ∈ Rn×n isnonsingular upper triangularand b ∈ Rn. This algorithm computes the solution ofU x = b.
For i = n, . . . , 1 tmp = 0
For j = i + 1, . . . , n
tmp = tmp + U (i, j) ∗ x(j) End for
x(i) = (b(i) − tmp)/U (i, i) End for
logo
Example 2
Solve system of linear equations.
6 −2 2 4
12 −8 6 10
3 −13 9 3
−6 4 1 −18
x1 x2
x3
x4
=
12 34 27
−38
Solution:
1st step Use 6 as pivot element, the first row as pivot row, and multipliers 2,12, −1are produced to reduce the system to
6 −2 2 4
0 −4 2 2
0 −12 8 1
0 2 3 −14
x1 x2
x3
x4
=
12 10 21
−26
13 / 89
logo
2nd step Use −4 as pivot element, the second row as pivot row, and multipliers 3, −12 are computed to reduce the system to
6 −2 2 4
0 −4 2 2
0 0 2 −5
0 0 4 −13
x1
x2 x3
x4
=
12 10
−9
−21
3rdstep Use 2 as pivot element, the third row as pivot row, and multipliers 2 is found to reduce the system to
6 −2 2 4
0 −4 2 2
0 0 2 −5
0 0 0 −3
x1 x2 x3
x4
=
12 10
−9
−3
logo
4thstep The backward substitution is applied:
x4 = −3
−3 = 1, x3 = −9 + 5x4
2 = −9 + 5 2 = −2, x2 = 10 − 2x4− 2x3
−4 = 10 − 2 + 4
−4 = −3, x1 = 12 − 4x4− 2x3+ 2x2
6 = 12 − 4 + 4 − 6
6 = 1.
This example is done since a(k)kk 6= 0 for all k = 1, 2, 3, 4.
How to do if a(k)kk = 0for some k?
15 / 89
logo
Example 3
Solve system of linear equations.
1 −1 2 −1 2 −2 3 −3
1 1 1 0
1 −1 4 3
x1 x2
x3
x4
=
−8
−20
−2 4
Solution:
1st step Use 1 as pivot element, the first row as pivot row, and multipliers 2, 1, 1 are produced to reduce the system to
1 −1 2 −1
0 0 −1 −1
0 2 −1 1
0 0 2 4
x1 x2
x3
x4
=
−8
−4 6 12
logo
2nd step Since a(2)22 = 0and a(2)32 6= 0, the operation
(E2) ↔ (E3)is performed to obtain a new system
1 −1 2 −1
0 2 −1 1
0 0 −1 −1
0 0 2 4
x1
x2 x3
x4
=
−8 6
−4 12
3rdstep Use −1 as pivot element, the third row as pivot row, and multipliers −2 is found to reduce the system to
1 −1 2 −1
0 2 −1 1
0 0 −1 −1
0 0 0 2
x1 x2
x3 x4
=
−8 6
−4 4
17 / 89
logo
4thstep The backward substitution is applied:
x4 = 4 2 = 2, x3 = −4 + x4
−1 = 2, x2 = 6 − x4+ x3
2 = 3,
x1 = −8 + x4− 2x3+ x2
1 = −7.
This example illustrates what is done if a(k)kk = 0for some k.
If a(k)pk 6= 0 for some p with k + 1 ≤ p ≤ n, then the operation (Ek) ↔ (Ep)is performed to obtain new matrix.
If a(k)pk = 0for each p, then the linear system does not have a unique solution and the procedure stops.
logo
Algorithm 2 (Gaussian elimination)
GivenA ∈ Rn×n andb ∈ Rn, this algorithm implements the Gaussian elimination procedure to reduceAtoupper triangular and modify the entries ofbaccordingly.
For k = 1, . . . , n − 1
Let p be the smallest integer with k ≤ p ≤ n and apk 6= 0.
If @ p, then stop.
If p 6= k, then perform (Ep) ↔ (Ek).
For i = k + 1, . . . , n t = A(i, k)/A(k, k) A(i, k) = 0
b(i) = b(i) − t × b(k) For j = k + 1, . . . , n
A(i, j) = A(i, j) − t × A(k, j) End for
End for
End for 19 / 89
logo
Number of floating-point arithmetic operations
Eliminate kth column For i = k + 1, . . . , n
t = A(i, k)/A(k, k); b(i) = b(i) − t × b(k).
For j = k + 1, . . . , n
A(i, j) = A(i, j) − t × A(k, j) End for
End for
Multiplications/divisions
(n − k) + (n − k) + (n − k)(n − k) = (n − k)(n − k + 2) Additions/subtractions
(n − k) + (n − k)(n − k) = (n − k)(n − k + 1)
logo
Total number of operations for multiplications/divisions
n−1
X
k=1
(n − k)(n − k + 2) =
n−1
X
k=1
(n2− 2nk + k2+ 2n − 2k)
= (n2+ 2n)
n−1
X
k=1
1 − 2(n + 1)
n−1
X
k=1
k +
n−1
X
k=1
k2
= (n2 + 2n)(n − 1) − 2(n + 1)(n − 1)n
2 +(n − 1)n(2n − 1) 6
= 2n3+ 3n2− 5n
6 .
Total number of operations for additions/subtractions
n−1
X
k=1
(n − k)(n − k + 1) =
n−1
X
k=1
(n2− 2nk + k2+ n − k)
= (n2+ n)
n−1
X
k=1
1 − (2n + 1)
n−1
X
k=1
k +
n−1
X
k=1
k2= n3− n 3 .
21 / 89
logo
Backward substitution
x(n) = b(n)/U (n, n).
For i = n − 1, . . . , 1
tmp = U (i, i + 1) × x(i + 1) For j = i + 2, . . . , n
tmp = tmp + U (i, j) × x(j) End for
x(i) = (b(i) − tmp)/U (i, i) End for
Multiplications/divisions 1 +
n−1
X
i=1
[(n − i) + 1] =n2+ n 2 Additions/subtractions
n−1
X
i=1
[(n − i − 1) + 1] =n2− n 2
logo
The total number of arithmetic operations in Gaussian elimination with backward substitution is:
Multiplications/divisions 2n3+ 3n2− 5n
6 + n2+ n 2 = n3
3 + n2−n 3 ≈ n3
3 Additions/subtractions
n3− n
3 +n2− n 2 = n3
3 + n2 2 −5n
6 ≈ n3 3
23 / 89
logo
Exercise
Page 368: 5, 10, 12, 15
logo
Pivoting Strategies
If a(k)kk is small in magnitude compared to a(k)jk , then
|mjk| =
a(k)jk a(k)kk
> 1.
Round-off error introduced in the computation of a(k+1)j` = a(k)j` − mjka(k)k`, for ` = k + 1, . . . , n.
Error can be increased when performing the backward substitution for
xk = bk−Pn
j=k+1a(k)kjxj a(k)kk
with a small value of a(k)kk.
25 / 89
logo
Example 4 The linear system
E1 : 0.003000x1 + 59.14x2 = 59.17, E2 : 5.291x1 − 6.130x2 = 46.78,
has the exact solution x1 = 10.00and x2 = 1.000. Suppose Gaussian elimination is performed on this system using four-digit arithmetic with rounding.
a11= 0.0030is small and m21= 5.291
0.0030 = 1763.6¯6 ≈ 1764.
Perform (E2− m21E1) → (E2):
0.0030x1 + 59.14x2 = 59.17
− 104309.37¯6x2 = −104309.37¯6.
logo
Rounding with four-digit arithmetic:
Coefficient of x2:
−6.130 − 1764 × 59.14 = −6.130 − 104322.96
≈ −6.130 − 104300 = −104306.13
≈ −104300.
Right hand side:
46.78 − 1764 × 59.17 = 46.78 − 104375.88
≈ 46.78 − 104400 = −104353.22
≈ −104400.
New linear system:
0.0030x1 + 59.14x2 = 59.17
− 104300x2 ≈ −104400.
27 / 89
logo
Approximated solution:
x2 = 104400
104300 ≈ 1.001, x1 = 59.17 − 59.14 × 1.001
0.0030 = 59.17 − 59.19914 0.0030
≈ 59.17 − 59.20
0.0030 = −10.00.
This ruins the approximation to the actual value x1 = 10.00.
logo
Partial pivoting
To avoid the pivot element small relative to other entries, pivoting is performed by selecting an element a(k)pq with a larger magnitude as the pivot.
Specifically, select pivoting a(k)pk with
|a(k)pk| = max
k≤i≤n|a(k)ik | and perform (Ek) ↔ (Ep).
This row interchange strategy is called partial pivoting.
29 / 89
logo
Example 5
Reconsider the linear system
E1 : 0.003000x1 + 59.14x2 = 59.17, E2 : 5.291x1 − 6.130x2 = 46.78.
Find pivoting with
max{|a11|, |a21|} = 5.291 = |a21|.
Perform (E2) ↔ (E1):
E1 : 5.291x1 − 6.130x2 = 46.78, E2 : 0.003000x1 + 59.14x2 = 59.17.
The multiplier for new system is m21= a21
a11
= 0.0005670.
logo
The operation (E2− m21E1) → (E2)reduces the system to 5.291x1 − 6.130x2 = 46.78,
59.14x2 ≈ 59.14.
The four-digit answers resulting from the backward substitution are the correct values x1 = 10.00and x2= 1.000.
31 / 89
logo
Example 6 The linear system
E1: 30.00x1 + 591400x2 = 591700, E2: 5.291x1 − 6.130x2 = 46.78,
is the same as that in previous example except that all the entries in the first equation have been multiplied by 104. The pivoting is a11= 30.00and the multiplier
m21= 5.291
30.00 = 0.1764 leads to the system
30.00x1 + 591400x2 = 591700
− 104300x2 ≈ −104400,
which has inaccurate solution x2 ≈ 1.001 and x1 ≈ −10.00.
logo
Scaled partial pivoting
Define a scale factor si as si = max
1≤j≤n|aij|, for i = 1, . . . , n.
If si= 0for some i, then the system has no unique solution.
In the ith column, choose the least integer p ≥ i with
|api|
sp = max
i≤k≤n
|aki| sk and perform (Ei) ↔ (Ep)if p 6= i.
The scale factors s1, . . . , snare computed only once and must also be interchanged when row interchanges are performed.
33 / 89
logo
Example 7
Apply scaled partial pivoting to the linear system E1: 30.00x1 + 591400x2 = 591700, E2: 5.291x1 − 6.130x2 = 46.78.
The scale factors s1and s2are
s1= max{|30.00|, |591400|} = 591400 and
s2= max{|5.291|, | − 6.130|} = 6.130.
Consequently,
|a11| s1
= 30.00
591400 = 0.5073 × 10−4,
|a21|
s2 = 5.291
6.130 = 0.8631, and the interchange (E1) ↔ (E2)is made.
logo
Applying Gaussian elimination to the new system 5.291x1 − 6.130x2 = 46.78, 30.00x1 + 591400x2 = 591700
produces the correct results: x1 = 10.00and x2 = 1.000.
35 / 89
logo
Exercise
Page 379: 2, 4, 6, 31
logo
Matrix factorization
This equation has a unique solutionx = A−1bwhen the coefficient matrix A isnonsingular.
Use Gaussian elimination to factor the coefficient matrix into a product of matrices. The factorization is called LU-factorizationand has the formA = LU, whereLisunit lower triangularandU isupper triangular.
The solution to the original problemAx = LU x = bis then found by a two-step triangular solve process:
Ly = b, U x = y.
LU factorization requiresO(n3)arithmetic operations.
Forward substitution for solving a lower-triangular system Ly = brequiresO(n2). Backward substitution for solving an upper-triangular system U x = y requiresO(n2) arithmetic operations.
37 / 89
logo
A =
1 1 0 3
2 1 −1 1
3 −1 −1 2
−1 2 3 −1
⇒ A1 := L1A ≡
1 0 0 0
−2 1 0 0
−3 0 1 0 1 0 0 1
A =
1 1 0 3
0 −1 −1 −5 0 −4 −1 −7
0 3 3 2
⇒ A2 := L2A1 ≡
1 0 0 0
0 1 0 0
0 −4 1 0
0 3 0 1
A1 =
1 1 0 3
0 −1 −1 −5
0 0 3 13
0 0 0 −13
= L2L1A
logo
We have
A = L−11 L−12 A2= LR.
where L and R are lower and upper triangular, respectively.
Question
How to compute L−11 and L−12 ?
L1 =
1 0 0 0
−2 1 0 0
−3 0 1 0 1 0 0 1
= I −
0 2 3
−1
1 0 0 0
L2 =
1 0 0 0
0 1 0 0
0 −4 1 0
0 3 0 1
= I −
0 0 4
−3
0 1 0 0
39 / 89
logo
Since
I −
0 2 3
−1
1 0 0 0
I +
0 2 3
−1
1 0 0 0
= I,
we have
L−11 =
1 0 0 0
−2 1 0 0
−3 0 1 0 1 0 0 1
−1
=
1 0 0 0 2 1 0 0 3 0 1 0
−1 0 0 1
logo
Since
I −
0 0 4
−3
0 1 0 0
I +
0 0 4
−3
0 1 0 0
= I,
we have
L−12 =
1 0 0 0
0 1 0 0
0 −4 1 0
0 3 0 1
−1
=
1 0 0 0
0 1 0 0
0 4 1 0
0 −3 0 1
41 / 89
logo
By the fact
L−11 L−12 =
1 0 0 0 2 1 0 0 3 0 1 0
−1 0 0 1
1 0 0 0
0 1 0 0
0 4 1 0
0 −3 0 1
=
1 0 0 0
2 1 0 0
3 4 1 0
−1 −3 0 1
,
it holds that
1 1 0 3
2 1 −1 1
3 −1 −1 2
−1 2 3 −1
=
1 0 0 0
2 1 0 0
3 4 1 0
−1 −3 0 1
1 1 0 3
0 −1 −1 −5
0 0 3 13
0 0 0 −13
.
logo
For a given vectorv ∈ Rnwithvk6= 0for some 1 ≤ k ≤ n, let
`ik = vi
vk, i = k + 1, . . . , n,
`k=
0 · · · 0 `k+1,k · · · `n,k T
, and
Mk= I − `keTk =
1 · · · 0 0 · · · 0 ... . .. ... ... ... 0 · · · 1 0 · · · 0 0 · · · −`k+1,k 1 · · · 0 ... ... ... . .. ...
0 · · · −`n,k 0 · · · 1
.
43 / 89
logo
Then one can verify that Mkv =
v1 · · · vk 0 · · · 0 T
.
Mk is called aGaussian transformation, the vector`kaGauss vector. Furthermore, one can verify that
Mk−1 = (I − `keTk)−1 = I + `keTk =
1 · · · 0 0 · · · 0 ... . .. ... ... ... 0 · · · 1 0 · · · 0 0 · · · `k+1,k 1 · · · 0 ... ... ... . .. ...
0 · · · `n,k 0 · · · 1
.
logo
Given a nonsingular matrix A ∈ Rn×n, denote A(1)≡ [a(1)ij ] = A.
Ifa(1)11 6= 0, then
M1 = I − `1eT1, where
`1=
0 `21 · · · `n1
T
, `i1= a(1)i1
a(1)11, i = 2, . . . , n, can be formed such that
A(2) = M1A(1) =
a(1)11 a(1)12 · · · a(1)1n 0 a(2)22 · · · a(2)2n ... ... . .. ... 0 a(2)n2 · · · a(2)nn
,
where
a(2)ij = a(1)ij − `i1× a(1)1j , for i = 2, . . . , n and j = 2, . . . , n.
45 / 89
logo
In general, at the k-th step, we are confronted with a matrix A(k) = Mk−1· · · M2M1A(1)
=
a(1)11 a(1)12 · · · a(1)1,k−1 a(1)1k · · · a(1)1n 0 a(2)22 · · · a(2)2,k−1 a(2)2k · · · a(2)2n ... ... . .. ... ... ... 0 0 · · · a(k−1)k−1,k−1 a(k−1)k−1,k · · · a(k−1)k−1,n 0 0 · · · 0 a(k)kk · · · a(k)kn
... ... ... ... . .. ... 0 0 · · · 0 a(k)kn · · · a(k)nn
.
If the pivota(k)kk 6= 0, then the multipliers
`ik = a(k)ik a(k)kk
, i = k + 1, . . . , n,
logo
can be computed and the Gaussian transformation Mk = I − `keTk, where `k=
0 · · · 0 `k+1,k · · · `nk
T
, can be applied to the left of A(k)to obtain
A(k+1) = MkA(k)
=
a(1)11 a(1)12 · · · a(1)1,k−1 a(1)1k a(1)1,k+1 · · · a(1)1n 0 a(2)22 · · · a(2)2,k−1 a(2)2k a(2)2,k+1 · · · a(2)2n ... ... . .. ... ... ... ... 0 0 · · · a(k−1)k−1,k−1 a(k−1)k−1,k a(k−1)k−1,k+1 · · · a(k−1)k−1,n 0 0 · · · 0 a(k)kk a(k)k,k+1 · · · a(k)kn
... ... ... 0 a(k+1)k+1,k+1 · · · a(k+1)k+1,n
... ... ... ... ... ...
0 0 · · · 0 0 a(k+1)n,k+1 · · · a(k+1)nn
,
47 / 89
logo
in which
a(k+1)ij = a(k)ij − `ika(k)kj , (2) for i = k + 1, . . . , n, j = k + 1, . . . , n. Upon the completion,
U ≡ A(n)= Mn−1· · · M2M1A isupper triangular. Hence
A = M1−1M2−1· · · Mn−1−1 U ≡ LU,
logo
where
L ≡ M1−1· · · Mn−1−1 = (I − `1eT1)−1(I − `2eT2)−1· · · (I − `n−1eTn−1)−1
= (I + `1eT1)(I + `2eT2) · · · (I + `n−1eTn−1)
= I + `1eT1 + `2eT2 + · · · + `n−1eTn−1
=
1 0 0 · · · 0
`21 1 0 · · · 0
`31 `32 1 · · · 0 ... ... ... . .. ...
`n1 `n2 `n3 · · · 1
isunit lower triangular. This matrix factorization is called the LU-factorizationofA.
49 / 89
logo
Algorithm 3 (LU Factorization)
Given anonsingular squarematrixA ∈ Rn×n, this algorithm computes aunit lower triangularmatrixLand anupper triangularmatrixU such thatA = LU. The matrix A is overwrittenby L and U .
For k = 1, . . . , n − 1 For i = k + 1, . . . , n
A(i, k) = A(i, k)/A(k, k) For j = k + 1, . . . , n
A(i, j) = A(i, j) − A(i, k) × A(k, j) End for
End for End for
logo
Forward Substitution
When a linear systemLx = bislower triangularof the form
`11 0 · · · 0
`21 `22 · · · 0 ... ... . .. ...
`n1 `n2 · · · `nn
x1 x2
... xn
=
b1 b2
... bn
,
where all diagonals`ii6= 0, xican be obtained by the following procedure
x1 = b1/`11,
x2 = (b2− `21x1)/`22,
x3 = (b3− `31x1− `32x2)/`33, ...
xn = (bn− `n1x1− `n2x2− · · · − `n,n−1xn−1)/`nn.
51 / 89
logo
The general formulation for computing xi is xi =
bi−
i−1
X
j=1
`ijxj
`ii, i = 1, 2, . . . , n.
Algorithm 4 (Forward Substitution)
Suppose thatL ∈ Rn×n isnonsingular lower triangularand b ∈ Rn. This algorithm computes the solution of Lx = b.
For i = 1, . . . , n tmp = 0
For j = 1, . . . , i − 1
tmp = tmp + L(i, j) ∗ x(j) End for
x(i) = (b(i) − tmp)/L(i, i) End for
logo
Example 8
E1 : x1 + x2 + 3x4 = 4,
E2 : 2x1 + x2 − x3 + x4 = 1, E3 : 3x1 − x2 − x3 + 2x4 = −3, E4 : −x1 + 2x2 + 3x3 − x4 = 4.
Solution:
The sequence {(E2− 2E1) → (E2), (E3− 3E1) → (E3), (E4− (−1)E1) → (E4), (E3− 4E2) → (E3),
(E4− (−3)E2) → (E4)}converts the system to the triangular system
x1 + x2 + 3x4 = 4,
− x2 − x3 − 5x4 = −7, 3x3 + 13x4 = 13,
− 13x4 = −13.
53 / 89
logo
LU factorization of A:
A =
1 1 0 3
2 1 −1 1
3 −1 −1 2
−1 2 3 −1
=
1 0 0 0
2 1 0 0
3 4 1 0
−1 −3 0 1
1 1 0 3
0 −1 −1 −5
0 0 3 13
0 0 0 −13
= LU.
logo
Solve Ly = b:
1 0 0 0
2 1 0 0
3 4 1 0
−1 −3 0 1
y1
y2 y3 y4
=
8 7 14
−7
which implies that y1 = 8,
y2 = 7 − 2y1 = −9, y3 = 14 − 3y1− 4y2= 26, y4 = −7 + y1+ 3y2 = −26.
55 / 89
logo
Solve U x = y:
1 1 0 3
0 −1 −1 −5
0 0 3 13
0 0 0 −13
x1 x2
x3
x4
=
8
−9 26
−26
which implies that x4 = 2,
x3 = (26 − 13x4)/3 = 0,
x2 = (−9 + 5x4+ x3)/(−1) = −1, x1 = 8 − 3x4− x2 = 3.
logo
Partial pivoting
At the k-th step, select pivoting a(k)pk with
|a(k)pk| = max
k≤i≤n|a(k)ik |
and perform (Ek) ↔ (Ep). That is, choose a permutation matrix
Pk=
Ik−1 0 0 0 0
0 0 0 1 0
0 0 Ip−k−1 0 0
0 1 0 0 0
0 0 0 0 In−p
so that
(PkA(k))kk
= max
k≤i≤n
(A(k))ik and
A(k+1)= M(k)PkA(k).
57 / 89
logo
LetP1, . . . , Pk−1 be thepermutationschosen andM1, . . . Mk−1 denote theGaussian transformationsperformed in the first k − 1steps. At the k-th step, a permutation matrixPk is chosen so that
|(PkMk−1· · · M1P1A)kk| = max
k≤i≤n|(Mk−1· · · M1P1A)ik| . As a consequence,|`ij| ≤ 1for i = 1, . . . , n, j = 1, . . . , i. Upon completion, we obtain an upper triangular matrix
U ≡ Mn−1Pn−1· · · M1P1A. (3) Since anyPkissymmetricandPkTPk= Pk2 = I, we have
Mn−1Pn−1· · · M2P2M1P2· · · Pn−1Pn−1· · · P2P1A = U, therefore,
Pn−1· · · P1A = (Mn−1Pn−1· · · M2P2M1P2· · · Pn−1)−1U.
logo
In summary,Gaussian elimination with partial pivotingleads to theLU factorization
P A = LU, (4)
where
P = Pn−1· · · P1 is a permutation matrix, and
L ≡ (Mn−1Pn−1· · · M2P2M1P2· · · Pn−1)−1
= Pn−1· · · P2M1−1P2M2−1· · · Pn−1Mn−1−1 . Since
Pj =
Ij−1 0 0 0 0
0 0 0 1 0
0 0 Ip−j−1 0 0
0 1 0 0 0
0 0 0 0 In−p
, `j =
0
... 0
`j+1,j ...
`nj
,
59 / 89
logo
it implies that for i < j,
eTi Pj = eTi , eTi `j = 0, Pj`i =
0 · · · 0 `˜i+1,i · · · `˜n,i T
≡ ˜`i,
⇒
P2M1−1P2 = P2(I + `1eT1)P2 = I + ˜`1eT1
⇒
P2M1−1P2M2−1 = (I + ˜`1eT1)(I + `2eT2) = I + ˜`1eT1 + `2eT2,
⇒
P3 P2M1−1P2M2−1 P3= I + ˆ`1eT1 + ˜`2eT2
⇒ · · ·
Therefore,Lisunit lower triangular.
logo
Algorithm 5 (LU -factorization with Partial Pivoting)
Given a nonsingular A ∈ Rn×n, this algorithm finds apermutationP, and computes aunit lower triangularLand anupper triangularU such thatP A = LU. A isoverwrittenby L and U , and P is not formed. Aninteger arraypis instead used for storing the row/column indices.
p(1 : n) = 1 : n For k = 1, . . . , n − 1
m = k
For i = k + 1, . . . , n
If|A(p(m), k)| < |A(p(i), k)|, thenm = i End For
` = p(k); p(k) = p(m); p(m) = ` For i = k + 1, . . . , n
A(p(i), k) = A(p(i), k)/A(p(k), k) For j = k + 1, . . . , n
A(p(i), j) = A(p(i), j) − A(p(i), k)A(p(k), j) End For
End For
End For 61 / 89
