In the previous post we saw that that floating point arithmetic behaves very differently from arithmetic with real numbers. The following Python code illustrates a problem called catastrophic cancellation:
x = 3.141592653589793
y = 6.022e23
x = (y + x) - y
print(x)
The correct answer is clearly pi but Python outputs 0 because the term (y + x) has the same floating point representation as y, and we cancel these identical terms the result is zero.
In Part 1 we saw that we nearly always have small errors in numerical calculations due to the finite precision of floating-point arithmetic. However, the error above is clearly not a small error! But what exactly do we mean by “small”?
In numerical calculations one way to quantify error is using absolute error. The absolute error is simply the absolute difference between the true result, x, and the calculated result, f:
In most of our calculations there will nearly always be a non-zero absolute error, which is why we have to be careful using the equality error on floating-point values.
If there are always errors, how can we tell whether they are large enough to matter? The answer is to look at the relative error:
In the case of our Python example above, x = pi and f = 0, and our relative error is 1, or as a percentage 100%. The art of numerical computing is to ensure that this relative error is very small, and this principle is clearly violated in the case of catastrophic cancellation.
In Python, we can check whether a result is within a relative error using numpy’s isclose
function, which takes an optional parameter rtol
specifying the relative tolerance. By default, the relative tolerance is 1e-5, and we can check the result of our catastrophic cancellation like this:
import numpy as np
x = 3.141592653589793
y = 6.022e23
x = (y + x) - y
f = 0
print(np.isclose(f, x))
The above program prints out False
indicating that our result is not within the desired relative error, which is what we expect.
But let’s see what happens if we increase the tolerance to 100% (we would never do this in practice, but let’s see what happens for the purposes of experimentation).
import numpy as np
x = 3.141592653589793
y = 6.022e23
x = (y + x) - y
f = 0
print(np.isclose(f, x, rtol=1.0))
Now we see that the result is True, as the relative error is within the specified tolerance of 1.
So in general in our numerical calculations we use isclose to test for equality. But have you spotted something strange about the definition of R? Let’s see what happens if we swap the arguments around when we call isclose
:
print(np.isclose(x, f, rtol=1.0))
Now the answer is False
! We can see why by looking at definition of R in the equation above. The denominator is x, which is the target value we are trying to represent, and if we swap x and f we will obtain a different result. Therefore, when we call np.isclose
, we always have to put the target value as the first argument, and therefore our definition of floating-point equality is not commutative. And after some careful thought you will also realise that floating-point addition and subtraction are not associative. Welcome to the topsy-turvy world of floating point arithmetic!
All of this suggests we have to be extremely careful when writing numerical programs. We cannot mindlessly take mathematical equations from a text-book and translate them into Python, and then expect our program to give the correct results.
How then can we avoid errors like catastrophic cancellation in our code? There is no simple answer to this question, just as there is no silver bullet for writing bug-free code in general (whether it is numerical or not), but the first step in ensuring that any code is free of errors is to ensure that is thoroughly tested.
Professional software engineers use a technique called Test-Driven Development (TDD) to help reduce the amount of bugs in their code. When we use TDD, before writing any code we write a test to check that our code is correct. Moreover, we write down this test in code, so that the test can be quickly run whenever we make any changes to the code.
Of course, at first, our test will fail because we have not even written the actual code to be tested! However, once the test is in place, we can easily check our program for errors as we start writing the code. This continuous feedback loop ensures that we are always building on a solid foundation as our program grows in size and complexity.
This iterative approach, where we continually and incrementally develop our system in response to feedback, is the key principle underpinning modern software engineering, and there is considerable evidence that it leads to more reliable software. One of the key pioneers of this approach to software engineering is Dave Farley, and you can learn more at his blog, and book, linked below.
In the next post in this series we will use Test-Driven Development (TDD) to write a Python function to calculate the sample variance. As we shall see, calculating the sample variance accurately is not as straightforward as it may first seem…
Paid subscribers can download a Jupyter notebook with explanations and examples of numerical computing in Python at the link below.
Keep reading with a 7-day free trial
Subscribe to The Life Algorithmic to keep reading this post and get 7 days of free access to the full post archives.