Introduction
Sensitivity and conditioning
Computer arithmetic
Bi7740: Scientific computing
Introductory considerations
Vlad Popovici
popovici@iba.muni.cz
Institute of Biostatistics and Analyses
Masaryk University, Brno
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
There is nothing more practical than a good theory.
Kurt Lewin (1890–1947)
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Outline
1 Introduction
2 Sensitivity and conditioning
3 Computer arithmetic
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Bibliography:
HEATH M.T. (2002). Scientific Computing. An introductory
survey. McGraw-Hill, 2nd edition. ISBN: 0-07-239910-4 Good
accompanying materials at
http://www.cs.illinois.edu/~heath/scicomp/,
including slides and demos! Used as basis for the first part of
the course.
KEPNER J. (2009). Parallel Matlab for Multicore and
Multinode Computers. SIAM Publishing. ISBN:
978-0-898716-73-3
GENTLE J.E. (2005). Elements of Computational Statistics.
Springer. ISBN:978-0387954899
H ˇREBÍ ˇCEK, J. et al. Vˇedecké výpoˇcty v matematické biologii
(Scientific computing in mathematical biology). Brno:
Akademické nakladatelství CERM, 2012. 117 pp. Neuveden.
ISBN 978-80-7204-781-9.
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Computing environments for the course:
Matlab, http://www.mathworks.com - commercial
GNU Octave, https://www.gnu.org/software/octave/ "quite
similar to Matlab"
R, http://www.r-project.org - "environment for statistical
computing and graphics"
WARNING: Some pieces of code shown during the course may
not represent the optimal implementation in the given language.
They are merely a device for demonstrating some principles.
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Scientific computing
Wikipedia:
"Computational science (also scientific computing or scientific
computation) is concerned with constructing mathematical models
and quantitative analysis techniques and using computers to
analyze and solve scientific problems."
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Scientific computing
Wikipedia:
"Computational science (also scientific computing or scientific
computation) is concerned with constructing mathematical models
and quantitative analysis techniques and using computers to
analyze and solve scientific problems."
Basically: find numerical solutions to mathematically-formulated
problems.
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
(J. Hadamard) A problem is well posed if its solution
exists
is unique
has a behavior that changes continuously with the initial
conditions;
otheriwse, it is ill posed.
Inverse problems are often ill posed.
Example: 3D to 2D projection.
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
continuous domain → discrete domain
well-posed but ill-conditioned problems: small errors in input
lead to large variations in the solution
improve conditioning by regularization
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
General computational approach
continuous domain → discrete domain
infinite → finite
differential → algebraic
nonlinear → (combination of) linear
accept approximate solutions, but control for the error
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Approximations
Modeling approximations:
"model" = approximation of the nature
data - inexact measurements or previous results
Implementation/computational approximations:
discretization of the continuous domain; truncation
rounding
errors in input data
errors propagated by the algorithm
accuracy of the final result
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Example: area of the Earth
model: sphere
A = 4πr2
r =?
π = 3.14159 . . .
rounded arithmetic
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Errors
Absolute error: approximate value - true value
Relative error:
absolute error
true value
→ approximate value = (1 + relative error) × (true value)
if the relative error is ∼ 10−d
, it means that ˆx has about d
exact digits: there exists τ = ±(0.0 . . . 0nd+1nd+2 . . . ) such
that ˆx = x + τ
true value is usually not known → use estimates or bounds on
the error
relative error can be taken relative to the approximate value
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Example/exercise - Implement!
Stirling’s approximation for factorials:
Sn =
√
2πn
n
e
n
≈ n!, n = 1, 2, . . .
where e = exp(1).
Relative error (Sn − n!)/n!:
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Errors: data and computational
compute f(x) for f : R → R
x ∈ R is the true value
f(x) true/desired result
ˆx approximate input
ˆf approximate result
total error:
ˆf(ˆx) − f(x) = (ˆf(ˆx) − f(ˆx)) + (f(ˆx) − f(x))
= computational error + propagated data error
the algorithm has no effect on propagated error
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Computational error
is sum of:
truncation error = (true result) - (result of the algorithm using
exact arithmetic)
Example: considering only the first terms of an infinite Taylor
series; stopping before convergence
rounding error = (result of the algorithm using exact
arithmetic) - (result of the algorithm using limited precision
arithmetic)
Example: π ≈ 3.14 or π ≈ 3.141593
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Finite difference approximation
f (x) = lim
h→0
f(x + h) − f(x)
h
≈
f(x + h) − f(x)
h
, for some small h > 0
truncation error: f (x) −
f(x+h)−f(x)
h ≤ Mh/2 where |f (t)| ≤ M
for t in a small neighborhood of x (HOMEWORK, 5p)
rounding error: 2 /h, for being the precision
total error is minimized for h ≈ 2
√
/M
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Figure : Total computational error as a tradeoff between truncation and
rounding error (from Heath - Scientific computing)
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Error analysis
For y = f(x), for f : R → R an approximate ˆy result is obtained.
forward error: ∆y = ˆy − y
backward error: ∆x = ˆx − x, for f(ˆx) = ˆy
x
f //
OO
backward error

ˆf

y
OO
forward error

= f(x)
ˆx
f // ˆy = ˆf(x) = f(ˆx)
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Compute f(x) = ex
for x = 1. Use the first 4 terms from Taylor
expansion:
ˆf(x) = 1 + x +
x2
2
+
x3
6
take "true" value: f(x) = 2.716262 and compute
ˆf(x) = 2.666667, then
forward error: |∆y| = 0.051615, or a relative f. error of about
2%
backward error: ˆx = lnˆf(x) = 0.989829 ⇒ |∆x| = 0.019171,
or a relative b. error of 2%
these are two perspectives on assessing the accuracy
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Exercise
Consider the general Taylor series with limit e:
∞
n=0
1
n!
= e
How many terms are needed for an approximation of e to three
decimal places?
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Backward error analysis
idea: approximate result is the exact solution of a modified
problem
how far from the original problem is the modified version?
how much error in the input data would explain all the error in
the result?
an approximate solution is good if it is an exact solution for a
nearby problem
backward analysis is usually easier
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Sensitivity and conditioning
insensitive (well-conditioned) problem: relative changes in
input data causes similar relative change in the result
large changes in solution for small changes in input data
indicate a sensitive (ill-conditioned) problem;
condition number:
cond =
absolute relative change in solution
absolute relative change in input
=
|∆y/y|
|∆x/x|
if cond >> 1 the problem is sensitive
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
condition number is a scale factor for the error:
relative forward err = cond × relative backward err
usually, only upper bounds of the cond. number can be
estimated, cond ≤ C, hence
relative forward err ≤ C × relative backward err
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
ˆx = x + ∆x
forward error: f(x + ∆x) − f(x) ≈ f (x)∆x, for small enough
∆x
relative forward error: ≈
f (x)∆x
f(x)
⇒ cond ≈
xf (x)
f(x)
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
ˆx = x + ∆x
forward error: f(x + ∆x) − f(x) ≈ f (x)∆x, for small enough
∆x
relative forward error: ≈
f (x)∆x
f(x)
⇒ cond ≈
xf (x)
f(x)
Example: tangent function is sensitive in neighborhood of π/2
tan(1.57079) ≈ 1.58058 × 105
; tan(1.57078) ≈ 6.12490 × 104
for x = 1.57079, cond ≈ 2.48275 × 105
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Stability
an algorithm is stable if is relatively insensitive to
perturbations during computation
stability of algorithms is analogous to conditioning of problems
backward analysis: an algorithm is stable if the result
produced is the exact solution of a nearby problem
stable algorithm: the effect of computational error is no worse
than the effect of small error in input data
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Accuracy
accuracy closeness of the result to the true solution of the
problem
depends on the conditioning of the problem AND on the
stability of the algorithm
stable algorithm + well-conditioned problem = accurate results
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
CPUs
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Number representation
internally, all data are represented in binary format (each digit
can be either 0 or 1, e.g. 1011001...)
bit, nybble, byte
word → specific to architecture: 1, 2, 4, or 8 bytes
integers:
unsigned (≥ 0): on n bits: 0, . . . , 2n
− 1. The stored
representation (for 1 byte) is b7b6b5b4b3b2b1b0 for a value
x = 7
i=0 bi2i
.
signed: 1 bit for sign, rest for the absolute value;
−2n−1
, . . . , 0, . . . , 2n−1
− 1. The stored representation (for 1
byte) is b7b6b5b4b3b2b1b0 for a value x = b7(−27
) + 6
i=0 bi2i
.
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Floating-point numbers
like in scientific notation: mantissa × radixexponent
, e.g.
2.35 × 103
formally
x = ± b0 +
b1
β
+
b2
β2
+ · · · +
bp−1
βp−1
× βE
where
β is the radix (or base)
p is the precision
L ≤ E ≤ U are the limits of the exponent
0 ≤ bk ≤ β
mantissa: m = b0b1 . . . bp−1; fraction: b1b2 . . . bp−1
the sign, mantissa and exponent are stored in fixed-sized
fields (the radix is implicit for a given system, β = 2 usually)
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Normalization:
b0 0 for all x 0
mantissa m satisfies 1 ≤ m < β
ensures unique representation, optimal use of available bits
Internal representation (on 64 bits - "double precision", binary
representation):
x = sign | exponent | fraction = b63 b62 . . . b52 b51 . . . b0
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Properties:
only a finite number of discrete values can be represented
total number of floating point numbers representable in
normalized format is
2(β − 1)βp−1
(U − L + 1) + 1
(Q: can you justify the result?)
undeflow level (smallest number): UFL = βL
overflow level (largest number): OFL = βU+1
(1 − β−p
)
not all real numbers can be represented exactly:
machine numbers
rounding → rounding error
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Example: let β = 2, p = 3, L = −1, U = 1, there are 25 distinct
numbers that can be represented:
UFL = 0.510; OFL = 3.510
note the non-uniform coverage
∀x ∈ R, fl(x) is the floating point representation; x − fl(x) is
the rounding error
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Rounding rules
chop = round toward zero: truncate the base−β
representation after p − 1st digit
round to nearest: fl(x) is the closest machine number to x
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Machine precision
machine precision, mach
with chopping: mach = β1−p
with rounding to nearest: mach = 1
2
β1−p
called also unit roundoff: the smallest number such that
fl(1 + ) > 1
maximum relative error of representation
fl(x) − x
x
≤ mach
usually 0 < UFL < mach < OFL
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Machine precision - example
For β = 2, p = 3, L = −1, U = 1,
mach = (0.01)2 = (0.25)10 with chopping
mach = (0.001)2 = (0.125)10 with rounding to nearest
The usual case (IEEE fp systems):
mach = 2−24
≈ 10−7
in single precision
mach = 2−53
≈ 10−16
in double precision
→ about 7 and 16 decimals of precision, respectively
(in R: p-value < 2.2e − 16)
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Gradual underflow
to improve representation of numbers around 0 - use
subnormal (or denormalized) numbers
when exponent is at minimum, alow leading digits to be 0
subnormals are less precise
→ gradual underflow
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Special values
IEEE standard:
Inf: infinity; the result of 1/0
NaN: the result of 0/0 or Inf/Inf
special representation of the exponent field
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Floating-point arithmetic
addition/subtraction: denormalization might be required:
3.52 × 103
+ 1.97 × 105
= 0.0352 × 105
+ 1.97 × 105
=
2.0052 × 105
→ might cause loss of some digits
multiplication/division: the result may not be representable
overflow is more serious than underflow: how to approximate
large numbers?
for underflow, the result may be approximated by 0
in FP arithm. addition and multiplication are commutative but
not associative: if is slightly smaller than mach, then
(1 + ) + = 1, but 1 + ( + ) > 1
ideally, x flop y = fl(xopy); IEEE standard ensures this for
within range results
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Example: divergent series
∞
n=1
1
n
in FP arithm, the sum of the series is finite;
depending on the system, this is because:
after a while, the sum overflows
1/n underflows
for all n such that
1
n
< mach
n−1
k=1
1
k
the sum does not change anymore
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Cancellation
subtracting 2 numbers of the same magnitude usually cancels
the most significant digits:
1.92403 × 102
− 1.92275 × 102
= 1.28000 × 10−1
→ only 3
significant digits
let > 0 be slightly smaller than mach, then (1 + ) − (1 − )
yields 0 in FP arithmetic, instead of 2 .
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Cancellation - example
For the quadratic equation, ax2
+ bx + c = 0, the two solutions are
given by
x1,2 =
−b ±
√
b2 − 4ac
2a
Problems:
for very large/small coefficients, the terms b2
or 4ac may
over-/underflow → rescale coeficients by max{a, b, c}.
cancellation between −b and
√
· can be avoided by computing
one root using x = 2c
−b
√
b2−4ac
Exercise: let x1 = 2000, x2 = 0.05 be the roots of a quadratic
equation. Compute the coefficients and then use the above
formulas to retrieve the roots. Try roots() function in Matlab.
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
Cancellation - another example
P(X) = (X − 1)6
= X6
− 6X5
+ 15X4
− 20X3
+ 15X2
− 6X + 1.
What happens around X = 1?
1 epsilon = [.01, .005, .001];
2 for k=1:3
3 x = linspace(1−epsilon(k), 1+epsilon(k), 100);
4 px = x.^6 − 6*x.^5 + 15*x.^4 − 20*x.^3 + 15*x.^2 ...
− 6*x + 1;
5 px0 = (x − 1).^6;
6 subplot(2, 3, k);
7 plot(x, px, '−b', x, zeros(1,100), '−r');
8 axis([1−epsilon(k), 1+epsilon(k), −max(abs(px)), ...
max(abs(px))]);
9 subplot(2, 3, k+3);
10 plot(x, px0, '−b', x, zeros(1, 100), '−r' );
11 axis([1−epsilon(k), 1+epsilon(k), −max(abs(px0)), ...
max(abs(px0))]);
12 end
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
...mathematically equivalent, but numerically different...
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
HOMEWORK
Let x = [a, b] ∈ R2
and let p be its Euclidean norm, p =
√
a2 + b2.
However, using this formula is prone to under- and over-flow errors.
show that min{|a|, |b|} ≤ p ≤
√
2 max{|a|, |b|}
implement in Matlab a procedure that would avoid
unnecessary under-/over-flows Hint: p = c (a/c)2 + (b/c)2.
Find a suitable c...
Vlad Bi7740: Scientific computing
Introduction
Sensitivity and conditioning
Computer arithmetic
In Matlab...
you can change the format of FP in output using format
option
mach is returned by the function eps():
single precision: eps(single(1)) gives
1.1921e − 07 = 2−23
double precision: eps(double(1)) gives
2.2204e − 16 = 2−52
to obtain the smallest or largest single/double precision
numbers, use realmin('single'), realmin('double'),
realmax('single'), realmax('double')
you have the special constants Inf and NaN
Vlad Bi7740: Scientific computing