¶

Sources of Error¶

Error can come from many sources when formulating problems and/or applying numerical methods:

• Model/Data Error

• Discretization Error

• Floating Point Error

• Convergence Error

Objectives¶

• Understand the different sources of error

• Explore some simple approaches to error analysis

• Quantify errors

  • absolute error

  • relative error

  • precision

• Long term goals

  • Use error estimates to control accuracy/reliability of solutions

  • Understand errors so you can believe and justify your solutions

Model and Data Error¶

Errors in fundamental formulation

• SIR model

  • simplistic averaged model of interactions

    $$\dot{I} = \alpha SI - \beta I$$

    (1)

  • basic model predicts a single peak

• Data Error - Inaccuracy in measurement or uncertainties in parameters

Unfortunately we cannot control model and data error directly but we can use methods that may be more robust in the presense of these types of errors.

Discretization or Truncation Error¶

Errors arising from approximating a function with a simpler function, e.g. Using the approximation $\sin(x) \approx x$ when $|x| \approx 0$.

Floating Point Error¶

Errors arising from approximating real numbers with finite-precision numbers and arithmetic.

Convergence Error¶

In some instances an algorithm is developed that will take a current approximation and then find an improvement on the current approximation. In some instances the errors generated in each indivudal step can accumulate or become magnified after repeating the algorithm a number of times.

Basic Definitions¶

Before exploring the different kinds of error, it is important to first define the ways that error is measured. Given a true value of a function $f$ and an approximate solution $F$ define:

Absolute Error:

Relative Error:

Note: these definitions assume $f$ and $F$ are scalar valued. However these definitions are readily extended to more complicate objects such as vectors or matrices with appropriate norms.

Decimal Precision¶

This definition of relative error provides a convenient estimate for the number of digits of decimal precision $p$

given a relative error $r$, the precision $p$ is the largest integer such that

Example

• if $r = 0.001 < 5\times10^{-3}$ has $p=3$ significant digits

• if $r = 0.006 < 5\times10^{-2}$ has $p=2$ significant digits (because this error would cause rounding up)

Example¶

let

Big-O Notation¶

In many situations an approximation will have a parameter associated with it, and the value of the parameter is often chosen to insure that the error is reasonable in a given situation. In such circumstances we often want to know the impact on the error if we change the value of the parameter. This leads to the definition of Big-O notation:

if and only if

In practice we use Big-O notation to say something about how the terms we may have left out of a series might behave. We saw an example earlier of this with the Taylor’s series approximations.

Example:¶

Consider approximating a differentiable function $f(x)$ by its Taylor polynomial (truncated Taylor series) expanded around $x_0=0$., i.e.

where

assuming the Taylor series converges

or for the case $f(x)=\sin(x)$ expanded around $x_0=0$

For $N=2$, we can then write $F(x)$ as

so our true function is

or the absolute error

We can also develop rules for error propagation based on Big-O notation:

In general, there are two theorems that do not need proof and hold when the value of x is large:

Let

then

and

On the other hand, if we are interested in small values of x, say $\Delta x$, the above expressions can be modified as follows:

then

and

Note: In this case we suppose that at least the polynomial with $k = \max(n, m)$ has the following form:

or

so that there is an $\mathcal{O}(1)$ term that guarantees the existence of $\mathcal{O}(\Delta x^r)$ in the final product.

To get a sense of why we care most about the power on $\Delta x$ when considering convergence the following figure shows how different powers on the convergence rate can effect how quickly we converge to our solution. Note that here we are plotting the same data two different ways. Plotting the error as a function of $\Delta x$ is a common way to show that a numerical method is doing what we expect and exhibits the correct convergence behavior. Since errors can get small quickly it is very common to plot these sorts of plots on a log-log scale to easily visualize the results. Note that if a method was truly of the order $n$ that they will be a linear function in log-log space with slope $n$.

Behavior of error as a function of $\Delta x$¶

Discretization Error¶

Taylor’s Theorem: Let $f(x) \in C^{N+1}[a,b]$ and $x_0 \in [a,b]$, then for all $x \in (a,b)$ there exists a number $c = c(x)$ that lies between $x_0$ and $x$ such that

where $T_N(x)$ is the Taylor polynomial approximation

and $R_N(x)$ is the residual (the part of the series we left off)

Note¶

The residual:

depends on the $N+1$ order derivative of $f$ evaluated at an unknown value $c\in[x,x_0]$.

If we knew the value of $c$ we would know the exact value of $R_N(x)$ and therefore the function $f(x)$. In general we don’t know this value but we can use $R_N(x)$ to put upper bounds on the error and to understand how the error changes as we move away from $x_0$.

Start by replacing $x - x_0$ with $\Delta x$. The primary idea here is that the residual $R_N(x)$ becomes smaller as $\Delta x \rightarrow 0$ (at which point $T_N(x) = f(x_0)$).

and $R_N(x)$ is the residual (the part of the series we left off)

where $M$ is an upper bound on

Example 1¶

$f(x) = e^x$ with $x_0 = 0$ on the interval $x\in(-1,1)$

Using this we can find expressions for the relative and absolute error as a function of $x$ assuming $N=2$.

Derivatives:

Taylor polynomials:

Remainders:

Accuracy:

We can also use the package sympy which has the ability to calculate Taylor polynomials built-in!

Lets plot this numerically for a section of $x$.

Example 2¶

Approximate

using $x_0 = 1$ to the 3rd Taylor series term on the inverval $x\in[1,\infty)$

plot this problem¶

Computational Issue #1: Accuracy... how many terms?¶

Given a Taylor Polynomial approximation of an arbitrary function $f(x)$, how do we determine how many terms are required such that $|R_N(x)|<tol$. And how do we determine the tolerance?

Computational Issue #2 Efficiency... Operation counts for polynomial evaluation¶

Given

or

what is the most efficient way to evaluate $P_N(x)$? (i.e. minimize number of floating point operations)

Consider two ways to write $P_3$

• The standard way:

• using nested multiplication (aka Horner’s Method):

Consider how many operations it takes for each...

Note: here we’re just counting multiplications as they will dominate the flop count

Adding up all the operations we can in general think of this as a pyramid (it’s really the triangle numbers)

Thus we can estimate that the algorithm written this way will take approximately $O(N^2 / 2)$ operations to complete.

Looking at nested iteration, however:

Here we find that the method is $O(N)$ compared to the first evaluation which $O(N^2)$ (we usually drop the 2 in these cases). That’s a huge difference for large $N$!

Algorithm¶

Fill in the function and implement Horner’s method:

def eval_poly(p, x):
    '''Evaluates a polynomial using Horner's method given coefficients p at x
    
      The polynomial is defined as
    
        P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]
        
    Parameters:
        p: list or numpy array of coefficients 
        x:  scalar float
        
    returns:
        P(x):  value of the polynomial at point x (float)
    '''
    pass

Convergence Error¶

In some circumstances a formula or algorithm is applied repeatedly as a way to obtain a final approximation. Usually, the errors that occur at each individual step are small. By repeating the algorithm, though, the errors can sometimes grow or become magnified.

As example of this phenomena is given below. The values of the terms in a difference equation are calculated, $$

$$

The true solution to the difference equation is $y_n = \left(\frac{1}{5}\right)^n$, where $n=$0, 1, 2, $\ldots$

And plot the result

Simply looking at the exact solution, the sequence of numbers generated by the difference equation above should get very close to zero. Instead, the numbers in the sequence initially get closer to zero, but at some point they begin to grow and get larger. An underlying problem is that the computer is not able to store the numbers exactly. The second number in the sequence, $y_1=\frac{1}{5}$ has a small error, and the computer stores it as $y_1 = \frac{1}{5}+\epsilon$ where $\epsilon$ is some small error.

Each time a new number in the loop is generated, the error is multiplied. For example, after the first iteration $y_2$ is

After the second time through the loop, the value of $y_3$ is

Even though the value of $\epsilon$ is very close to zero, every iteration makes the error grow. Repeated multiplication will result in a very large number.

The error associated with the initial representation of the number $\frac{1}{5}$ is a problem with the way a digital computer stores floating point numbers. In most instances the computer cannot represent a number exactly, and the small error in approximating a given number can give rise to other problems.

Floating Point Error¶

Errors arising from approximating real numbers with finite-precision numbers

or $\frac{1}{3} \approx 0.333333333$ in decimal, results from a finite number of registers to represent each number.

Floating Point Systems¶

Numbers in floating point systems are represented as a series of bits that represent different pieces of a number. In normalized floating point systems there are some standard conventions for what these bits are used for. In general the numbers are stored by breaking them down into the form

where

1. ± is a single bit representing the sign of the number

2. $d_1 . d_2 d_3 d_4 \ldots d_p$ is called the mantissa. Note that technically the decimal could be moved but generally, using scientific notation, the decimal can always be placed at this location. The digits $d_2 d_3 d_4 \ldots d_p$ are called the fraction with $p$ digits of precision. Normalized systems specifically put the decimal point in the front like we have and assume $d_1 \neq 0$ unless the number is exactly 0.

3. $\beta$ is the base. For binary $\beta = 2$, for decimal $\beta = 10$, etc.

4. $E$ is the exponent, an integer in the range $[E_{\min}, E_{\max}]$

The important points on any floating point system is that

1. There exist a discrete and finite set of representable numbers

2. These representable numbers are not evenly distributed on the real line

3. Arithmetic in floating point systems yield different results from infinite precision arithmetic (i.e. “real” math)

Properties of Floating Point Systems¶

All floating-point systems are characterized by several important numbers

• Smalled normalized number (underflow if below - related to subnormal numbers around zero)

• Largest normalized number (overflow if above)

• Zero

• Machine $\epsilon$ or $\epsilon_{\text{machine}}$

• inf and nan, infinity and Not a Number respectively

Example: Toy System¶

Consider the toy 2-digit precision decimal system (normalized)

with $E \in [-2, 0]$.

Number and distribution of numbers

1. How many numbers can we represent with this system?

2. What is the distribution on the real line?

3. What is the underflow and overflow limits?

4. What is the smallest number $\epsilon_{mach}$ such that $1+\epsilon_{mach} > 1$?

How many numbers can we represent with this system?

• sign bit: 2

• $d_1$: 9 (normalized numbers $d_1\neq 0$)

• $d_2$: 10

• $E$: 3

• zero: 1

total:

What is the distribution on the real line?

What is the underflow and overflow limits?

• Smallest number that can be represented is the underflow: $1.0 \times 10^{-2} = 0.01$

• Largest number that can be represented is the overflow: $9.9 \times 10^0 = 9.9$

What is the smallest number $\epsilon_{mach}$ such that $1+\epsilon_{mach} > 1$?

• $\epsilon_{mach} = 0.1$

Binary Systems¶

Consider the 2-digit precision base 2 system:

Number and distribution of numbers¶

1. How many numbers can we represent with this system?

2. What is the distribution on the real line?

3. What is the underflow and overflow limits?

4. What is $\epsilon_{mach}$?

How many numbers can we represent with this system?

What is the distribution on the real line?

• Smallest number that can be represented is the underflow: $1.0 \times 2^{-1} = 0.5$

• Largest number that can be represented is the overflow: $1.1 \times 2^1 = 3$

• $\epsilon_{mach} = 0.1 = 2^{-1}= 1/2$

Note: these numbers are in a binary system.

Quick rule of thumb:

correspond to 8s, 4s, 2s, 1s . halves, quarters, eighths, ...

Real Systems - IEEE 754 Binary Floating Point Systems¶

Single Precision¶

• Total storage alloted is 32 bits

• Exponent is 8 bits $\Rightarrow E \in [-126, 128]$

• Fraction 23 bits ($p = 24$)

s EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF
0 1      8 9                     31

• Overflow $= 2^{128}\approx3.40\times10^{38}$

• Underflow $= 2^{-126} \approx 1.17 \times 10^{-38}$

• $\epsilon_{\text{machine}} = 2^{-23} \approx 1.19 \times 10^{-7}$

Double Precision¶

• Total storage alloted is 64 bits

• Exponent is 11 bits $\Rightarrow E \in [-1022, 1024]$

• Fraction 52 bits ($p = 53$)

s EEEEEEEEEE FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FF
0 1       11 12                                                      63

• Overflow $= 2^{1024} \approx 1.8 \times 10^{308}$

• Underflow $= 2^{-1022} \approx 2.2 \times 10^{-308}$

• $\epsilon_{\text{machine}} = 2^{-52} \approx 2.2 \times 10^{-16}$

Python Access to IEEE Numbers¶

Access many important parameters, such as machine epsilon:

import numpy
numpy.finfo(float).eps

Examples¶

Show that $(1 + \epsilon_{mach}) > 1$

Why should we care about this?¶

• Floating point arithmetic is not commutative or associative

• Floating point errors compound, do not assume even double precision is enough!

• Mixing precision can be dangerous

Example 1: Simple Arithmetic¶

Simple arithmetic $\delta < \epsilon_{\text{machine}}$.

Compare

Example 2: Catastrophic cancellation¶

Let us examine what happens when we add two numbers $x$ and $y$ where $x + y \neq 0$. We can actually estimate these bounds by doing some error analysis. Here we need to introduce the idea that each floating point operation introduces an error such that

where $\text{fl}(\cdot)$ is a function that returns the floating point representation of the expression enclosed, $\text{op}$ is some operation (e.g. $+, -, \times, /$), and $\delta$ is the floating point error due to $\text{op}$.

Back to our problem at hand. The floating point error due to addition is

Comparing this to the true solution using a relative error we have

so that if $\delta = \mathcal{O}(\epsilon_{\text{machine}})$ we are not too concerned.

What if instead we consider a floating point error on the representations of $x$ and $y$, $x \neq y$, and say $\delta_x$ and $\delta_y$ are the magnitude of the errors in their representation. Here we will assume this constitutes the floating point error rather than being associated with the operation itself.

Now consider the difference between the two numbers

Again computing the relative error we then have

The important distinction here is that now the error is dependent on the values of $x$ and $y$ and more importantly, their difference. Of particular concern is the relative size of $x - y$. As it approaches zero relative to the magnitudes of $x$ and $y$ the error could be arbitrarily large. This is known as catastrophic cancellation.

Example 3: Function Evaluation¶

Consider the function

with $x\in[-10^{-4}, 10^{-4}]$.

Taking the limit as $x \rightarrow 0$ we can see what behavior we would expect to see from evaluating this function:

What does floating point representation do?

Example 4: Evaluation of a Polynomial¶

Note: $f(1) = 0$ (and will be close to zero for $x\approx 1$)

Here we compare polynomial evaluation using naive powers compared to Horner’s method as implemented in eval_poly(p,x) defined above.

Example 5: Rational Function Evaluation¶

Compute $f(x) = x + 1$ by the function

Do you expect there to be issues?

Combination of Errors¶

In general we need to concern ourselves with the combination of both discretization error and floating point error.

Reminder:¶

• Discretization error: Errors arising from approximation of a function, truncation of a series...

• Floating-point Error: Errors arising from approximating real numbers with finite-precision numbers

or $\frac{1}{3} \approx 0.333333333$ in decimal cannot be represented exactly as a binary floating point number

Example 1¶

Consider a finite difference approximation to the first derivative of a function

Note: in the limit $\Delta x\rightarrow 0$, this is the standard definition of the first derivative. However we’re interested in the error for a finite $\Delta x$.

Moreover, (as we will see in future notebooks), there are many ways to approximate the first derivative. For example we can write the “centered first derivative” as

Here we will just compare the the error for the two different finite-difference formulas given

at $x=1$ for decreasing values of $\Delta x$. We will also introduce the idea of an ‘inline’ or lambda function in python.

Example 2¶

Evaluate $e^x$ with its Taylor series.

Can we pick $N < \infty$ that can approximate $e^x$ over a give range $x \in [a,b]$ such that the relative error $E$ satisfies $E < 8 \cdot \varepsilon_{\text{machine}}$?

We can try simply evaluating the Taylor polynomial directly for various $N$

And test this

Can we do better?¶

Note:

the largest value of $x$ such that $e^x <$ MAX is:

and numpy.exp handles that just fine

Your homework: the great Exp Challenge

One final example (optional): How to calculate Relative Error¶

Say we wanted to compute the relative error between two values $x$ and $y$ using $x$ as the normalizing value. Algebraically the forms

and

are equivalent. In finite precision what form might be expected to be more accurate and why?

Example based on a blog post by Nick Higham

Using this model the original definition contains two floating point operations such that

For the other formulation we have

If we assume that all $\text{op}$s have similar error magnitudes then we can simplify things by letting

To compare the two formulations we again use the relative error between the true relative error $e_i$ and our computed versions $E_i$.

Original definition:

Manipulated definition:

We see then that our floating point error will be dependent on the relative magnitude of $e$

Comparison of Relative Errors of estimates of Relative Error ;^)¶

Some other links that might be helpful regarding IEEE Floating Point:

• What Every Computer Scientist Should Know About Floating-Point Arithmetic

• IEEE 754 Floating Point Calculator

• Numerical Computing with IEEE Floating Point Arithmetic

Future issues with fp64 and High-Performance Computing¶

The Issues¶

• In traditional High-Performance computing IEEE fp64 has become the standard precision necessary for accurate, reproducible calculations for a wide range of scientific computing (e.g. climate models, fusion, solid mechanics)

• Until recently, the needs for HPC drove the development of Chips/Hardware such that Commodity Computers and Super Computers benefited from the same technology.

• However, with the rise of general purpose GPU’s and AI, the landscape is changing rapidly

A brief history¶

Dongarra et al., 2024

A brief history of floating point hardware¶

Dongarra et al., 2024

• 1980’s: Dedicated seperate floating point co-processors (e.g. Intel 8087, Motorola 68881 co-processors)

• 1987: Introduction of Intel x486 CPU’s with built in floating point registers

• 1999: Introduction of Nvidia GeForce 256 separate “Graphical Processing Unit” GPU low precision, fast parallel graphics.

• mid-2000s: Adoption of programmable General Purpose GPU’s for fp acceleration, addition of fp64 on GPU’s

• 2006: Introduction of Nvidia CUDA language for programmable GPU’s: shift to GPU’s for high-performance computing and ML/AI

• ~2020+: ML/AI revolution: Deep learning algorithms driven by matrix multiplications that tolerate low precision

Current Floating Point fp64 performance for CPU’s and GPU’s¶

Current Floating Point formats for CPU’s and GPU’s¶

Near Future NVIDIA GPU FloatingPoint roadmap StorageReview.com¶

Beyond Blackwell¶

Interesting times indeed¶

...the landscape of high-performance computation is increasingly complex...but there are important classes of problems that still need high-precision floating point. Some Options:

• fp64 Emulation leveraging low-precision hardware

• clever mixed precision algorithms

Operation Counting¶

Discretization Error: Why not use more terms in the Taylor series?

Floating Point Error: Why not use the highest precision possible?

Example 1: Matrix-Vector Multiplication¶

Let $A, B \in \mathbb{R}^{N \times N}$ and $x \in \mathbb{R}^N$.

1. Count the approximate number of operations it will take to compute $A x$.

2. Do the same for $A B$.

Matrix-vector product: Defining $[A]_i$ as the $i$th row of $A$ and $A_{ij}$ as the $i$, $j$th entry then

Take an explicit case, say $N = 3$, then the operation count is

This leads to 15 operations (6 additions and 9 multiplications).

Take another case, say $N = 4$, then the operation count is

This leads to 28 operations (12 additions and 16 multiplications).

Generalizing this there are $N^2$ multiplications and $N (N -1)$ additions for a total of

Matrix-Matrix product ($AB$): Defining $[B]_j$ as the $j$th column of $B$ then

The inner product of two vectors is represented by

leading to $\mathcal{O}(3N)$ operations. Since there are $N^2$ entries in the resulting matrix then we would have $\mathcal{O}(N^3)$ operations.

There are methods for performing matrix-matrix multiplication faster. In the following figure we see a collection of algorithms over time that have been able to bound the number of operations in certain circumstances. Here