Roundoff error and finite precision representations - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Introduction to finite precision¶

Round-off errors are an unavoidable consequence of representing numbers with a finite number of digits. This affects both irrational numbers and rational numbers that cannot be perfectly represented in a given base.

Because memory is finite, computers must round or truncate numbers. While a single rounding error is tiny, these small errors can accumulate over millions of computations and lead to massive discrepancies.

The average human is capable of around one mistake per second, but computers can make millions of mistakes a second!

The magnitude of round-off error is influenced by the precision (the number of digits used) and the base of the number system. The combination of precision and base sets the finite bounds of numbers we are able to represent.

print(f"Pi: {np.pi}")
print(f"e: {np.e}")
print(f"1/3: {1/3}")
print(f"sqrt(2): {np.sqrt(2)}")

Binary Representation (Base-2)¶

Humans typically use a base-10 (decimal) numbering system. Computers use a base-2 (binary) system, where each bit is either 0 or 1.

An exception to this is the Mayans, who used a base-20 (vigesimal) system.

Binary integers¶

In the decimal system, 1305 is expressed as:

$1305_{10} = 5 \times 10^0 + 0 \times 10^1 + 3 \times 10^2 + 1 \times 10^3$

In binary, each digit’s position corresponds to a power of 2:

$10100011001_2 = 1 \times 2^0 + 0 \times 2^1 + 0 \times 2^2 + 1 \times 2^3 + 1 \times 2^4 + 0 \times 2^5 + 0 \times 2^6 + 0 \times 2^7 + 1 \times 2^8 + 0 \times 2^9 + 1 \times 2^{10} = 1305_{10}$

print(np.binary_repr(1305))

Binary fractions¶

Binary representation can also include fractional parts using powers of 2 (e.g. 2^-1, 2^-2). For example, 54.75₁₀ can be written in binary.

$110110.11_2 = (1 \times 2^5 + 1 \times 2^4 + 0 \times 2^3 + 1 \times 2^2 + 1 \times 2^1 + 0 \times 2^0) + (1 \times 2^{-1} + 1 \times 2^{-2}) = 54.75_{10}$

print(f"54.75: {exe.decimal_to_binary(54.75)}")

Unexpectedly non-terminating fractions¶

Some numbers that are terminating decimals in one base aren’t in binary!

Why is this important?

print(f"0.1: {exe.decimal_to_binary(0.1)}")

Integer Representation (Exact but Bounded)¶

Integers are used for signed numbers and do not suffer from round-off error. However, they have a limited absolute range bounded by the number of bits used to store them.

The range is bounded by the maximum value representable: $\text{base}^{\text{precision}}$ and split between positive and negative values. Due to the representation of zero, the negative range is slightly larger.

Negatives can be generally be represented in different ways:

Use a ‘sign bit’ (e.g.: 0=positive, 1 = negative)
Apply a bias (half the full precision range)
The Two’s Complement method which works specifically in binary arithmatic.

Integer Bounds¶

The minimum and maximum values for a signed integers in binary are:

$min = -2^{\text{bits}-1}$

$max = 2^{\text{bits}-1} - 1$

What is the largest integer a double-precision (64-bit) variable can store? We can check with the built-in numpy examiner:

print(f"min: {-2**63}")
print(f"max: {2**63-1}")
print()
print(np.iinfo(np.int64))

Overflow error¶

What happens when you store a number too large? It will cause an overflow error. Let’s try some:

# works because 2**62 is within the range of a 64-bit integer
print(np.int64(2**62))

# This will cause an overflow error because 2**63 is too large
try:
    print(np.int64(2**63))
except OverflowError as e:
    print(f"OverflowError: {e}")

# This also works
print(np.int64(1000000000000000000))

# This will also cause an overflow error
try:
    print(np.int64(10000000000000000000))
except OverflowError as e:
    print(f"OverflowError: {e}")

Integer Operations¶

Some operations between data types can give unintuitive answers. What do you think the answers to these operations are?

a = np.int64(4)
b = np.int64(3)

print(f"a + b = {a+b}")
print(f"a - b = {a-b}")
print(f"a * b = {a*b}")
print(f"a / b = {a/b}")

Surprised?¶

Were you expecting $a/b$ to result in 1? This is a case of package tools hiding numerical complications. In this case, python is sophisticated enough to know that division of integers aren’t always integers, but other languages optimized for speed might truncate it to 1. Beware!

Floating-Point Numbers (Approximate but Wide Range)¶

Writing out very large or small numbers is impractical. It is much more efficient to use scientific notation to represent the magnitude as an exponent: $10,000,000,000,000,000,000 = 10^{19}$

Scientific Notation (Base-10)¶

We can remove placeholder zeros by using a floating point to separate the fractional part (mantissa) from the order of magnitude (exponent).

Scientific Notation: $mantissa \times 10^{exponent}$

Decimal	Scientific Notation	Mantissa	Exponent
265.73	$2.6573 \times 10^2$	2.6573	2
0.0001	$1 \times 10^{-4}$	1	-4
-0.0034123	$-3.4123 \times 10^{-3}$	-3.4123	-3
$1500^*$	$1.5 \times 10^3$	1.5	3

*Assuming the trailing zeros are not significant.

Note:

The mantissa is a fraction. If we could assume the decimal point, it becomes an integer. (Why is this important?)
The exponent is the power of the number system’s base (in this case, 10).

Binary floating-point numbers¶

The same floating-point concept can be applied to binary numbers using base 2: $mantissa \times 2^{exponent}$

To normalize $54.75_{10} = 110110.11_2$ , we move the binary point so that it is after the first non-zero digit: $1.1011011_2 \times 2^5$

The exponent is 5, which in binary is 101₂. So the full floating-point representation is: $1.1011011_2 \times 2^{101_2}$

In binary, the mantissa always starts with a 1. Since this is true we can drop it from the representation entirely and represent the fractional part as an integer (why would we do this?)

IEEE 754 Standard and Precision for binary floating point¶

Since we are limited by a finite number of bits, floating-point numbers are divided into precision parts according to the IEEE 754 standard:

Precision	Total Bits	Sign	Exponent	Mantissa
Single	32	1	8	23
Double	64	1	11	52
Quad	128	1	15	112

The sign bit determines if the number is positive or negative. The exponent bits store the magnitude (shifted by the bias to account for negatives), and the mantissa stores the fractional part (following the leading 1).

Python by default using double precision floats.

# Look at how numbers are stored internally using our helper function
print(f"54.75 in binary fraction is, {exe.decimal_to_binary(54.75)} but in memory is: {exe.python_internal_binary(54.75)}")
# Sign = 0 (positive), exponent = 5+ (2^10-1) = 1028 = 10000000100, mantissa after leading 1 = 1011011 + trailing zeros.
print(f"-54.75 in memory: {exe.python_internal_binary(-54.75)}")

Sources of Round-off Error¶

Roundoff error is the difference between the true number and the finite-precision representation. This error can be mitigated systematically by using higher precision at the cost of computational speed and memory.

Although individual errors are small, algorithms involve many steps and errors can easily accumulate. Compilers and simplification steps may mitigate these errors, but it is often best to rely on packaged subroutines that are numerically stabilized - specifically to avoid such problems.

NB: “Numerical stability” has two meanings in numerical methods. The first is robustness to roundoff error, the second (discussed with PDEs) is concerned with propogation of errors.

Unrepresentable Numbers¶

Some numbers cannot be represented cleanly in base-2. The binary representation of 0.1 is actually a repeating fraction, which means it cannot be perfectly represented with a finite number of bits.

print(f"0.1 in binary is: {exe.decimal_to_binary(0.1)}")
print(f"0.1 in memory:    {exe.python_internal_binary(0.1)}\n")

print("Because 0.1 is a repeating fraction in binary, it is not perfectly 0.1 internally:")
print(format(0.1, '.55f'))

Disparate Magnitudes (Loss of Associativity)¶

Finite precision means small numbers can be completely lost when added to very large numbers. Because of this, we cannot always rely on the associative property of addition.

print(f"-1 + (1 + 1e-20) = {-1+(1+1e-20)}")
print(f"(-1 + 1) + 1e-20 = {(-1+1)+1e-20}")

# Observe what happens here:
print(f"-1 + 1 + 1e-20 = {-1 + 1 + 1e-20}")
print(f"-1 + 1e-20 + 1 = {-1 + 1e-20 + 1}")

Subtractive Cancellation¶

Subtractive cancellation happens when two nearly equal numbers are subtracted, causing a dramatic loss of significant digits.

A common task is to find the rate of change of a measurement. This is discussed later in finite difference but naively one would expect: $\Delta T = T_{final} - T_{initial}$ which generates spurious results when $T_{final} \approx T_{initial}$ .

What does this mean for data collection and analysis?

a = np.float32(1.23456789)
b = np.float32(1.23456780)
res = a - b

print(f"a = {a:.20f}")
print(f"b = {b:.20f}")
print(f"a - b = {res}")
print(f"Expected ~0.00000009")

Machine epsilon¶

Whereas subtractive cancellation exposes the binary ‘noise’, at some point a computer just isn’t able to distinguish between numbers with too many significant digits. Machine epsilon is defined as the smallest increment you can add to 1.0 that results in a number that is different from 1.0.

It is closely related to the mantissa (which captures the significant digits in the binary representation). For double precision, this is approximately $2^{52} = 10^{-16}$ (not exactly since the mantissa is in binary with teh leading 1 dropped.)

print(2**52)
[print(i, ":", 1.0+10**i) for i in range(-10,-20,-1)];

Example: The Quadratic Formula¶

Recall the quadratic formula: $x = \frac{-b\pm \sqrt{b^2-4ac}}{2a}$

If $4ac \ll b^2$ this will result in subtractive cancellation in the numerator. Let’s solve the roots of: $x^2 + 10^8x + 1$

a, b, c = 1, 1e8, 1
desc = np.sqrt(b**2 - 4*a*c)

print('Naive roots are: ', (-b+desc)/(2*a), 'and', (-b-desc)/(2*a))

Numerically Stable Implementation¶

The problem occurs when $b \gg ac$ . The discriminant becomes $\approx b$ and the numerator can suffer subtractive cancellation on either the plus or negative, depending on the sign of $b$ .

A numerically stabilized method accounts for the sign of $b$ ; $sign(b)= \{1, -1\}$ , if $b$ is positive or negative, and uses the relationship $x_1 \cdot x_2 = c/a$ .

q = -0.5 * (b + np.sign(b) * np.sqrt(b**2 - 4*a*c))

print('Stable roots are: ', q/a, 'and', c/q)