Transcript An Introduction to the Quality of Computed Solutions

An Introduction to the Quality
of Computed Solutions
Sven Hammarling
NAG Ltd, Oxford
[email protected]
1
Plan of Talk
• Introduction
• Floating point numbers and IEEE arithmetic
• Condition, stability and error analysis with
examples
• Implications for software.
– LAPACK and NAG
• Other approaches
• Summary
2
Introduction to NAG
3
NAG History
• 1970 • 1971 • 1973 -
•
•
•
•
•
1976 1978 1980 1990 1998 -
Nottingham Algorithms Group
Mark 1 NAG Library
NAG moved to Oxford, renamed
Numerical Algorithms Group
NAG Ltd, a non-profit company
NAG Inc established in USA
NAG Ltd financially self-sufficient
NAG GmbH established
Nihon NAG KK established
4
What does non-profit mean?
No shareholders, no owner
surplus is re-invested in the Company
5
How many people work at
NAG?
 85 worldwide,  65 in Oxford
 12 in Downers Grove
Distributors worldwide
6
NAG Products and Services
• Product Lines
–
–
–
–
–
Numerical Libraries
Statistical Systems
NAGWare: Compiler and tools
Visualization and Graphics
PDE Solutions
• Consultancy
• Customer Support
www.nag.co.uk or www.nag.com
7
NAG Numerical Libraries
F77
FL90Plus
F90
C
Parallel
SMP
MPI
8
Quality of Computed Solutions
9
Quality of Computed Solutions
The quality of computed solutions is concerned
with assessing how good a computed solution is
in some appropriate measure
10
Software Quality
Quality software should implement reliable
algorithms and should provide measures of
solution quality
11
Floating Point Numbers
Floating point numbers are a subset of the
real numbers that can be conveniently
represented in the finite word length of a
computer, without unduly restricting the
range of numbers represented
For example, the IEEE standard uses 64 bits
to represent double precision numbers in the
approximate range
12
Floating Point Numbers - Representation
13
Floating Point Numbers - Example
14
Floating Point Numbers - Example
(cont’d)
0
1
2
3
15
IEEE Arithmetic Standard
ANSI/IEEE Standard 754-1985 is a standard for
binary arithmetic (b = 2). The standard specifies:
•
•
•
•
Floating point number formats
Results of the basic floating point operations
Rounding modes
Signed zero, infinity (
) and not-a-number
(NaN)
• Floating point exceptions and their handling
• Conversion between formats
Most modern machines use IEEE arithmetic.
16
IEEE Arithmetic Formats
Format
Single
Double
Extended
Pr ecision
24 bits
53 bits
 64
Exponent
8 bits
11 bits
 15
Approx Range
1038
10308
104932
Approx precision
108
1016
1020
17
Why Worry about Computed
Solutions?
• Tacoma bridge collapse
• North Sea oil rig collapse
• Vancouver stock exchange index, Jan 1982
to Nov 1983
• Ariane 5 rocket, flight 501, failure, 4 June
1996
• Patriot missile, 25 Feb 1991
• Auckland Bridge
• London millennium bridge
18
Tacoma Bridge
19
Auckland Bridge, 1975
20
Millennium Bridge, 2000
21
Web Sites
• Disasters attributable to bad numerical
computing:
http://www.math.psu.edu/dna/disasters/
• Numerical problems: RISKS-LIST:
http://catless.ncl.ac.uk/Risks/
• London millennium bridge:
http://www.arup.com/MillenniumBridge/
22
Example - Means
Mean of a set of numbers can be outside
their range. Using 3 figure arithmetic:
 5.01  5.03
2  10.0 2  5.00
Does this matter? If so, can we guarantee
an answer within range?
23
Example - Sample Variance
n
1
2
2
sn 
 xi  x 

n  1 i 1
1
2
n
n

1
1
 
2
2
sn 
2
  xi    xi  

n  1  i 1
n  i 1  
For xT  10000,10001,10002  using 8 figure arithmetic
1 gives:
1.0,
 2  gives: 0.0
Are either of these reasonable? (Actually (1) is the
correct answer)
24
Excel Example:
Standard Deviation
25
Overflow/Underflow Example:
Modulus of a Complex Number
x  xr2  xi2
2
b
 a 1    , where a  max  xr , xi  , b  min  xr , xi
a
 x  0 if a  0 
26

Excel Example:
Standard Deviation - Overflow
27
Stability
The stability of a method for solving a
problem is concerned with the sensitivity of
the method to (rounding) errors in the
solution process
A method that guarantees as accurate a
solution as the data warrants is said to be
stable, otherwise the method is unstable
28
Condition
The condition of a problem is concerned with
the sensitivity of the problem to perturbations
in the data
A problem is ill-conditioned if small changes
in the data cause relatively large changes in
the solution. Otherwise a problem is wellconditioned
29
Condition Examples
x 3  21x 2  120 x  100  0
x1  1, x2  x3  10
0.99 x 3  21x 2  120 x  99.99  0
x1  1, x2  1117
. , x3  9.04
101
. x  21x  120 x  100.01  0
x1  1, x2 , x3  9.90  104
. i
3
2
30
Condition Examples (cont’d)
I
z
10
10
x
(ae  be )dx
x
When a  b  1,
I 0
When a  1, b  101
.
I  220
33
Condition Examples (cont’d)
99 x1  98 x2  197
100 x1  99 x2  199
x1  x2  1
98.99 x1  98 x2  197
100 x1  99 x2  199
x1  100, x2  99
35
Condition Number for Linear
Equations
Ax  b, ( A  E ) x  b
A( x  x)  Ex, so that x  x  A1Ex
giving
xx
x
 A
1
 E   ( A)
E
A
,
where  ( A)  A  A1 is the condition number of
A with respect to matrix inversion
37
Condition Number Example
 99 98 
A
 , A 1  199
100 99 
 99 98 
1
1
A 
 , A 1  199
 100 99 
 1  A  199
2
4 10
4
38
Stability Examples
1.6 x 2  100.1x  1.251  0
Four significant figure arithmetic on the standard formula
x  (b  b 2  4ac) (2a)
gives
x1  62.53,
x2  0.03125,
but using
x2  c (ax1 )
gives
x2  0.01251
Exact solution is
x1  62.55,
x2  0.0125
39
Stability Examples (Contd)
1
yn  (1 e)  x n e x dx
0
Easy to show that
0  yn  1 ( n  1)
Integrating by parts gives
yn  1  nyn 1 , y0  1  1 e  0.632121
recursive.m
40
Stability Examples (Contd)
~
y0  y0  
~
~
y1  1  y0  y1  
~
~
y2  1  2 y1  y2  2
~
y3  1  3 ~
y2  y3  6
In general
n
~
y n  y n  1 n ! 
bg
41
Stability Examples (cont’d)
yn 1  1  yn  n
yn  m  yn  m  
yn  m 1  1  yn  m   n  m 
 yn  m 1    n  m 
yn  m  2  yn  m  2  
  n  m  n  m  1 
In general
yn  yn   1
brecursive.m
mn

  n  m  n  m  1  n  1 
42
Error Analysis
Error analysis is concerned with establishing
whether or not an algorithm is stable for the
problem in hand
A forward error analysis is concerned with how
close the computed solution is to the exact
solution
A backward error analysis is concerned with
how well the computed solution satisfies the
problem to be solved
43
Example
Ax  b, r  b  Ax
 99 98 
 1
1
A
, b   , x   
100 99 
 1
 1
 2.97 
Consider x  

 2.99 
 1.97 
 0.01
xx 
 , but r  


1.99
0.01




 1.01 
Now consider x  
,
 0.99 
 0.01
 1.97 
xx 
 , but r  

0.
01

1.97




The Purpose of Error Analysis
“The clear identification of the factors determining the
stability of an algorithm soon led to the development of
better algorithms. The proper understanding of inverse
iteration for eigenvectors and the development of the
QR algorithm by Francis are the crowning
achievements of this line of research.
“For me, then, the primary purpose of the rounding
error analysis was insight.”
Wilkinson, 1986
Bulletin of the IMA, Vol 22, p197
45
46
Backward Error and
Perturbation Analysis
Ax  b,
bA  E g~x  b
If we know how perturbations in A affect the
solution x , then we can estimate the accuracy of the
computed solution ~
x . That is, an estimate of the
backward error allows us to estimate the forward error.
x~
x
E
 A
~
x
A
bg
47
Condition and Error Analysis
forward error 
condition number  backward error
LAPACK
Linear Algebra PACKage for high-performance computers
• Systems of linear equations
• Linear least squares problems
• Eigenvalue and singular value problems,
including generalized problems
• Matrix factorizations
• Condition and error estimates
• The BLAS as a portability layer
Dense and banded linear algebra for Shared Memory
49
LAPACK and Error Bounds
“In Addition to providing faster routines than
previously available, LAPACK provides more
comprehensive and better error bounds. Our goal is to
provide error bounds for most quantities computed by
LAPACK.”
LAPACK Users’ Guide, Chapter 4 – Accuracy and
Stability
Error Bounds in LAPACK
The Users’ Guide gives details of error
bounds and code fragments to
compute those bounds. In many cases
the routines return the bounds directly
Error Bounds in LAPACK:
Example
DGESVX is an 'expert' driver for solving
AX  B
DGESVX( ..., RCOND, FERR, BERR, WORK,...)
RCOND
- Estimate of 1/   A 
FERR(j) - Estimated forward error for X j
BERR(j) - Componentwise relative backward error for X j
(smallest relative change in any element of A and B
that makes X j an exact solution)
WORK(1) - Reciprocal of pivot growth factor (1/  )
ODE Software Example
NAG routine D02PZF returns global error estimates
for the Runge-Kutta solvers D02PCF and D02PDF:
D02PCZ( RMSERR, ERRMAX, TERRMX,…)
RMSERR: Approximate Root mean square error
ERRMAX: Max approximate true error
TERRMX: Point at which max approximate true
error occurred
(These routines are based upon RKSUITE)
54
What to do if Software does not
Provide Error Estimates?
• Run the problem with perturbed data
• Better still, use a software tool such as
PRECISE which allows you to perform a
stochastic analysis
• Put pressure on developers to provide
estimates
PRECISE
• Provides a module for statistical backward
error analysis
• Provides a module for sensitivity analysis
• http://www.cerfacs.fr/algor/Softs/PRECISE/
index.html
Quality of Reliable Software
Chaitin-Chatelin & Frayssé define the quality index
of a reliable algorithm at the computed solution x as
J ( x) 
B( x)

,
where B( x) is the backward error and  is the machine
precision.
(Lectures on Finite Precision Computations)
CADNA
• CADNA is another package with the same
aims as PRECISE
http://www-anp.lip6.fr/cadna/Accueil.php
• It is now free to academics
58
Cancellation
Using four figure decimal arithmetic
s  1.000  1.000 10  1.000 10
4
 1.000 10  1.000 10
0
4
4
4
59
Interval Analysis
• Interval arithmetic works on intervals that contain
the exact solution
1.000
1.000  1.000 104 1.000 104   1.000 104 1.000 104 
 1.000 104 1.0010 104   1.000 104 1.000 104 
 0 10
• Interval analysis can provide automatic error
bounds for a number of problems
• See www.cs.utep.edu/interval-comp/main.html
60
Summary
• Error analysis gives insight
• In linear algebra we usually aim for a
backward error analysis
• Be concerned about the quality of
computed solutions
• Use quality software, complain about low
quality software
• Write quality software and be proud of
your results
61
Quote
“You have been solving these damn
problems better than I can pose them.”
Sir Edward Bullard, Director NPL, in a
remark to Wilkinson. (mid 1950s.)
See NAG Newsletter 2/83, page 11
62
References
[1] E. Anderson, Z. Bai, C. H. Bischof, J. Demmel, J. J. Dongarra,
J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S.
Blackford & D. C. Sorensen. LAPACK Users’ Guide, 3rd Edition,
SIAM, 1999.
[2] F. Chaitin-Chatelin & V. Frayssé. Lectures on Finite Precision
Computations, SIAM, 1996.
[3] G. H. Golub & C. F. Van Loan. Matrix Computations, 3rd
Edition, The Johns Hopkins University Press, 1996.
[4] N. J. Higham. Accuracy and Stability of Numerical Algorithms,
2nd edition, SIAM, 2002.
[5] www.cs.utep.edu/interval-comp/main.html
[6] www.nag.co.uk/numeric/numerical_libraries.asp
References (cont’d)
[7] F. S. Acton. Numerical Methods that Usually Work. Harper and
Row, 1970.
[8] G. E. Forsythe. What is a Satisfactory Quadratic Equation
Solver. In Constructive Aspects of the Fundamental Theorem of
Algebra, B. Dejon and P. Henrici, P., editors, pages 53-61, Wiley,
1969.
[9] G. E. Forsythe. Pitfalls in Computation, or Why a Math Book
isn't enough. Amer. Math. Monthly, 9: 931-995, 1970.
[10] J. H. Wilkinson. The Perfidious Polynomial. In Studies in
Numerical Analysis, G. H. Golub, editor, volume 24 of Studies in
Mathematics, Mathematical Association of America, Washington
D.C., pages 1-24, 1984.