9) Function errors - Computational Methods for Scientists Fall 26

Today¶

Conditioning and well posedness
Relative and absolute errors
Condition number

using Plots
default(linewidth=4, legendfontsize=12, xtickfontsize=12, ytickfontsize=12)

1. Conditioning and well-posedness¶

We say that a mathematical function $f(x)$ is well conditioned if small changes in $x$ produce small changes in $f(x)$ . (What we mean by “small” will be made more precise.)

The function $f(x)$ may represent a simple expression such as

$f(x) := 2 x$
$f(x) := \sqrt{x}$
$f(x) := \log{x}$
$f(x) := x - 1$

A function may also represent something more complicated, implementable on a computer or by physical experiment.

Find the positive root of the polynomial $t^2 + (1-x)t - x.$
Find the eigenvectors of the matrix
$A(x) = \begin{bmatrix} 1 & 1 \\ 0 & x \end{bmatrix} .$
(1)
Find how much the bridge flexes when the truck of mass $x$ drives over it.
Find the length of the rubber band when it finally snaps, as a function of temperature $x$ during manufacturing.
Find the time at which the slope avalanches as a function of the wind speed $x$ during the storm.
Find the probability that the slope avalanches in the next 48 hours as a function of the wind speed $x$ during the storm.
Find the probability that the hurricane makes landfall as a function of the observations $x$ .

Specification¶

Some of these problems are fully-specified
Others involve sophisticated models and ongoing community research problems.
In some cases, the models that are computable may incur greater uncertainty than the underlying system. In such cases, an analog experiment might produce smoother variation of the output as the problem data $x$ are varied.
In others, the model might be better behaved than what it seeks to model.
Some of these problems may be ill-posed

Well-posedness¶

A problem is said to be well-posed if

a solution exists,
the solution is unique, and
the solution depends continuously on the problem specification.

Mathematically, continuous variation in part 3 can be arbitrarily fast, but there may be measurement error, manufacturing tolerances, or incomplete specification in real-world problems, such that we need to quantify part 3. This is the role of conditioning.

Example:¶

Computing $e^x$

e^x = \sum_{k=0}^{\infty} x^k/k! = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dotsb

(2)

function myexp(x)
    sum = 1
    for k in 1:100
        sum += x^k / factorial(big(k))
    end
    return sum
end
myexp(1) - exp(1)

1.445646891729250136554224997792468345749669676277240766303534163549513721620773e-16

function myexp(x)
    sum = 0
    term = 1
    n = 1
    while sum + term != sum
        sum += term
        term *= x / n
        n += 1
    end
    sum
end
myexp(1) - exp(1)

4.440892098500626e-16

How accurate is it?

plot(myexp, xlims=(-2, 2))

Loading...

plot(exp, xlims=(-1e-15, 1e-15), linestyle=:solid, label = "exp")
plot!(myexp, xlims=(-1e-15, 1e-15), linestyle=:dash, label="myexp") # the bang operator overlays the second plot in the same figure

Loading...

GKS: Possible loss of precision in routine SET_WINDOW

What’s happening?

We’re computing $f(x) = e^x$ for values of $x$ near zero.
This function is well approximated by $1 + x$ .
Values of $y$ near 1 cannot represent every value.
After rounding, the error in our computed output $\tilde f(x)$ is of order $\epsilon_{\text{machine}}$ .

2. Absolute and Relative Errors¶

Absolute Error¶

\lvert \tilde f(x) - f(x) \rvert

(3)

Relative Error¶

\frac{\lvert \tilde f(x) - f(x) \rvert}{\lvert f(x) \rvert}

(4)

Suppose I want to compute $e^x - 1$

plot([x -> myexp(x) - 1 , x -> x],
     xlims=(-1e-15, 1e-15), label=["myexp" "y=x"])

Loading...

What now?

We’re capable of representing outputs with 16 digits of accuracy
Yet our algorithm myexp(x) - 1 can’t find them
We can’t recover without modifying our code

Modify the code

function myexpm1(x)
    sum = 0
    term = x
    n = 2
    while sum + term != sum
        sum += term
        term *= x / n
        n += 1
    end
    sum
end

myexpm1 (generic function with 1 method)

plot(myexpm1, xlims=(-1e-15, 1e-15), label="myexpm1")

Loading...

Plot relative error¶

function relerror(x, f, f_ref)
    fx = f(x)
    fx_ref = f_ref(x)
    max(abs(fx - fx_ref) / abs(fx_ref), 1e-17)
end
badexpm1(t) = exp(t) - 1
plot(x -> relerror(x, badexpm1, expm1), yscale=:log10, xrange=(-1e-15, 1e-15), label = "rel error")

Loading...

Floating point representation is relative (see float.exposed)¶

IEEE 754 double-precision floating-point format

Let $\operatorname{fl}$ round to the nearest floating point number.

\operatorname{fl}(x) = x (1 + \epsilon), \quad \text{where} |\epsilon| \le \epsilon_{\text{machine}}

(5)

This also means that the relative error in representing $x$ is small:

\frac{|\operatorname{fl}(x) - x|}{|x|} \le \epsilon_{\text{machine}}

(6)

plot(x -> (1 + x) - 1, xlims=(-1e-15, 1e-15), label = "y = 1 + x - 1")
plot!(x -> x, label="y = x")

Loading...

Exact arithmetic, correctly rounded¶

Take an elementary math operation $*$ (addition, subtraction, multiplication, division), and the discrete operation that our computers perform, $\circledast$ . Then

x \circledast y := \operatorname{fl}(x * y)

(7)

with a relative accuracy $\epsilon_{\text{machine}}$ ,

\frac{|(x \circledast y) - (x * y)|}{|x * y|} \le \epsilon_{\text{machine}} .

(8)

How do operations compose?¶

Is this true?

\frac{\Big\lvert \big((x \circledast y) \circledast z\big) - \big((x * y) * z\big) \Big\rvert}{|(x * y) * z|} \le^? \epsilon_{\text{machine}}

(9)

f(x; y=1, z=-1) = (x+y)+z # The best arbitrary numbers are 0, 1, and -1
plot(x -> abs(f(x) - x)/abs(x), xlims=(-1e-15, 1e-15), label = "abs(f(x) - x)/abs(x)")

Loading...

3. Condition number¶

What sort of functions cause small errors to become big?

Consider a function $f: X \to Y$ and define the absolute condition number

\hat\kappa = \lim_{\delta \to 0} \max_{|\delta x| < \delta} \frac{|f(x + \delta x) - f(x)|}{|\delta x|} = \max_{\delta x} \frac{|\delta f|}{|\delta x|}.

(10)

If $f$ is differentiable, then $\hat\kappa = |f'(x)|$ .

Floating point offers relative accuracy, so it’s more useful to discuss relative condition number,

\kappa = \max_{\delta x} \frac{|\delta f|/|f|}{|\delta x|/|x|} = \max_{\delta x} \Big[ \frac{|\delta f|/|\delta x|}{|f| / |x|} \Big]

(11)

or, if $f$ is differentiable,

\kappa = |f'(x)| \frac{|x|}{|f|} .

(12)