10) Conditioning

Today:

Recap from last time
Condition numbers

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

1. Recap from last time¶

Last time we saw:

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)

(1)

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

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

(2)

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}}

(3)

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 = "rel error")

Loading...

Activity¶

Use Jupyter and Julia
Look at how fast series converge when taking only finitely many terms
Explore instability, as is occuring for large negative x above, but not for the standard library expm1

function myexp(x, k)
    sum = 0
    term = 1
    n = 1
    # modify so at most k terms are summed
    while sum + term != sum
        sum += term
        term *= x / n
        # stop at first k terms
        if n == k
            break
        end
        n += 1
    end
    sum
end

myexp(-30, 1000), exp(-30) # stopping at the 1000th term Vs the Julia built-in

(6.1030424788918156e-6, 9.357622968840175e-14)

rel_error(x, k) = abs(myexp(x, k) - exp(x)) / exp(x)
ks = 2 .^ (0:10)  # [1, 2, 4, ..., 1024];

11-element Vector{Int64}:
    1
    2
    4
    8
   16
   32
   64
  128
  256
  512
 1024

plot(ks, k -> rel_error(-10, k), xscale=:log10, yscale=:log10, label = "rel_error", markershape = :circle, markercolor = "brown")

Loading...

plot(x -> rel_error(x, 1000) + 1e-17, xlims=(-20, 20), yscale=:log10)

Loading...

What happened? Let’s look at the terms for positive and negative $x$

function expterms(x, k=50)
    term = 1.
    terms = [term]
    for n in 1:k
        term *= x / n
        push!(terms, term)
    end
    terms
end
x = -10
@show sum(expterms(x)) - exp(x)
@show (sum(expterms(x)) - exp(x)) / exp(x)
expterms(x)

sum(expterms(x)) - exp(x) = 1.2092180894291739e-12
(sum(expterms(x)) - exp(x)) / exp(x) = 2.6634800885273225e-8

51-element Vector{Float64}:
     1.0
   -10.0
    50.0
  -166.66666666666669
   416.66666666666674
  -833.3333333333335
  1388.8888888888891
 -1984.1269841269846
  2480.1587301587306
 -2755.7319223985896
  2755.7319223985896
 -2505.2108385441725
  2087.6756987868107
     ⋮
    -4.902469756513546e-8
     1.2256174391283866e-8
    -2.9893108271424062e-9
     7.117406731291443e-10
    -1.655210867742196e-10
     3.761842881232264e-11
    -8.359650847182808e-12
     1.81731540156148e-12
    -3.866628513960596e-13
     8.055476070751242e-14
    -1.64397470831658e-14
     3.28794941663316e-15

x = -10
@show exp(x)
bar(expterms(x), label = "expterms")

exp(x) = 4.5399929762484854e-5

Loading...

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|}.

(4)

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]

(5)

or, if $f$ is differentiable,

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

(6)

f(x) = x - 1; fprime(x) = 1
plot(x -> abs(fprime(x)) * abs(x) / abs(f(x)), xlims=(0, 2), label = "kappa for x-1")

Loading...

Example:¶

Back to $f(x) = e^x - 1$

f(x) = exp(x) - 1
fprime(x) = exp(x)
plot(x -> abs(fprime(x)) * abs(x) / abs(f(x)), label = "kappa for exp(x)")

Loading...

What does it mean?

The function $f(x) = e^x - 1$ is well-conditioned
The function $f_1(x) = e^x$ is well-conditioned
The function $f_2(x) = x - 1$ is ill-conditioned for $x \approx 1$

The algorithm is unstable¶

f(x) = exp(x) - 1 is unstable
Algorithms are made from elementary operations
Unstable algorithms do something ill-conditioned

A stable algorithm¶

We used the series expansion previously.

accurate for small $x$
less accurate for negative $x$ (see activity)
we could use symmetry to fix
inefficient because we have to sum lots of terms

Standard math libraries define a more efficient stable variant, $\texttt{expm1}(x) = e^x - 1$ .

expm1(-40) + 1 # this way, using the built-in, stable variant, it is exactly 0

0.0

exp(-40) # almost zero

4.248354255291589e-18

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

Loading...

Another example:¶

$\log(1 + x)$

What is the condition number of $f(x) = \log(1 + x)$ for $x \approx 0$ ?

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

Loading...

condlog1p(x) = abs(1/(1+x) * x / log1p(x))
condlog(x) = abs(1/x * x / log(x))
plot([condlog1p condlog], xlims=(-1, 2), ylims=(0, 100), label = ["cond of log1p" "cond of log"])

Loading...

Reliable = well-conditioned and stable algorithm¶

Mathematical functions $f(x)$ can be ill-conditioned (big $\kappa$ )¶

Modeling is how we turn an abstract question into a mathematical function
We want well-conditioned models (small $\kappa$ )
Some systems are intrinsically sensitive (chaotic): fracture, chaotic systems, combustion

Algorithms `f(x)` can be unstable¶

Unreliable, though sometimes practical
Unstable algorithms are constructed from ill-conditioned parts

Exercise 10.1:¶

Find a function $f(x)$ that models something you’re interested in
Plot its condition number $\kappa$ as a function of $x$
Plot the relative error (using single or double precision; compare using Julia’s big)
Is the relative error ever much bigger than $\kappa \epsilon_{\text{machine}}$ ?
If so, can you find what caused the instability in your algorithm?

Exercise 10.2:¶

How big are these condition numbers?

$f(x) = x+c$
$f(x) = cx$
$f(x) = x^p$
$f(x) = e^x$
$f(x) = \log x$
$f(x) = \sin x$

1. Recap from last time¶

Exact arithmetic, correctly rounded¶

How do operations compose?¶

Activity¶

Condition number¶

Example:¶

The algorithm is unstable¶

A stable algorithm¶

Another example:¶

Reliable = well-conditioned and stable algorithm¶

Mathematical functions f(x)f(x)f(x) can be ill-conditioned (big κ\kappaκ)¶

Algorithms f(x) can be unstable¶

Exercise 10.1:¶

Further reading: FNC Introduction¶

Exercise 10.2:¶

Mathematical functions $f(x)$ can be ill-conditioned (big $\kappa$ )¶

Algorithms `f(x)` can be unstable¶