2025-08-28 Conditioning and errors#

Last Time#

  • Intro to Julia

  • Intro to floating point

  • Verification and validation

Today#

  • Discussion of JFE Editorial Policy

  • Conditioning and well posedness

  • Relative and absolute errors

  • Condition numbers

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

Why is reliability so important?#

Modeling and simulation is used to plan at the greatest scale.

  • Dams, bridges

  • Toxic/nuclear waste management

  • Nuclear stockpile stewardship

  • Climate policy/water resource management/reinsurance

  • Airplanes, etc.

We can’t just 3d print a bunch of widgets and load test them. Tests are expensive ($ billions over many years) or worse.

“Will you sign the blueprint?”#

– Ivo Babuška

Conditioning#

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\)

Measurement and conditioning#

Geodesy#

  • Dual-band differential GPS can be used to measure locations on earth with 1cm accuracy.

  • It’s used to track motion of tectonic plates.

Smaller length scales#

  • Would you use GPS to survey for a construction project? Why or why not?

Conditioning#

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

    \[\begin{split} A(x) = \begin{bmatrix} 1 & 1 \\ 0 & x \end{bmatrix} .\end{split}\]

  • Find much the bridge flexes when the truck of mass \(x\) drives over it.

  • Find how long it takes to clean up when the toddler knocks the glass off the counter, as a function of the chair location \(x\).

  • 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

  1. a solution exists,

  2. the solution is unique, and

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

Computing \(e^x\)#

\[ e^x = \sum_{k=0}^{\infty} x^k/k! \]
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))
plot([exp myexp], xlims=(-1e-15, 1e-15))
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}}\).

Absolute Error#

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

Relative Error#

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

Suppose I want to compute \(e^x - 1\)#

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

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

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

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

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

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

This also means that the relative error in representing \(x\) is small:

\[ \frac{|\operatorname{fl}(x) - x|}{|x|} \le \epsilon_{\text{machine}} \]
plot(x -> (1 + x) - 1, xlims=(-1e-15, 1e-15))
plot!(x -> x)

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)\]

with a relative accuracy \(\epsilon_{\text{machine}}\),

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

Seems easy, 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}} \]
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))

Which operation caused the error?#

  1. \(\texttt{tmp} = \operatorname{fl}(x + 1)\)

  2. \(\operatorname{fl}(\texttt{tmp} - 1)\)

Use Julia’s BigFloat

@show typeof(big(.1))
@show big(.1)          # Or BigFloat(.1); parsed as Float64, then promoted
@show BigFloat(".1");  # Parse directly to BigFloat
typeof(big(0.1)) = BigFloat
big(0.1) = 0.1000000000000000055511151231257827021181583404541015625
BigFloat(".1") = 0.1000000000000000000000000000000000000000000000000000000000000000000000000000002
tmp = 1e-15 + 1
tmp_big = big(1e-15) + 1 # Parse as Float64, then promote
abs(tmp - tmp_big) / abs(tmp_big)
1.102230246251563524952071662733800140614440125894379682676737388538642032894338e-16
r = tmp - 1
r_big = big(tmp) - 1
abs(r - r_big) / abs(r_big)
0.0

This last experiment shows that step 2 is exact (there is no rounding). The only approximation occurs when rounding in step 1, and yet I blame step 2 for the error.

Conditioning#

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|}. \]
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] \]
or, if \(f\) is differentiable,
\[ \kappa = |f'(x)| \frac{|x|}{|f|} . \]

Condition numbers#

\[ \kappa = |f'(x)| \frac{|x|}{|f|} \]
f(x) = x - 1; fp(x) = 1
plot(x -> abs(fp(x)) * abs(x) / abs(f(x)), xlims=(0, 2))

Back to \(f(x) = e^x - 1\)#

f(x) = exp(x) - 1
fp(x) = exp(x)
plot(x -> abs(fp(x)) * abs(x) / abs(f(x)))

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

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))
cond1(x) = abs(1/(1+x) * x / log1p(x))
cond2(x) = abs(1/x * x / log(x))
plot([cond1 cond2], xlims=(-1, 2), ylims=(0, 100))

Reliable = well-conditioned and stable#

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: fracture, chaotic systems, combustion

Algorithms f(x) can be unstable#

  • Unreliable, though sometimes practical

  • Unstable algorithms are constructed from ill-conditioned parts

An ill-conditioned problem from Paul Olum#

From Surely You’re Joking, Mr. Feynman (page 113)

So Paul is walking past the lunch place and these guys are all excited. “Hey, Paul!” they call out. “Feynman’s terrific! We give him a problem that can be stated in ten seconds, and in a minute he gets the answer to 10 percent. Why don’t you give him one?” Without hardly stopping, he says, “The tangent of 10 to the 100th.” I was sunk: you have to divide by pi to 100 decimal places! It was hopeless.

What’s the condition number?

\[ \kappa = |f'(x)| \frac{|x|}{|f|} \]
  • \(f(x) = \tan x\)

  • \(f'(x) = 1 + \tan^2 x\)

\[ \kappa = \lvert x \rvert \Bigl( \lvert \tan x \rvert + \bigl\lvert \frac{1}{\tan x} \bigr\rvert \Bigr)\]
tan(1e100)
-0.4116229628832498
tan(BigFloat("1e100", precision=100)) # check the precision!
2.738770191194700792231947565702313204942361296608129856300353713193607246577622