Nicholas Higham's tips for designing stable numerical algorithms

I was looking through Nicholas Higham’s Accuracy and Stability of Numerical Algorithms and was pleasantly surprised by an informative summary of general tips to bear in mind when designing numerical algorithms or writing programs that depend on numerical algorithms.

Here they are:

1) Try to avoid subtracting quantities contaminated by error.

The problem of cancellation and the many woes of subtraction are a recurring theme in his advice. Although subtractions may be unavoidable, you are best off trying to avoid them unless you are doing integer (whole number) arithmetic where the input data and output solutions are exact.

2) Minimize the size of intermediate quantities relative to the final solution

Higham’s reasoning: If intermediate quantities are very large (relative to the final solution) then the final answer may be the result of “damaging subtractive cancellation”. In other words, large intermediate numbers swamp the initial data, resulting in loss of information.

Since we don’t have mathjax or any kind of LaTeX plugin for discourse, it is difficult for me to describe this more rigorously, but essentially bigger numbers mean bigger absolute errors (as opposed to relative errors which are scale invariant) and, as stated in tip #1, subtraction can be very very mean and error prone, especially when numbers are close together.

3) Look for different formulations of a computation that are mathematically but not numerically equivalent.

Higham’s example: f(x) = (1 - cos(x)) / x^2 with x = 1.2 x 10^-5 and cos x rounded to 10 significant figures.

Software that isn’t smart enough to detect pathological math problems will screw this problem up. f(x) cannot be greater than or equal to 0.5 or less than 0, but plugging this problem in naively may yield you nonsense answers depending on the tech you are using. Instead, you can use the equivalent version

f(x) = (1/2)[sin(x/2)/(x/2)]^2 , ( this comes from the identity cos(x) = 1 - 2sin^2(x/2) )

This dodges the pathological behavior of the earlier f(x) because the issue with the earlier f(x) is that cos x = 0.999999999999 (rounded to ten significant figures) and you are subtracting it from 1. This type of computational error is an example of the cancellation phenomenon. The cancellation phenomenon occurs when two nearly equal numbers are subtracted from one-another.

4) Use only well-conditioned transformations of the problem. In matrix computations this means multiplying by orthogonal matrices instead of nonorthogonal matrices when possible.

Here’s what this means:

In numerical computing, there is a value known as the condition number.

For those who don’t yet have calculus down, the condition number essentially shows the sensitivity of a solution to perturbations (errors, uncertainties, etc.) in the input data. Data almost always has some error or uncertainties of some kind, so we say that a function or function-equivalent (like a matrix) with a very high condition number is ill-conditioned since small errors in the input data can greatly magnify errors in the output values. When the condition number is relatively small, we say that the approximation problem is well-conditioned. What counts as well-conditioned or ill-conditioned depends on the problem since different kinds of problems require their solutions to have different levels of accuracy.

For the folks who are good at calculus, the condition number is motivated by a first-order taylor approximation to the relative error between an approximation and a true value. Assuming the function is at least twice continuously differentiable, you find the first order taylor polynomial of f(x+h), substract f(x) from both sides to get f(x+h) - f(x) = f’(x)(x+h) + R1, where R1 is the remainder term for the first order taylor approximation. You find R1 using the Taylor Remainder Theorem and then apply the Extreme Value Theorem to show that there exists a maximum value z on the interval [x, x+h] such f’’(z) is the maximum value of f on the interval [f(x), f(x+h)]. You can then use big-O notation to show that

approximated_value - true_value = [(x * f’(x))/ f(x)] * (x+h)/(x) + O((x+h)^2)

The [(x * f’(x)) / f(x)] term is the relative condition number of f.

5) Take precautions to avoid unnecessary overflow and underflow.

My favorite example: Do you find yourself doing a lot of multiplication and division in your programs? Multiplication in particular can sneak up on you and cause serious overflow and underflow. Use logarithms to convert the multiplication into addition. You can undo logarithmic transformations with e^(ln(x)) = x. Logarithms are the most overpowered thing that you learned in school, but school probably never showed you just how much easier your life becomes when you use logarithms. Multiplication and division are a pain in the ass for both math folks and non-math folks; always transform the problem into addition and multiplication using logarithms except when the problem is trivial or the computation is short. The whole reason logarithms were invented in the renaissance was because it is much, much, much, much, much faster to do addition and subtraction rather than multiplication and division. The industrial revolution and the entire world of comfort we live in is entirely owed to logarithms. You have no idea how much your life relies on logarithms. You should build an altar to logarithms in your house and have your children memorize all of its properties.

1 Like

One of my electrical engineering professors specifically emphasized the danger with “small differences between large numbers” and the need to restructure calculations to avoid them.

1 Like

Oh yes, very! Did you know that solving the quadratic equation (numerically) is a highly non-trivial computer science problem that we didn’t have a reliable answer to until the late 20th century?

The easiest issue there is to deal with the case where b is much greater than |4ac| so that sqrt(b^2 - 4ac) is roughly equal to |b|. In this case, you get cancellation errors and the trick to getting around this issue to only solve for the larger of the two roots and then use the larger root to algebraically solve for the other root using the fact that x1*x2 = c/a, where x1 and x2 are the roots of the quadratic equation.

Within the radical, there is an even worse problem that you actually cannot fix through numerical tricks. To account for the situation in which b^2 is roughly 4ac (meaning the roots of the quadratic equation are close together), you have to increase the machine precision or use some trick that amounts to extending the precision of your machine.

The most non-trivial part of reliably solving the quadratic equation comes with trying to account for potential underflow and overflow. You have to implement a clever scaling algorithm to account for this situation.

1 Like