Transcript Engineering Computation Part 3 E. T. S. I. Caminos, Canales y Puertos Learning Objectives for Lecture 1.

Engineering
Computation
Part 3
E. T. S. I. Caminos, Canales y Puertos
1
Learning Objectives for Lecture
1. Motivate Study of Systems of Equations and particularly
Systems of Linear Equations
2. Review steps of Gaussian Elimination
3. Examine how roundoff error can enter and
be magnified in Gaussian Elimination
4. Introduce Pivoting and Scaling as defenses against
roundoff.
5. Consider what an engineer can do to generate well
formulated problems.
E. T. S. I. Caminos, Canales y Puertos
2
Systems of Equations
•
In Part 2 we have tried to determine the value x, satisfying f(x)=0. In
this part we try to obtain the values x1,x2, xn, satisfying the system of
equations:
•
These systems can be linear or nonlinear, but in this part we deal with
linear systems:
E. T. S. I. Caminos, Canales y Puertos
3
Systems of Equations
•
•
where a and b are constant coefficients, and n is the number of
equations.
Many of the engineering fundamental equations are based on
conservation laws. In mathematical terms, these principles lead to
balance or continuity equations relating the system behavior with
respect to the amount of the magnitude being modelled and the
extrenal stimuli acting on the system.
E. T. S. I. Caminos, Canales y Puertos
4
Systems of Equations
•
Matrices are rectangular sets of elements represented by a single
symbol. If the set if horizontal it is called row, and if it is vertical, it is
called column.
Column 3
Row 2
Column
vector
Row vector
E. T. S. I. Caminos, Canales y Puertos
5
Systems of Equations
•
There are some special types of matrices:
Symmetric matrix
Diagonal matrix
E. T. S. I. Caminos, Canales y Puertos
Identity matrix
Upper triangular matrix
6
Systems of Equations
Lower triangular matrix
Banded matrix
Half band width
All elements are null with the
exception of thoise in a band
centered around the main diagonal.
This matrix has a band width of 3
and has the name of tridiagonal.
E. T. S. I. Caminos, Canales y Puertos
7
Systems of Equations
Linear Algebraic Equations
a11x1 + a12x2 + a13x3 + … + a1nxn = b1
a21x1 + a22x2 + a23x3 + … + a2nxn = b2
…..
an1x1 + an2x2 + an3x3 + … + anxn = bn
where all aij's and bi's are constants.
In matrix form:  a11
 a 21
 a 31

 a n1
a12
a 22
a 32
a13
a 23
a 33
a n2
a n3
a1n 
a 2n 
a 3n 

a nn 
nxn
or simply [A]{x} = {b}
E. T. S. I. Caminos, Canales y Puertos
 x1 
 b1 


x2 

b2 

=
x
b
 3
 3
 
 


x n 

b n 

nx1
nx1
8
Systems of Equations
Matrix product:
Resulting dimensions
•
Matrix representation of a system
E. T. S. I. Caminos, Canales y Puertos
9
Systems of Equations
•
Graphic Solution: Systems of equations are hyperplanes (straight
lines, planes, etc.). The solution of a system is the intersection of
these hyperplanes.
Compatible and determined
system. Vectors are linearly
independent. Unique solution.
Determinant of A is non-null.
E. T. S. I. Caminos, Canales y Puertos
10
Systems of Equations
Incompatible system, Linearly
dependent vectors. Null determinant
of A. There is no solution.
Compatible but undetermined
system. Linearly dependent vectors.
Null determinant of A. There exists
an infinite number of solutions.
E. T. S. I. Caminos, Canales y Puertos
11
Systems of Equations
Compatible and determined system.
Linearly independent vectors. Nonnull
determinant of A, but close to zero. There
exists a solution but it is difficult to find
precisely. It is an ill conditioned system
leading to numerical errors.
E. T. S. I. Caminos, Canales y Puertos
12
Gauss elimination
•
Naive Gauss elimination method: The Gauss’ method has two phases:
Forward elimination and backsustitution. In the first, the system is reduced to
an upper triangular system:
Pivot
equation
substract
pivot
•
First, the unknown x1 is eliminated. To this end, the first row is multiplied by a21/a11 and added to the second row. The same is done with all other
succesive rows (n-1 times) until only the first equation contains the first
unknown x1.
E. T. S. I. Caminos, Canales y Puertos
13
Gauss elimination
•
This operation is repeated with all variables xi, until an upper triangular matrix
is obtained.
•
Next, the system is solved by backsustitution.
•
The number of operations (FLOPS) used in the Gauss method is:
Pass 1
E. T. S. I. Caminos, Canales y Puertos
Pass 2
14
Gauss elimination
1. Forward Elimination (Row Manipulation):
a. Form augmented matrix [A|b]:
 a11
a 21

 a n1
a12
a 22
a n2
 a11
a1n   x1   b1 
a 2n   x 2    b 2   a 21

   
 x  b 
a n1
a nn  
n  n
a12
a 22
a n2
a1n b1 
a 2n b 2 

a nn b n 
b. By elementary row manipulations, reduce [A|b] to
[U|b'] where U is an upper triangular matrix:
DO i = 1 to n-1
DO k = i+1 to n
Row(k) = Row(k) - (aki/aii)*Row(i)
ENDDO
ENDDO
E. T. S. I. Caminos, Canales y Puertos
15
Gauss elimination
2. Back Substitution
Solve the upper triangular system [U]{x} = {b´}
 u11 u12
 0 u
22

0
 0


 0
0
u13
u 23
u 33
0
u1n   x1   b1 
u 2n   x 2   b2 
    
u 3n   x 3    b3 
   
   
u nn   x n  bn 
xn = b'n / unn
DO i = n-1 to 1 by (-1)
n
xi 
bi   u ij x j
ji 1
u ii
END
E. T. S. I. Caminos, Canales y Puertos
16
Gauss elimination (example)
Consider the system of equations
50
1
 2
1
2
40 4 
6 30 



 x1 
1 

x
=
2
 2
 
3 



 x3 

To 2 significant figures, the exact solution is:
0.016 


x true    0.041
 0.091


We will use 2 decimal digit arithmetic with rounding.
E. T. S. I. Caminos, Canales y Puertos
17
Gauss elimination (example)
Start with the augmented matrix:
50
1
2

1
2
40 4
6 30
1
2
3
Multiply the first row by –1/50
and add to second row.
Multiply the first row by –2/50
and add to third row:
50
0
0

1
2
40 4
6 30
1
2
3
50
0
0

1
40
0
Multiply the second row by –6/40
and add to third row:
E. T. S. I. Caminos, Canales y Puertos
2
4
29
1
2
2.7 
18
Gauss elimination (example)
Now backsolve:
50
0
0

x3 
x2
x1
1
40
0
2.7
 0.093
29
2  4x 3 


 0.040
40
1  2x 3  x 2 


 0.016
50
E. T. S. I. Caminos, Canales y Puertos
2
4
29
1
2
2.7 
(vs. 0.091, et = 2.2%)
(vs. 0.041, et = 2.5%)
(vs. 0.016, et = 0%)
19
Gauss elimination (example)
Consider an alternative solution
interchanging rows:
After forward elimination, we obtain:
2
50
1

2

0
 0
6 30
1
2
40 4
6
150
0
3
1
2 
30
3 
750 74 

200 19 
Now backsolve:
x3 = 0.095
x2 = 0.020
x1 = 0.000
(vs. 0.091, et = 4.4%)
(vs. 0.041, et = 50%)
(vs. 0.016, et = 100%)
Apparently, the order of the equations matters!
E. T. S. I. Caminos, Canales y Puertos
20
Gauss elimination (example)
WHAT HAPPENED?
• When we used 50 x1 + 1 x2 + 2 x3 = 1 to solve for x1, there
was little change in other equations.
• When we used 2 x1 + 6 x2 + 30 x3 = 3 to solve for x1 it made
BIG changes in the other equations. Some coefficients for
other equations were lost!
The second equation has little to do with x1.
It has mainly to do with x3.
As a result we obtained LARGE numbers in the table,
significant roundoff error occurred and information was lost.
Things didn't go well!
• If scaling factors | aji / aii | are  1 then the effect of roundoff
errors is diminished.
E. T. S. I. Caminos, Canales y Puertos
21
Gauss elimination (example)
Effect of diagonal dominance:
As a first approximation roots are:
xi  bi / aii
Consider the previous examples:
50
1
2

2
50
1

1
2
40 4
6 30
6 30
1
2
40 4
E. T. S. I. Caminos, Canales y Puertos
1
2
3
3
1
2 
x true 

0.016 

=  0.041

 0.091

x1  1/50 =0.02
x 2  2/40 =0.05
x 3  3/30 =0.10
x1  3/2 =1.5
x 2  1/1 =1.0
x 3  2/4 = 0.50
22
Gauss elimination (example)
Goals: 1. Best accuracy (i.e. minimize error)
2. Parsimony (i.e. minimize effort)
Possible Problems:
A. Zero on diagonal term  ÷ by zero.
B. Many floating point operations (flops) cause numerical
precision problems and propagation of errors.
C. System may be ill-conditioned: det[A]  0.
D. No solution or an infinite # of solutions: det[A] = 0.
Possible Remedies:
A. Carry more significant figures (double precision).
B. Pivot when the diagonal is close to zero.
C. Scale to reduce round-off error.
E. T. S. I. Caminos, Canales y Puertos
23
Gauss elimination (pivoting)
PIVOTING
A. Row pivoting (Partial Pivoting) In any good routine, at each step i, find
maxk | aki | for k = i, i+1, i+2, ..., n
Move corresponding row to pivot position.
(i) Avoids zero aii
(ii) Keeps numbers small & minimizes round-off,
(iii) Uses an equation with large | aki | to find xi
 Maintains diagonal dominance.
 Row pivoting does not affect the order of the variables.
 Included in any good Gaussian Elimination routine.
E. T. S. I. Caminos, Canales y Puertos
24
Gauss elimination (pivoting)
B. Column pivoting Reorder remaining variables xj
for j = i, . . . ,n so get largest | aji |
Column pivoting changes the order of the
unknowns, xi, and thus leads to complexity in the
algorithm. Not usually done.
C. Complete or Full pivoting
Performing both row pivoting and column pivoting.
(If [A] is symmetric, needed to preserve symmetry.)
E. T. S. I. Caminos, Canales y Puertos
25
Gauss elimination (pivoting)
How to fool pivoting:
Multiply the third equation by 100 and then performing
pivoting will yield:
 200 600 3000 300

 50
 1
1
40
2
4

1 

2 
Forward elimination then yields (2-digit arithmetic):
 200

 0
 0
600
150
0
3000
750
200
Backsolution yields:
x3 = 0.095 (vs. 0.091, et = 4.4%)
x2 = 0.020 (vs. 0.041, et = 50.0%)
x1 = 0.000 (vs. 0.016, et = 100%)
The order of the rows is still poor!!
E. T. S. I. Caminos, Canales y Puertos
300 
74 

19 
26
Gauss elimination (scaling)
SCALING
A. Express all equations (and variables) in comparable units so all
elements of [A] are about the same size.
B. If that fails, and maxj |aij| varies widely across the rows, replace each
row i by:
a ij
aij 
max j | a ij |
This makes the largest coefficient |aij| of each equation equal to 1 and
the largest element of [A] equal to 1 or -1
NOTE: Routines generally do not scale automatically; scaling can cause
round-off error too!
SOLUTIONS
• Don't actually scale, but use hypothetical scaling factors to determine
what pivoting is necessary.
• Scale only by powers of 2: no roundoff or division required.
E. T. S. I. Caminos, Canales y Puertos
27
Gauss elimination (scaling)
How to fool scaling:
A poor choice of units can undermine
the value of scaling.
Begin with our original example:
If the units of x1 were expressed
in µg instead of mg the matrix
might read:
50

1
 2
1
2 1
40 4 2 

6 30 3 
50000

 1000
 2000
1
40
6
2 1
4 2

3 3 
Scaling then yields:
1 0.00002 0.00001 0.00001


1
0.04
0.004
0.002


1 0.003
0.015
0.0015 
Which equation is used to determine x1 ? Why bother to scale ?
E. T. S. I. Caminos, Canales y Puertos
28
Gauss elimination (operation counting)
OPERATION COUNTING
In numerical scientific calculations, the number of multiplies &
divides often determines CPU time. (This represents the
numerical effort!)
One floating point multiply or divide (plus any associated
adds or subtracts) is called a FLOP. (The adds/subtracts use
little time compared to the multiplies/divides.)
FLOP = FLoating point OPeration.
Examples:
E. T. S. I. Caminos, Canales y Puertos
a *x + b
a / x – b
29
Gauss elimination (operation counting)
Useful identities in counting FLOPS:
m
1)
m
 c f (i)  c f (i)
i 1
m
m
i 1
m
i  1 2  3  4 
i 1
m
6)

i 1
i 2  12  22 
1  m
i 1
[f (i)  g(i)]   f (i)   g(i)
i 1
5)
1  1  1 
3)
i 1
m
2)
m
i 1
m
4)
1  m  k  1
ik
m(m  1) m 2
m 

 O(m)
2
2
3
m(m

1)(2m

1)
m
 m2 

 O(m 2 )
6
3
O(mn) means that there are terms of order mn and lower.
E. T. S. I. Caminos, Canales y Puertos
30
Gauss elimination (operation counting)
Simple Example of Operation Counting:
DO i = 1 to n
Y(i) = X(i)/i – 1
ENDDO
X(i) and Y(i) are arrays whose values change when i
changes. In each iteration
X(i)/i – 1
represents one FLOP because it requires one division (&
one subtraction).
The DO loop extends over i from 1 to n iterations:
n
1  n FLOPS
E. T. S. I. Caminos, Canales y Puertos
i 1
31
Gauss elimination (operation counting)
Another Example of Operation Counting:
DO i = 1 to n
Y(i) = X(i) X(i) + 1
DO j = i to n
Z(j) = [ Y(j) / X(i) ] Y(j) + X(i)
ENDDO
ENDDO
With nested loops, always start from the innermost loop.
[Y(j)/X(i)] * Y(j) + X(i) represents 2 FLOPS
n
n
 2  21  2(n  i  1)
ji
FLOPS
ji
E. T. S. I. Caminos, Canales y Puertos
32
Gauss elimination (operation counting)
For the outer i-loop: X(i) • X(i) + 1 represents 1 FLOP
n
n
n
i 1
i 1
i 1
 1  2  n  i  1)   3  2n  1  2 i
  3  2n  n  2
n(n  1)
2
= 3n +2n2 - n2 - n
= n2 + 2n = n2 + O(n)
E. T. S. I. Caminos, Canales y Puertos
33
Gauss elimination (operation counting)
Forward Elimination:
DO k = 1 to n–1
DO i = k+1 to n
r = A(i,k)/A(k,k)
DO j = k+1 to n
A(i,j)=A(i,j) – r*A(k,j)
ENDDO
B(i) = B(i) – r*B(k)
ENDDO
ENDDO
E. T. S. I. Caminos, Canales y Puertos
34
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Back Substitution:
X(n) = B(n)/A(n,n)
DO i = n–1 to 1 by –1
SUM = 0
DO j = i+1 to n
SUM = SUM + A(i,j)*X(j)
ENDDO
X(i) = [B(i) – SUM]/A(i,i)
ENDDO
E. T. S. I. Caminos, Canales y Puertos
35
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Forward Elimination
n
Inner loop:
 1 = n - (k +1) +1 = n - k
j=k+1
n
Second loop:
 2 + (n - k)  (2  n)  k  (n  k)
i=k+1
= (n2 + 2n) – 2(n + 1)k + k2
E. T. S. I. Caminos, Canales y Puertos
36
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Forward Elimination (cont'd)
n 1
Outer loop
=
2
2
[(
n

2
n
)

2(
n

1)
k

k
]

k 1
n 1
 (n 2  2n)

k 1
n 1
1 2(n  1)
n 1
 
k
k 1
 (n 2  2n)(n  1)  2(n  1)
k2
k 1
(n  1))(n)
2
(n  1)(n)(2n  1) n3

=
+ O(n2 )
6
3
E. T. S. I. Caminos, Canales y Puertos
37
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Back Substitution
n
Inner Loop:
 1  n  (i  1)  1  n - i
ji 1
Outer Loop:
n 1
n 1
n 1
i 1
i 1
i 1
1  (n  i)  (1  n)1   i
(n  1)n
 (1  n)(n  1) 
2
n2
=
+ O(n)
2
E. T. S. I. Caminos, Canales y Puertos
38
Gauss elimination (operation counting)
Total flops = Forward Elimination + Back Substitution
=
n3/3 + O (n2)
+ n2/2 + O (n)

n3/3 + O (n2)
To convert (A,b) to (U,b') requires n3/3, plus terms of order n2 and smaller,
flops.
To back solve requires:
1 + 2 + 3 + 4 + . . . + n = n (n+1) / 2 flops;
Grand Total: the entire effort requires n3/3 + O(n2) flops altogether.
E. T. S. I. Caminos, Canales y Puertos
39
Gauss-Jordan Elimination
Diagonalization by both forward and backward elimination in
each column.  a11 a12 a13 ... a1n   x1   b1 
a
 21
 a 31


a n1
a 22
a 32
a 23
a 33
a n2
a n3
 b 
... a 2n  
x
2
2



   
... a 3n   x 3  =  b3 
   
...
   
... a nn  
 x nn 
 
b n 

Perform elimination both backwards and forwards until:
1
0

0


0
0 0
1
0
0 1
0 0
0
0 
0


1 
 x1 
 x1 
x 
x 
2
2


 

 

=
x
x
 3
 3
 
 
 
 


x n 

x n 

3
Operation count for Gauss-Jordan is: n  O(n 2 )
(slower than Gauss elimination)
2
E. T. S. I. Caminos, Canales y Puertos
40
Gauss-Jordan Elimination
Example (two-digit arithmetic):
50
1
2

1
2
40 4
6 30
1
2
3
1 0 0.038
0 1
0.1
0 0
29

x1 = 0.015
x2 = 0.041
x3 = 0.093
E. T. S. I. Caminos, Canales y Puertos
1 0.02 0.04
0 40
4
0
6
30

0.019 
0.05 
2.7 
1 0 0
0 1 0
0 0 1

0.02 
2
3
0.015
0.041
0.093
(vs. 0.016, et = 6.3%)
(vs. 0.041, et = 0%)
(vs. 0.091, et = 2.2%)
41
Gauss-Jordan Matrix Inversion
The solution of: [A]{x} = {b} is:
{x} = [A]-1{b}
where [A]-1 is the inverse matrix of [A]
Consider:
[A] [A]-1 = [ I ]
1) Create the augmented matrix: [ A | I ]
2) Apply Gauss-Jordan elimination:
==> [ I | A-1 ]
E. T. S. I. Caminos, Canales y Puertos
42
Gauss-Jordan Matrix Inversion
Gauss-Jordan Matrix Inversion (with 2 digit arithmetic):
50
 A I  =  1
2

1 0 0  1 0.02 0.04
4
0 1 0   0 40
6
30
0 0 1  0
1
2
40 4
6 30
0.02 0 0 
-0.02 1 0 
-0.04 0 1 
1 0 0.038
0.02
0.005 0 


0
1
0.1

0.0005
0.025
0


0 0
28
0.037
0.15 1 
1 0 0
0.02

0 1 0 0.00037
0 0 1 0.0013
0.00029
0.026
0.0054
0.0014 
0.0036 

0.036 
MATRIX INVERSE [A-1]
E. T. S. I. Caminos, Canales y Puertos
43
Gauss-Jordan Matrix Inversion
CHECK:
[ A ] [ A ]-1 = [ I ]
50
2
 2
1
40
6
2
4
30 

 0.020
-0.00037
 -0.0013
-0.00029
0.026
-0.0054
-0.0014   0.997
-0.0036    0.000
0.036 
  0.001
0.13
1.016
0.012
0.002 
0.001
1.056 

[ A ]-1 { b } = { x }
 0.020
-0.00037
 -0.0013
-0.0029
0.026
-0.0054
-0.0014 
-0.0036 
0.036 

E. T. S. I. Caminos, Canales y Puertos
1  0.015 
 2    0.033 
 3  0.099 
x true
0.016 
  0.041
 0.091
 0.016 
x   0.040 
 0.093 
Gaussian
Elimination
44
LU decomposition
•
LU decomposition - The LU decomposition is a method that uses the
elimination techniques to transform the matrix A in a product of triangular
matrices. This is specially useful to solve systems with different vectors b,
because the same decomposition of matrix A can be used to evaluate in an
efficient form, by forward and backward sustitution, all cases.
A = L U
E. T. S. I. Caminos, Canales y Puertos
45
LU decomposition
A = L U
Decomposition
U X = D
A X = B
Initial system
L U X = B
Transformed system 1
Substitution
L D = B
Transformed system 2
X
Backward sustitution
D
Forward sustitution
E. T. S. I. Caminos, Canales y Puertos
46
LU decomposition
•
LU decomposition is very much related to Gauss method, because the upper
triangular matrix is also looked for in the LU decomposition. Thus, only the
lower triangular matrix is needed.
•
Surprisingly, during the Gauss elimination procedure, this matrix L is obtained,
but one is not aware of this fact. The factors we use to get zeroes below the
main diagonal are the elements of this matrix L.
Substract
E. T. S. I. Caminos, Canales y Puertos
47
LU decomposition
resto
resto
E. T. S. I. Caminos, Canales y Puertos
48
LU decomposition (Complexity)
Basic Approach
Consider [A]{x} = {b}
a) Gauss-type "decomposition" of [A] into [L][U]
n3/3 flops
[A]{x} = {b} becomes [L] [U]{x} = {b}; let [U]{x}  {d}
b) First solve [L] {d} = {b} for {d} by forward subst.
n2/2 flops
c) Then solve [U]{x} = {d} for {x} by back substitution
n2/2 flops
E. T. S. I. Caminos, Canales y Puertos
49
LU Decompostion: notation
a11
 L =  a ij

0 
a nn 
a11
 U =  0

1
L
=
 1  a ij

0
1 
1 a ij 
U
=
 1  0 1 


0

 L0  = a ij

0
0
0 a ij 

 U0  = 0 0 


1 0

D
=
  0 1 
E. T. S. I. Caminos, Canales y Puertos
a ij 
a nn 
[A] = [L] + [U0]
[A] = [L0] + [U]
[A] = [L0] + [U0] + [D]
[A] = [L1] [U]
[A] = [L] [U1]
50
LU decomposition
LU Decomposition Variations
Doolittle [L1][U]
General [A]
Crout
[L][U1]
General [A]
Cholesky [L][L] T
Pos. Def. Symmetric [A]
Cholesky works only for Positive Definite Symmetric matrices
Doolittle versus Crout:
• Doolittle just stores Gaussian elimination factors where Crout
uses a different series of calculations (see C&C 10.1.4).
• Both decompose [A] into [L] and [U] in n3/3 FLOPS
• Different location of diagonal of 1's
• Crout uses each element of [A] only once so the same array
can be used for [A] and [L\U] saving computer memory!
E. T. S. I. Caminos, Canales y Puertos
51
LU decomposition
Matrix Inversion
Definition of a matrix inverse:
[A] [A]-1 = [ I ]
==> [A] {x} = {b}
[A]-1 {b} = {x}
First Rule:
Don’t do it.
(numerically unstable calculation)
E. T. S. I. Caminos, Canales y Puertos
52
LU decomposition
Matrix Inversion
If you really must -1) Gaussian elimination: [A | I ] –> [U | B'] ==> A-1
2) Gauss-Jordan: [A | I ] ==> [I | A-1 ]
Inversion will take
n3 + O(n2)
flops if one is careful about where zeros are (taking advantage of
the sparseness of the matrix)
Naive applications (without optimization) take 4n3/3 + O(n2) flops.
For example, LU decomposition requires n3/3 + O(n2) flops. Back
solving twice with n unit vectors ei:
2 n (n2/2) = n3 flops.
Altogether: n3/3 + n3 = 4n3/3 + O(n2) flops
E. T. S. I. Caminos, Canales y Puertos
53
FLOP Counts for Linear Algebraic Equations
Summary
FLOP Counts for Linear Algebraic Equations, [A]{x} =
{b}
Gaussian Elimination (1 r.h.s)
n3/3 + O (n2)
Gauss-Jordan (1 r.h.s)
n3/2 + O (n2)
LU decomposition
n3/3 + O (n2)
Each extra LU right-hand-side
n2
Cholesky decomposition (symmetric A) n3/6 + O (n2)
Inversion (naive Gauss-Jordan)
4n3/3 +O (n2)
Inversion (optimal Gauss-Jordan)
Solution by Cramer's Rule
E. T. S. I. Caminos, Canales y Puertos
n3 + O (n2)
n!
54
Errors in Solutions to Systems of Linear Equations
System of Equations
Errors in Solutions to Systems of Linear Equations
Objective:
Solve [A]{x} = {b}
Problem:
Round-off errors may accumulate and even be
exaggerated by the solution procedure. Errors are
often exaggerated if the system is ill-conditioned
Possible remedies to minimize this effect:
1. Partial or complete pivoting
2. Work in double precision
3. Transform the problem into an equivalent system of
linear equations by scaling or equilibrating
E. T. S. I. Caminos, Canales y Puertos
55
Errors in Solutions to Systems of Linear Equations
Ill-conditioning
• A system of equations is singular if det|A| = 0
• If a system of equations is nearly singular it is ill-conditioned.
Systems which are ill-conditioned are extremely sensitive to
small changes in coefficients of [A] and {b}. These systems
are inherently sensitive to round-off errors.
Question:
Can we develop a means for detecting these situations?
E. T. S. I. Caminos, Canales y Puertos
56
Errors in Solutions to Systems of Linear Equations
Ill-conditioning of [A]{x} = {b}:
Consider the graphical interpretation for a 2-equation system:
 a11
a
 21
a12   x1 
 b1 
  =  

a 22   x 2 
b 2 
(1)
(2)
We can plot the two linear equations on a graph of x1 vs. x2.
a11x1+ a12x2 = b1
x2
b2/a22
b1/a12
b1/a11
E. T. S. I. Caminos, Canales y Puertos
b2/a21
x1
a21x1+ a22x2 = b2
57
Errors in Solutions to Systems of Linear Equations
Ill-conditioning of [A]{x} = {b}:
Consider the graphical interpretation for a 2-equation system:
 a11
a
 21
a12   x1 
 b1 
  =  

a 22   x 2 
b 2 
(1)
(2)
We can plot the two linear equations on a graph of x1 vs. x2.
x1
x1
x2
Uncertainty
in x2
Well-conditioned
E. T. S. I. Caminos, Canales y Puertos
x2
Uncertainty
in x2
Ill-conditioned
58
Errors in Solutions to Systems of Linear Equations
Ways to detect ill-conditioning:
1. Calculate {x}, make small change in [A] or {b} and
determine change in solution {x}.
2. After forward elimination, examine diagonal of upper
triangular matrix. If aii << ajj, i.e. there is a relatively
small value on diagonal, then this may indicate illconditioning.
3. Compare {x}single with {x}double
4. Estimate "condition number" for A.
Substituting the calculated {x} into [A]{x} and checking this
against {b} will not always work!!!
E. T. S. I. Caminos, Canales y Puertos
59
Errors in Solutions to Systems of Linear Equations
Ways to detect ill-conditioning:
• If det|A| = 0 the matrix is singular
==> the determinant may be an indicator of conditioning
• If det|A| is near zero is the matrix ill-conditioned?
• Consider:
0
0 
0.01
det  0
0.01
0   106
0
0.01
 0
After scaling:
1 0 0 
det 0 1 0   1
0 0 1 
==> det|A| will provide an estimate of conditioning if it is
normalized by the "magnitude" of the matrix.
E. T. S. I. Caminos, Canales y Puertos
60
Norms
Norms and the Condition Number
We need a quantitative measure of ill-conditioning.
This measure will then directly reflect the possible magnitude
of round off effects.
To do this we need to understand norms:
Norm: Scalar measure of the magnitude of a matrix or vector
("how big" a vector is).
Not to be confused with the dimension of a matrix.
E. T. S. I. Caminos, Canales y Puertos
61
Vector Norms
Vector Norms: Scalar measure of the magnitude of a vector
Here are some vector norms for n x 1
vectors {x} with typical elements xi.
Each is in the general form of a p norm
defined by the general relationship:
x
p




1/ p
n
x
i 1
p
i




 n

x 1 
xi 
1. Sum of the magnitudes:


 i 1

x   max x i
2. Magnitude of largest element:
i
(infinity norm)

3. Length or Euclidean norm:
E. T. S. I. Caminos, Canales y Puertos
x
2




1/ 2
n
x
i 1
2
i




62
Norms
Vector Norms
Required Properties of vector norm:
1. ||x||  0 and ||x|| = 0 if and only if [x]=0
2 ||kx|| = k ||x|| where k is any positive scalar
3. ||x+y||  ||x|| + ||y||
Triangle Inequality
For the Euclidean vector norm we also have
4. ||x•y||  ||x|| ||y||
because the dot product or inner product property
satisfies:
||xy|| = ||x||•||y|| |cos()|  ||x|| • ||y||.
E. T. S. I. Caminos, Canales y Puertos
63
Matrix Norms
Matrix Norms: Scalar measure of the magnitude of a matrix.
Matrix norms corresponding to vector norms above are defined
by the general relationship:
A
p
 max
x 1
 Ax
p
n
1. Largest column sum:
(column sum norm)
2. Largest row sum:
(row sum norm)
(infinity norm)
E. T. S. I. Caminos, Canales y Puertos
A 1 max
j
a
ij
i 1
n
A

 max
i
a
ij
j1
64
Matrix norms
3. Spectral norm: ||A|| 2 = (µmax)1/2
where µmax is the largest eigenvalue of [A]T[A]
If [A] is symmetric, (µmax)1/2 = max , is the largest
eigenvalue of [A].
(Note: this is not the same as the Euclidean or Frobenius
norm, seldom used:
Ae
E. T. S. I. Caminos, Canales y Puertos
n
n
2
A
  ij
i 1 j 1
65
Matrix norms
Matrix Norms
For matrix norms to be useful we require that
0. || Ax ||  || A || ||x ||
General properties of any matrix norm:
1. || A ||  0 and || A || = 0 iff [A] = 0
2. || k A || = k || A || where k is any positive scalar
3. || A + B ||  || A || + || B ||
"Triangle Inequality"
4. || A B ||  || A || || B ||
Why are norms important?
 Norms permit us to express the accuracy of the solution {x} in
terms of || x ||
 Norms allow us to bound the magnitude of the product [ A ] {x}
and the associated errors.
E. T. S. I. Caminos, Canales y Puertos
66
Error Analysis
Forward and backward error analysis can estimate the effect of truncation
and roundoff errors on the precision of a result. The two approaches are
alternative views:
1. Forward (a priori) error analysis tries to trace the accumulation of
error through each process of the algorithm, comparing the
calculated and exact values at every stage.
2. Backward (a posteriori) error analysis views the final solution as the
exact solution to a perturbed problem. One can consider how
different the perturbed problem is from the original problem.
Here we use the condition number of a matrix [A] to specify the amount
by which relative errors in [A] and/or {b} due to input, truncation, and
rounding can be amplified by the linear system in the computation of {x}.
E. T. S. I. Caminos, Canales y Puertos
67
Error Analysis
Backward Error Analysis of [A]{x} = {b} for errors in {b}
Suppose the coefficients {b} are not precisely represented. What might
be the effect on the calculated value for {x + dx}?
Lemma: [A]{x} = {b} yields ||A|| ||x||  ||b||
x 
b
A

or
A
1

x
b
Now an error in {b} yields a corresponding error in {x}:
[A ]{x + dx} = {b + db}
[A]{x} + [A]{ dx} = {b} + {db}
Subtracting [A]{x} = {b} yields:
[A]{dx} = {db} ––> {dx} = [A]-1{db}
E. T. S. I. Caminos, Canales y Puertos
68
Error Analysis
Backward Error Analysis of [A]{x} = {b} for errors in {b}
Taking norms we have:
x  A 1
And using the lemma:
A
1

x
b
we then have :
x
 A
b
1
A
b

b
x
b
b
Define the condition number as
k = cond [A]  ||A-1|| ||A||  1
If k  1 or k is small, the system is well-conditioned
If k >> 1, system is ill conditioned.
1 = || I || = || A-1A ||  || A-1 || || A || = k = Cond(A)
E. T. S. I. Caminos, Canales y Puertos
69
Error Analysis
Backward Error Analysis of [A]{x} = {b} for errors in [A]
If the coefficients in [A] are not precisely represented, what
might be effect on the calculated value of {x+ dx}?
[A + dA ]{x + dx} = {b}
[A]{x} + [A]{ dx} + [dA]{x+dx} = {b}
Subtracting [A]{x} = {b} yields:
[A]{ dx} = – [dA]{x+dx}
or
{dx} = – [A]-1 [dA] {x+dx}
Taking norms and multiplying by || A || / || A || yields :
x
x  x
 A
E. T. S. I. Caminos, Canales y Puertos
A
1
A
A

A
A
70
E. T. S. I. Caminos, Canales y Puertos
71
Linfield & Penny, 1999
E. T. S. I. Caminos, Canales y Puertos
72
Linfield & Penny, 1999
Error Analysis
Estimate of Loss of Significance:
Consider the possible impact of errors [dA] on the precision of {x}.
If
A
A
~ 10
p
,
x
implies that if
x  x
then
~ 10s ,
Or, taking log of both sides:
x
x  x

A
A
then 10s   10 p
s > p - log10()
• log10() is the loss in decimal precision; i.e., we start with p
decimal figures and end-up with s decimal figures.
• It is not always necessary to find [A]-1 to estimate k = cond[A].
Instead, use an estimate based upon iteration of inverse matrix
using LU decomposition.
E. T. S. I. Caminos, Canales y Puertos
73
Iterative Solution Methods
Impetus for Iterative Schemes:
1. May be more rapid if coefficient matrix is "sparse"
2. May be more economical with respect to memory
3. May also be applied to solve nonlinear systems
Disadvantages:
1. May not converge or may converge slowly
2. Not appropriate for all systems
Error bounds apply to solutions obtained by direct and iterative
methods because they address the specification of [dA] and
{db}.
E. T. S. I. Caminos, Canales y Puertos
74
Iterative Solution Methods
Basic Mechanics:
Starting with:
a11x1 + a12x2 + a13x3 + ... + a1nxn
a21x1 + a22x2 + a23x3 + ... + a2nxn
a31x1 + a32x2 + a33x3 + ... + a3nxn
:
an1x1 + an2x2 + an3x3 + ... + annxn
=
=
=
=
b1
b2
b3
:
bn
Solve each equation for one variable:
x1 = [b1 – (a12x2 + a13x3 + ... + a1nxn )} / a11
x2 = [b2 – (a21x1 + a23x3 + ... + a2nxn )} / a22
x3 = [b3 – (a31x1 + a32x2 + ... + a3nxn )} / a33
:
xn = [bn – (an1x2 + an2x3 + ... + an,n-1xn-1 )} / ann
E. T. S. I. Caminos, Canales y Puertos
75
Iterative Solution Methods
 Start with initial estimate of {x}0.
 Substitute into the right-hand side of all the equations.
 Generate new approximation {x}1.
 This is a multivariate one-point iteration: {x}j+1 = {g({x}j)}
 Repeat process until the maximum number of iterations
reached or until:
|| xj+1 – xj ||  d + e || xj+1 ||
E. T. S. I. Caminos, Canales y Puertos
76
Convergence
To solve
[A]{x} = {b}
Separate [A] into:
[A] = [Lo] + [D] + [Uo]
[D]
= diagonal (aii)
[Lo]
= lower triangular with 0's on diagonal
[Uo] = upper triangular with 0's on diagonal
Rewrite system:
[A]{x} = ( [Lo] + [D] + [Uo] ){x} = {b}
[D]{x} + ( [Lo] + [Uo] ){x} = {b}
Iterate:
[D]{x}j+1 = {b} – ( [Lo]+[Uo] ) {x}j
{x}j+1 = [D]-1{b} – [D]-1 ( [Lo]+[Uo] ) {x}j
Iterations converge if:
|| [D]-1 ( [Lo] + [Uo] ) || < 1
(sufficient if equations are diagonally dominant)
E. T. S. I. Caminos, Canales y Puertos
77
Iterative Solution Methods – the Jacobi Method


x1j1   b1  a12 x 2j  a13 x 3j 

x 2j1   b 2  a 21 x1j  a 23 x 3j 


x nj1   b n  a n1 x1j  a n2 x 2j 

E. T. S. I. Caminos, Canales y Puertos


 a1n x nj  a11

 a 2n x nj  a 22


 a n,n 1 x nj 1  a nn

78
Iterative Solution Methods -- Gauss-Seidel
In most cases using the newest values within the right-hand
side equations will provide better estimates of the next value. If
this is done, then we are using the Gauss-Seidel Method:
( [Lo]+[D] ){x}j+1 = {b} – [Uo] {x}j
or explicitly:

 a
x 2j1   b 2  a 21 x1j1  a 23 x 3j 

x 3j1   b3

31
x1j1  a 32 x 2j1 

x nj1   b n  a n1 x1j1  a n2 x 2j1 


 a 2n x nj  a 22


 a 3n x nj  a 33


 a n,n 1 x nj11  a nn

If this is done, then we are using the Gauss-Seidel Method
E. T. S. I. Caminos, Canales y Puertos
79
Iterative Solution Methods -- Gauss-Seidel
If either method is going to converge,
Gauss-Seidel will converge
faster than Jacobi.
Why use Jacobi at all?
Because you can separate the n-equations into n
independent tasks, it is very well suited computers with
parallel processors.
E. T. S. I. Caminos, Canales y Puertos
80
Convergence of Iterative Solution Methods
Rewrite given system:
[A]{x} = { [B] + [E] } {x} = {b}
where [B] is diagonal, or triangular so we can solve
[B]{y} = {g} quickly. Thus,
[B] {x}j+1= {b}– [E] {x}j
which is effectively:
{x}j+1 = [B]-1 ({b} – [E] {x}j )
True solution {x}c satisfies: {x}c = [B]-1 ({b} –
Subtracting yields:
So
[E] {x}c)
{x}c – {x}j+1= – [B]-1 [E] [{x}c – {x}j]
||{x}c – {x}j+1 ||  ||[B]-1 [E]|| ||{x}c – {x}j ||
Iterations converge linearly if || [B]-1 [E] || < 1
=> || ([D] + [Lo])-1 [Uo] || < 1 For Gauss-Seidel
=> || [D] -1 ([Lo] + [Uo]) || < 1 For Jacobi
E. T. S. I. Caminos, Canales y Puertos
81
Convergence of Iterative Solution Methods
Iterative methods will not converge for all systems of
equations, nor for all possible rearrangements.
If the system is diagonally dominant,
i.e., | aii | > | aij | where i  j then
bi a i1
a i2
xi 

x1 
x2 
a ii a ii
a ii
with all
a ij
a ii
E. T. S. I. Caminos, Canales y Puertos
a in

xn
a ii
< 1.0, i.e., small slopes.
82
Convergence of Iterative Solution Methods
A sufficient condition for convergence exists:
n
a ii 
a
ij
j1
ji
Notes:
1. If the above does not hold, still may converge.
2. This looks similar to infinity norm of [A]
E. T. S. I. Caminos, Canales y Puertos
83
Improving Rate of Convergence of G-S Iteration
Relaxation Schemes:
trial
old
xnew
=

x
+
(1
)
x
i
i
i
where 0.0 <   2.0
(Usually the value of l is close to 1)
Underrelaxation ( 0.0 <  < 1.0 )
More weight is placed on the previous value. Often used to:
- make non-convergent system convergent or
- to expedite convergence by damping out oscillations.
Overrelaxation ( 1.0 <   2.0 )
More weight is placed on the new value. Assumes that the
new value is heading in right direction, and hence pushes
new value close to true solution.
The choice of  is highly problem-dependent and is empirical,
so relaxation is usually only used for often repeated calculations
of a particular class.
E. T. S. I. Caminos, Canales y Puertos
84
Why Iterative Solutions?
We often need to solve [A]{x} = {b} where n = 1000's
• Description of a building or airframe,
• Finite-Difference approximations to PDE's. Most of A's
elements will be zero; a finite-difference approximation to
Laplace's equation will have five aij0 in each row of A.
Direct method (Gaussian elimination)
• Requires n3/3 flops (say n = 5000; n3/3 = 4 x 1010 flops)
• Fills in many of n2-5n zero elements of A
Iterative methods (Jacobi or Gauss-Seidel)
• Never store [A] (say n = 5000; [A] would need 4n2 = 100 Mb)
• Only need to compute [A-B] {x}; and to solve [B]{xt+1} = {b}
E. T. S. I. Caminos, Canales y Puertos
85
Why Iterative Solutions?
• Effort:
Suppose [B] is diagonal,
solving
[B] {v} = {b}
Computing
[A-B] x
n flops
4n flops
For m iterations
5mn flops
For n = m = 5000,
5mn = 1.25x108
At worst O(n2).
E. T. S. I. Caminos, Canales y Puertos
86