Linear Equations-Direct Methods

Download as pdf or txt
Download as pdf or txt
You are on page 1of 50

Chapter 5

Systems of Linear Equations: Direct Methods

Computational Methods
(SECE)

Yonas Y.
Computational Methods (SECE) Chapter 5 Yonas Y. 1 / 50
Outline

1 Linear Systems: Direct Methods


Introduction

2 Gaussian elimination and backward substitution


Backward substitution
Forward substitution
Gaussian elimination
LU decomposition

3 Pivoting strategies
Partial pivoting
GEPP stability
Matrices requiring no pivoting

4 The Cholesky decomposition

5 Estimating errors and the condition number


Sensitivity on matrix coefficients

Computational Methods (SECE) Chapter 5 Yonas Y. 2 / 50


Linear Systems: Direct Methods Introduction

Introduction
Kirchhoff’s laws of electrical circuits state that both the net flow of current
through each junction and the net voltage drop around each closed loop of a
circuit are zero.
Suppose we have a circuit shown in the diagram below.

Using G as a reference point, Kirchhoff’s laws imply that the currents satisfy
the following system of linear equations:
5i1 + 5i2 = V ,
i3 − i4 − i5 = 0,
2i4 − 3i5 = 0,
i1 − i2 − i3 = 0,
5i2 − 7i3 − 2i4 = 0.
Computational Methods (SECE) Chapter 5 Yonas Y. 3 / 50
Linear Systems: Direct Methods Introduction

Linear Systems: Problem Statement

We consider linear systems of equation of the form


n
X
aik xk = bi , i = 1, . . . , n.
k=1

or
Ax = b.

The matrix elements aik and the right-hand side elements bi are given.
We are looking for the unknowns xk .

Computational Methods (SECE) Chapter 5 Yonas Y. 4 / 50


Linear Systems: Direct Methods Introduction

Direct vs. iterative methods

For a linear system of n equations in n variables,

Ax = b

A is given, real, nonsingular matrix, n × n, and b is given, real vector.

Two types of solution approaches:


1 Direct methods: yield exact solution in absence of roundoff error.
Variations of Gaussian elimination.

2 Iterative methods: iterate in a similar fashion to what we do for nonlinear


problems.
Used only when direct methods are ineffective.

Computational Methods (SECE) Chapter 5 Yonas Y. 5 / 50


Gaussian elimination and backward substitution Backward substitution

Backward substitution
Occasionally, A has a special structure
If A is diagonal  
a11
 a22 
A=
 
.. 
 . 
ann
then the linear equations are uncoupled, and the solution is
bi
xi = , i = 1, 2, . . . , n.
aii
If A is an upper triangular matrix
 
a11 a12 · · · a1n
. .. 
a22 . .

 . 
A= .
,
 .. .. 
 . 
ann
Where aij = 0 for all i > j.
Computational Methods (SECE) Chapter 5 Yonas Y. 6 / 50
Gaussian elimination and backward substitution Backward substitution

Solve an upper triangular system of equations backwards.


bn
Last row reads ann xn = bn so xn = ann
bn−1 −an−1,n xn
Next xn−1 = an−1,n−1

Then find xn−2 , etc.

Algorithm: Backward Substitution.


Given an upper triangular matrix A and a right-hand-side b,
1 for k = n : −1 : 1
Pn
bk − j=k+1 akj xj
2 xk = akk
3 end

Note that the successful completion of the algorithm is possible only if the
diagonal elements of A, aii , are nonzero.

Computational Methods (SECE) Chapter 5 Yonas Y. 7 / 50


Gaussian elimination and backward substitution Backward substitution

Cost of backward substitution

What is the cost of this algorithm?

In a simplistic way we just count each floating point operation (such as +


and ∗) as a flop.

The number of flops required here is


n−1
X n
X
1+ ((n − k) + (n − k) + 1) = (2(n − k) + 1) = n2 .
k=1 k=1

Thus the cost of the algorithm is represented as O(n2 ).

In addition, data locality and data movement are crucial to the execution of
an algorithm.

Computational Methods (SECE) Chapter 5 Yonas Y. 8 / 50


Gaussian elimination and backward substitution Forward substitution

Forward substitution

If A is a lower triangular matrix


 
a11
a21 a22 
A= . ,
 
 .. .. ..
. . 
an1 an2 ··· ann

Where aij = 0 for all i < j.

Algorithm: Forward Substitution.


Given an lower triangular matrix A and a right-hand-side b,
1 for k = 1 : n
Pk−1
bk − j=1 akj xj
2 xk = akk
3 end
Computational Methods (SECE) Chapter 5 Yonas Y. 9 / 50
Gaussian elimination and backward substitution Forward substitution

Example 5.1. Consider solving the equations

5x1 = 15,
x1 + 2x2 = 7,
−x1 + 3x2 + 2x3 = 5.

The matrix A is lower triangular, and in our general notation we have


   
5 0 0 15
A= 1 2 0 , b =  7 .
−1 3 2 5

Applying the forward substitution algorithm we get


15
x1 = 5 = 3,
7−3
then x2 = 2 = 2,
5+3−6
then x3 = 2 = 1.

Computational Methods (SECE) Chapter 5 Yonas Y. 10 / 50


Gaussian elimination and backward substitution Gaussian elimination

Gaussian elimination
Elementary row transformations are used to reduce the given linear system to an
equivalent upper triangular form.
The reason this works is that the solution to our problem Ax = b is invariant to
multiplication of a row by a constant,
subtraction of a multiple of one row from another, and
row interchanges.

Performing each of the operations can be recorded directly on the augmented


matrix  
a11 a12 · · · a1n b1
 a21 a22 . . . a2n b2 
 
(A|b) =  .

.. .. ..  .

 .. ..
. . . . 
an1 an2 · · · ann bn

Computational Methods (SECE) Chapter 5 Yonas Y. 11 / 50


Gaussian elimination and backward substitution Gaussian elimination

Eliminating one column at a time

x1 x2 x3 x4
a11 a12 a13 a14 b1
a21 a22 a23 a24 b2
a31 a32 a33 a34 b3
a41 a42 a43 a44 b4

1 Permute rows i = 1, . . . , 4 (if necessary) such that a11 6= 0. This element is


called the pivot.
2 Subtract multiples li1 = ai1 /a11 of row 1 from row i, i = 2, . . . , 4.
0
3 Set aik = aik − li1 a1k , i, k = 2, . . . , 4.
0
4 Set bi = bi − li1 b1 , i = 2, ..., 4.

Computational Methods (SECE) Chapter 5 Yonas Y. 12 / 50


Gaussian elimination and backward substitution Gaussian elimination

x1 x2 x3 x4
a11 a12 a13 a14 b1
0 0 0 0
0 a22 a23 a24 b2
0 0 0 0
0 a32 a33 a34 b3
0 0 0 0
0 a42 a43 a44 b4

0
1 Permute rows i = 2, . . . , 4 (if necessary) such that a22 6= 0. This element is
the next pivot.
0 0 0
2 Subtract multiples li2 = ai2 /a22 of row 2 from row i, i = 3, . . . , 4.
00 0 0 0
3 Set aik = aik − li2 a2k , i, k = 3, . . . , 4.
00 0 0 0
4 Set bi = bi − li2 b2 , i = 3, ..., 4.

Computational Methods (SECE) Chapter 5 Yonas Y. 13 / 50


Gaussian elimination and backward substitution Gaussian elimination

x1 x2 x3 x4
a11 a12 a13 a14 b1
0 0 0 0
0 a22 a23 a24 b2
00 00 00
0 0 a33 a34 b3
00 00 00
0 0 a43 a44 b4

00
1 Permute rows i = 3, . . . , 4 (if necessary) such that a33 6= 0. This is the next
pivot.
00 00 00
2 Subtract multiples l43 = a43 /a33 of row 3 from row 4.
000 00 00 00
3 Set a44 = a44 − l43 a34 .
000 00 00 00
4 Set b4 = b4 − l43 b3 .

Computational Methods (SECE) Chapter 5 Yonas Y. 14 / 50


Gaussian elimination and backward substitution Gaussian elimination

Example 5.2. Consider the simple example

x1 − 4x2 = −10
0.5x1 − x2 = −2

We have n = 2; hence the loops in the general algorithm collapse into


defining just one update for k = 1, i = 2, j = 2:
 
1 −4 −10
(A|b) = ,
0.5 −1 −2
a21
l21 = = 0.5,
a11
0
a22 = −1 − 0.5 · (−4) = 1,
0
b2 = −2 − 0.5 · (−10) = 3,
 
0 0 1 −4 −10
(A |b ) = ,
0 1 3

Computational Methods (SECE) Chapter 5 Yonas Y. 15 / 50


Gaussian elimination and backward substitution Gaussian elimination

Algorithm: Gaussian Elimination.


Given a real, nonsingular n × n matrix A and and a vector b of size n, first
transform into upper triangular form,
1 for k = 1 : n − 1
2 for i = k + 1 : n
3 lik = aakkik
4 for j = k + 1 : n
5 aij = aij − lik akj
6 end
7 bi = bi − lik bk
8 end
9 end
Next, apply the algorithm of backward substitution.

Computational Methods (SECE) Chapter 5 Yonas Y. 16 / 50


Gaussian elimination and backward substitution Gaussian elimination

Example 5.3. Use Gaussian elimination with backward substitution to solve the
following linear system.

4x1 − x2 + x3 = 8,
2x1 + 5x2 + 2x3 = 3,
x1 + 2x2 + 4x3 = 11.

The augmented matrix will be


 
4 −1 1 8
(A|b) =  2 5 2 3 ,
1 2 4 11

We have n = 3; hence the loops in the general algorithm collapse into defining
after two update’s for k = 1, 2, i = 2, 3, j = 2, 3:
For k = 1,
a21 a31
l21 = = 0.5, l31 = = 0.25.
a11 a11
0 0 0
a22 = 5−0.5·(−1) = 5.5, a23 = 2−0.5·(1) = 1.5, b2 = 3−0.5·(8) = −1.
0 0 0
a32 = 2−0.25·(−1) = 2.25, a33 = 4−0.25·(1) = 3.75 b3 = 11−0.25·(8) = 9.
Computational Methods (SECE) Chapter 5 Yonas Y. 17 / 50
Gaussian elimination and backward substitution Gaussian elimination

Thus
 
0 0
4 −1 1 8
(A |b ) =  0 5.5 1.5 −1  ,
0 2.25 3.75 9
For k = 2,
0
0 a 2.25
l32 = 32
0 = = 0.41.
a22 5.5
00 00
a33 = 3.75 − 0.41 · (1.5) = 3.135, b3 = 9 − 0.41 · (−1) = 9.41.

Thus  
00 00
4 −1 1 8
(A |b ) =  0 5.5 1.5 −1  ,
0 0 3.135 9.41

Then apply backward substitution


00
b3 9.41 −1−1.5×3 8−1×3+1×−1
x3 = 00
a33
= 3.135 = 3, x2 = 5.5 = −1, x1 = 4 = 1.

Computational Methods (SECE) Chapter 5 Yonas Y. 18 / 50


Gaussian elimination and backward substitution LU decomposition

LU decomposition
Gaussian elimination decomposes matrix A into a product L × U of a lower
triangular matrix L and an upper triangular matrix U.

Elementary lower triangular matrices


First, elementary row operations are applied to matrix A to zero out its first
column.  
1
−l21 1 
 
(1) −l31 1
M = .

 .. . .. 
 . 
−ln1 1

This matrix has zeros everywhere except in the main diagonal and the first
column.
Then,
A(1) = M (1) A and b(1) = M (1) b

Computational Methods (SECE) Chapter 5 Yonas Y. 19 / 50


Gaussian elimination and backward substitution LU decomposition

Likewise, to zero out the second column of A(1) below the main diagonal

A(2) = M (2) A(1) = M (2) M (1) A and b(1) = M (2) b(1) = M (2) M (1) b

where  
1

 1 

..
M (2)
 
=
 −l32 . .

 .. .. 
 . . 
−ln2 1

After n − 1 stages, matrix A is transformed into an upper triangular matrix,

U = A(n−1) = M (n−1) · · · M (2) M (1) A

and likewise
b(n−1) = M (n−1) · · · M (2) M (1) b

Computational Methods (SECE) Chapter 5 Yonas Y. 20 / 50


Gaussian elimination and backward substitution LU decomposition

(n−1) −1
Multiplying U by [M ] , then by [M (n−2) ]−1 , etc., we obtain
A = LU
where L = [M (1) ]−1 · · · [M (n−2) ]−1 [M (n−1) ]−1 .

Thus, the elements of L will take the folowing form


 
1
l21 1 
 
L = l31 l32 1 .
 
 .. .. .. .. 
. . . . 
ln1 ln2 ··· ln,n−1 1

Figure: 5.1 LU decomposition for the case n = 10.


Computational Methods (SECE) Chapter 5 Yonas Y. 21 / 50
Gaussian elimination and backward substitution LU decomposition

Example 5.4. Decompose the following matrix


 
1 −1 3
A = 1 1 0 .
3 −2 1

1
1. The Gaussian elimination process for the first column yields l21 = 1 = 1,
l31 = 13 = 3, so
   
1 0 0 1 −1 3
M (1) = −1 1 0 , A(1) = M (1) A = 0 2 −3 .
−3 0 1 0 1 −8

2. The Gaussian elimination process for the second column yields l32 = 12 , so
   
1 0 0 1 −1 3
M (2) = 0 1 0 , U = A(2) = M (2) A(1) = 0 2 −3  .
0 −0.5 1 0 0 −6.5

Computational Methods (SECE) Chapter 5 Yonas Y. 22 / 50


Gaussian elimination and backward substitution LU decomposition

3. Note that
   
1 0 0 1 0 0
[M (1) ]−1 = 1 1 0 , [M (2) ]−1 = 0 1 0 ,
3 0 1 0 0.5 1
 
1 0 0
[M (1) ]−1 [M (2) ]−1 = 1 1 0 = L,
3 0.5 1

and
    
1 0 0 1 −1 3 1 −1 3
LU = 1 1 0 0 2 −3  = 1 1 0 = A
3 0.5 1 0 0 −6.5 3 −2 1

Gaussian elimination procedure decomposes A into a product of a unit lower


triangular matrix L and an upper triangular matrix U.

Computational Methods (SECE) Chapter 5 Yonas Y. 23 / 50


Gaussian elimination and backward substitution LU decomposition

Algorithm: Solving Ax = b by LU Decomposition.

Given a real nonsingular matrix A, apply LU decomposition first:

A = LU.

Given also a right-hand-side vector b:

1 Forward substitution: solve


Ly = b.

2 Backward substitution: solve


Ux = y.

Computational Methods (SECE) Chapter 5 Yonas Y. 24 / 50


Gaussian elimination and backward substitution LU decomposition

Example 5.5. Determine the LU factorization for matrix A in the linear system
Ax = b, where
   
1 1 0 3 8
2 1 −1 1  7
A= 3 −1 −1 2 
 and b=  14  .

−1 2 3 −1 −7
and solve the system.

The LU decomposition of matrix A is:


    
1 1 0 3 1 0 0 0 1 1 0 3
2 1 −1 1  2 1 0 0 0 −1 −1 −5 
A=  3 −1 −1 2  =  3
    = LU.
4 1 0 0 0 3 13 
−1 2 3 −1 −1 3 0 1 0 0 0 −13

To solve
     
1 0 0 0 1 1 0 3 x1 8
2 1 0 0 0 −1 −1 −5  x2   7 
Ax = LUx = 
3
   =  .
4 1 0 0 0 3 13  x3   14 
−1 3 0 1 0 0 0 −13 x4 −7
Computational Methods (SECE) Chapter 5 Yonas Y. 25 / 50
Gaussian elimination and backward substitution LU decomposition

We first introduce the substitution y = Ux. Then b = L(Ux) = Ly. That is,
    
1 0 0 0 y1 8
 2 1 0 0 y2   7 
Ax = LUx =   3 4 1 0 y3  =  14  .
   

−1 3 0 1 y4 −7

This system is solved for y by a simple forward-substitution process:

y1 = 8, y2 = −9, y3 = 26 and y4 = −26.

We then solve Ux = y for x, the solution of the original system; that is,
    
1 1 0 3 x1 8
0
 −1 −1 −5   x2  =  −9  .
   
0 0 3 13  x3   26 
0 0 0 −13 x4 −26

Using backward substitution we obtain x4 = 2, x3 = 0, x2 = 1, x1 = 3.


Computational Methods (SECE) Chapter 5 Yonas Y. 26 / 50
Pivoting strategies

Pivoting strategies
The process of Gaussian elimination (or LU decomposition) cannot be applied
unmodified to all nonsingular matrices.
For instance
    
0 1 x1 4
Ax = b ⇐⇒ =
1 1 x2 7

x1 x2
0 1 4
1 1 7

Pivot a11 = 0 −→ we exchange rows 1 and 2.


x1 x2
1 1 7
0 1 4
It is rare to hit a precisely zero pivot, but common to hit a very small one.
Undesirable accumulation of roundoff error may arise in such cases.
Computational Methods (SECE) Chapter 5 Yonas Y. 27 / 50
Pivoting strategies

Example 5.6. Consider the system


x1 + x2 + x3 = 1,
x1 + x2 + 2x3 = 2,
x1 + 2x2 + 2x3 = 1.

Matrix A is nonsingular and has a unique solution of x = (1, −1, 1)T .


But applying Gaussian elimination yields after the first stage
   
1 1 1 1 1 1 1 1
 1 1 2 2  ⇒  0 0 1 1 
1 2 2 1 0 1 1 0
0
Next, a22 = 0 and we cannot proceed further.
One obvious solution is:
   
1 1 1 1 1 1 1 1
 1 2 2 1  ⇒  0 1 1 0 
1 1 2 2 0 0 1 1

and backward substitution proceeds as usual.


Computational Methods (SECE) Chapter 5 Yonas Y. 28 / 50
Pivoting strategies

Example 5.7. Consider the numerical example (5 decimal digits)

0.00035x1 + 1.2654x2 = 3.5267


1.2547x1 + 1.3182x2 = 6.8541

x1 x2 x1 x2
0.00035 1.2654 3.5267 0.00035 1.2654 3.5267
1.2547 1.3182 6.8541 0 -4535.0 -12636

with
l21 = 1.2547/0.00035 ≈ 3584.9
0
a22 = 1.3182 − 3584.9 × 1.2654 ≈ 1.3182 − 4536.3 ≈ −4535.0
0
b2 = 6.8541 − 3584.9 × (3.5267) ≈ 6.8541 − 12643 ≈ −12636.

Computational Methods (SECE) Chapter 5 Yonas Y. 29 / 50


Pivoting strategies

Backsubstitution gives

x2 = −12636/(−4535.0) ≈ 2.7863
x1 = (3.5267 − 1.2654 × 2.7863)/0.00035
≈ (3.5267 − 3.5258)/0.00035 = 0.0009/0.00035 = 2.5714.

The exact solution (up to 5 decimal digits) is

x1 = 2.5354, x2 = 2.7863

We have cancellation in the computation of x1 !


0
The very small pivot in the first elimination step induced large numbers a22
0
and b2 . Information, e.g. a22 , got lost.

Computational Methods (SECE) Chapter 5 Yonas Y. 30 / 50


Pivoting strategies Partial pivoting

Partial pivoting

(k−1)
If a pivotal element akk is small in magnitude
Roundoff errors are amplified, both in subsequent stages of the Gaussian
elimination process and (even more apparently) in the backward substitution
phase.

The most common strategy for pivoting is to search for the maximal element in
modulus:
The index p of the pivot in the k-th step of GE is determined by
(k−1) (k−1)
|apk | = max |aik |.
i≥k

If p > k then rows p and k are exchanged.


This strategy implies that |lik | ≤ 1. Then proceed with the elimination
process.

Computational Methods (SECE) Chapter 5 Yonas Y. 31 / 50


Pivoting strategies Partial pivoting

Let’s consider example 5.7 using Gaussian elimination with partial pivoting(GEPP)

x1 x2 x1 x2
1.2547 1.3182 6.8541 1.2547 1.3182 6.8541
0.00035 1.2654 3.5267 0 1.2650 3.5248

with
l21 = 0.00027895
x2 = 2.7864
x1 = (6.8541 − 1.3182 × 2.7864)/1.2547
≈ (6.8541 − 3.6730)/1.2547 = 3.1811/1.2547 ≈ 2.5353.

There is a deviation from the exact solution in the last digit only.

Computational Methods (SECE) Chapter 5 Yonas Y. 32 / 50


Pivoting strategies Partial pivoting

How does partial pivoting affect the LU decomposition?

Note that for the matrix of


 
1 1 1
A = 1 1 2 .
1 2 2
Is it possible to write A = LU with L unit lower triangular and U nonsingular
upper triangular matrix?
NO!

However, the operation of row interchange is captured in general by a


permutation matrix P.
For the above matrix
PA = LU,
where  
1 0 0
P = 0 0 1 .
0 1 0
Computational Methods (SECE) Chapter 5 Yonas Y. 33 / 50
Pivoting strategies GEPP stability

GEPP stability

Question: Does the incorporation of a partial pivoting procedure guarantee


stability of the Gaussian elimination algorithm, in the sense that roundoff errors do
not get amplified by factors that grow unboundedly as the matrix size n increases?

We have seen that |lik | ≤ 1. What about |uik |?

We would like to be sure that


(k)
gn (A) = max |ai,j |,
i,j,k

i.e., the elements in U do not grow to fast (exponentially as a function of n)


in the course of GE.

Unfortunately, it is not possible in general to guarantee stability of GEPP.

Computational Methods (SECE) Chapter 5 Yonas Y. 34 / 50


Pivoting strategies GEPP stability

Complete pivoting

In complete pivoting we look for (one of) the largest element in modulus
(k−1) (k−1)
|aqr | = max |aij |,
k≤i,j≤n

at each stage of GE and interchange both row q and column r with row i and
column j, respectively.

This strategy does guarantee stability.

However, searching in each stage makes it significantly more expensive than


partial pivoting.

Instances in which (scaled) partial pivoting really fails appear to be extremely


rare in practice ⇒ GEPP is the commonly used approach.

Computational Methods (SECE) Chapter 5 Yonas Y. 35 / 50


Pivoting strategies Matrices requiring no pivoting

Matrices requiring no pivoting

Families of matrices for which it is known that pivoting is not required.

1 Symmetric positive matrices.


2 Diagonally dominant matrices.
A matrix A ∈ Rn×n is diagonally dominant, if
X
|aii | ≥ |aik |, i = 1, . . . , n.
k6=i

Theorem
If A is nonsingular and diagonally dominant then the LU factorization can be
computed without pivoting.

Computational Methods (SECE) Chapter 5 Yonas Y. 36 / 50


The Cholesky decomposition

The Cholesky decomposition

Recall that a matrix A is symmetric if AT = A and positive definite if

xT Ax > 0, ∀x 6= 0.

If a symmetric matrix A ∈ Rn,n is positive definite, then the following conditions


hold
1 aii > 0 for i = 1, . . . , n.
2
2 aik < aii akk for i 6= k, i, k = 1, . . . , n.
3 There is a k with maxi,k |aik | = akk .

Gaussian elimination without pivoting is a stable algorithm for symmetric positive


definite matrices.

Computational Methods (SECE) Chapter 5 Yonas Y. 37 / 50


The Cholesky decomposition

Symmetrizing the LU decomposition


Consider the LU decomposition of a symmetric positive definite matrix, written as
A = LU.
We can write
u12 u1n 
··· ···
 
u11 1 u11 u11
u23 u2n 

 u22 
 1 u22 ··· u22 
 ..  .. .. 
U=
 . 
 . .  = D Ũ,

 ..  .. .. 
 .  . . 
unn 1
where ukk > 0.

So the LU decomposition reads


A = LD Ũ.
Transposing it
A = AT = Ũ T DLT ,
Computational Methods (SECE) Chapter 5 Yonas Y. 38 / 50
The Cholesky decomposition

Since this decomposition is unique we must have Ũ T = L, i.e., the decomposition


reads

A = LDLT ,

Note that L is a unit lower triangular matrix and D is a diagonal matrix with
positive diagonal entries.

It is also possible to write D = D 1/2 D 1/2 with


√ √
D 1/2 = diag { u11 , . . . , unn },

Then,
A = GG T ,

with G = LD 1/2 a lower triangular matrix. This is called the Cholesky


decomposition.

Computational Methods (SECE) Chapter 5 Yonas Y. 39 / 50


The Cholesky decomposition

Example 5.8. In the 2 × 2 case we can require


    
a a12 g11 0 g11 g21
A = 11 =
a21 a22 g21 g22 0 g22

We have three unknown quantities g11 , g12 , and g22 .


2 √
For a11 we get a11 = g11 =⇒ g11 = a11 .
a12
For a12 =⇒ g21 = g11 .
p
Finally, for a22 =⇒ g22 = 2 .
a22 − g21

This way of decomposing, can be straightforwardly applied to larger matrices in


the same fashion.
A key here, and the reason attempting to follow the same strategy for general
matrices (that is, not necessarily symmetric positive definite) will fail, is that the
arguments of the square root are guaranteed to be positive if A is symmetric
positive definite.

Computational Methods (SECE) Chapter 5 Yonas Y. 40 / 50


The Cholesky decomposition

Algorithm: Cholesky Decomposition.


Given a symmetric positive definite n × n matrix A, this algorithm overwrites its
lower part with its Cholesky factor.
1 for k = 1 : n − 1

2 akk = akk
3 for i = k + 1 : n
4 aik = aakkik
5 end
6 for j = k + 1 : n
7 for i = j : n
8 aij = aij − aik ajk
9 end
10 end
11 end

12 ann = ann

Computational Methods (SECE) Chapter 5 Yonas Y. 41 / 50


Estimating errors and the condition number

Estimating errors and the condition number


The error in the numerical solution
Suppose that, using some algorithm, we have computed an approximate solution
x.
b
Absolute error =⇒ k x−bxk
kx−b
xk
Relative error =⇒ kxk

Two questions regarding the accuracy of b x as an approximaton to the solution x


of the linear system of equations Ax = b.
1 First we investigate what we can derive from the size of the residual
br = b − Abx.
Note that ideally, r = b − Ax = 0, but in practice there will always be some
nonzero residual r because of roundoff errors.
2 Then, what is the effect of errors in the initial data (b, A) on the solution x?
That is, how sensitive is the solution to perturbations in the initial data?
Computational Methods (SECE) Chapter 5 Yonas Y. 42 / 50
Estimating errors and the condition number

Example 5.9. Let


   
1.2969 0.8648 0.8642
A= , b= .
0.21612 0.1441 0.1440

Suppose that, using some algorithm, the approximate solution is


 
0.9911
x= .
−0.4870
b

Then
−10−8
 
br = b − Ab
x= =⇒ k br k∞ = 10−8 .
10−8

However  
2
x= , =⇒ k x−b
x k∞ = 1.513.
−2

Thus, the error is of the same size as the solution itself, roughly 108 times that of
the residual!

Computational Methods (SECE) Chapter 5 Yonas Y. 43 / 50


Estimating errors and the condition number

Condition number and a relative error estimate


So, how does the residual br = b − Ab x − x?
x affect the error z := b

We assume that k A k is any matrix norm and k x k a compatible vector


norm, i.e., we have k Ax k≤k A kk x k for all k x k.
We have
x − x) = Ab
Az = A(b x − b = −br.

Therefore

k b k=k Ax k ≤ k A kk x k, k z k=k −A−1br k ≤ k A−1 kk br k

We get an estimate for the relative error of x:

kzk kb
x−xk k br k k br k
= ≤ k A kk A−1 k =: κ(A)
kxk kxk kbk kbk

Computational Methods (SECE) Chapter 5 Yonas Y. 44 / 50


Estimating errors and the condition number

Definition
We therefore define the condition number of the matrix A as
κ(A) =k A kk A−1 k
is called the condition number of A.
The relative error in the solution is bounded by the condition number of the
matrix A times the relative error in the residual.

The condition number is at least 1 (ideally conditioned):


1 =k I k=k AA−1 k ≤ k A kk A−1 k= κ(A).

For a singular matrix at the other end of the scale, κ(A) = ∞.


If A is nonsingular and E is the matrix with smallest norm such that A + E is
singular, then
k E k2 1
=
k A k2 κ2 (A)
So, the condition number measures how close to singularity in the relative sense a
matrix is.
Computational Methods (SECE) Chapter 5 Yonas Y. 45 / 50
Estimating errors and the condition number

Example 5.10. Continuing with Example 5.9, we have


 
−1 8 0.1441 −0.8648
A = 10 .
−0.2161 1.2969
This yields
k A−1 k∞ = 1.513 × 108 , hence κ∞ (A) = 2.1617 · 1.513 × 108 ≈ 3.27 × 108 .
1.513 3.27
The numbers 2 < 0.8642 confirm the estimate
k z k∞ k br k∞
≤ κ∞ (A) .
k x k∞ k b k∞
a12 a21
We now make A singular by replacing a22 = a11

Then A + E is singular with


   
0 0 0 0
E= ≈ ,
0 0.8648×0.2161
1.2969 − 0.144 0 −7.7 × 10−9
So, indeed,
kE k 1
≈ .
kAk κ(A)
(This estimate holds in `2 and `∞ norms.)
Computational Methods (SECE) Chapter 5 Yonas Y. 46 / 50
Estimating errors and the condition number Sensitivity on matrix coefficients

Sensitivity on matrix coefficients

Input data A and b are often perturbed due, e.g., to rounding.

How big is the change δx if the matrix A is altered by δA and b is perturbed


by δb.

The LU decomposition can be considered as the exact decomposition of a


slightly perturbed matrix à = A + δA.

Forward/backward substitition add (rounding) errors.

All together: The computed solution b


x is the solution of a perturbed problem.

(A + δA)b
x = b + δb, x = x + δx.
b

How do these perturbations affect the solution?


This is called backward error analysis.

Computational Methods (SECE) Chapter 5 Yonas Y. 47 / 50


Estimating errors and the condition number Sensitivity on matrix coefficients

A simple substitution shows that

br = b − Ab x − δb.
x = (δA)b

How large are δA and δb, and thus br?

An LU decomposition procedure with pivoting can be shown to yield a


perturbation bounded by

k δA k∞ ≤ ηφ(n)gn (A),

where
φ(n) is a low order polynomial in n (cubic at most) and
η is the rounding unit.

The extra perturbations which the forward and backward substitutions yield,
in particular the bound on k δb k∞ , are significantly smaller.

Computational Methods (SECE) Chapter 5 Yonas Y. 48 / 50


Estimating errors and the condition number Sensitivity on matrix coefficients

Plugging the expression for br into the relative error expression

k x−b
xk k br k
≤ κ(A)
kxk kbk

we finally obtain

k x−b
xk k δb k + k δA kk b
xk
≤ κ(A)
kxk kbk

This can further be expressed as, provided k δA k < 1/ k A−1 k


 
k x−b
xk κ(A) k δb k k δA k
≤ + .
kxk 1 − κ(A)(k δA k / k A k) kbk kAk

In summary, a stable algorithm is responsible for producing a small residual.

Computational Methods (SECE) Chapter 5 Yonas Y. 49 / 50


Estimating errors and the condition number Sensitivity on matrix coefficients

End of Chapter 5

Questions?

Computational Methods (SECE) Chapter 5 Yonas Y. 50 / 50

You might also like