The choice of algorithm for solving a problem can have a dramatic effect on the accuracy.

Introduction

Perhaps the only formula we are likely to remember from our school days is that for the roots of a quadratic equation: if b² ≥ 4ac then the function q(x) = ax² + bx + c is zero at x = α and x = β, where

α =(-b + √(b² - 4ac))/2a,

and

β =(-b - √(b² - 4ac))/2a.

The quadratic can be factored as

q(x) = a(x - α)(x - β)

from which we see that

aαβ = c.

Using a calculator or programming a computer to calculate the roots α and β would seem a straightforward task.

The Calculation

However, even for this computation, the way we organise the calculation, i.e., the algorithm we employ, can have a dramatic effect on the accuracy. Our initial approach may be to follow the formulae above: first calculate d = √(b² - 4ac) and then set α = (-b + d)/(2a), β = (-b - d)/(2a).

Problems arise when b is large in magnitude relative to 4ac for in this situation d will be close to either b or -b and subtractive cancellation error will occur in calculating one of the roots. For example if a = 1, b = 1000 and c = 1, then b and d are the same to the first five digits so that there are five less digits available to store their difference d - b. In this case β is calculated accurately (there is no subtractive cancellation error) but approximately three digits of accuracy are lost in calculating α. A better algorithm is to calculate the root larger in magnitude first and then use the relationship aαβ = c to calculate the other. So if b is positive we calculate β =(-b - (d)/(2a) and set α = c/(βa). Otherwise, we calculate α first and then β. Using this approach both roots are calculated accurately. (If we chose always to calculate β then α, we would lose approximately six digits accuracy α for the case a = 1, b = -1000 and c = 1, for example.)

This simple example illustrates a number of points:

• numerical instability can occur in the simplest of calculations;
• an algorithm stable for one set of inputs can be unstable for another;
• it is not always obvious that an algorithm is likely to be unstable.

Modern Numerical Software

If an algorithm is discovered to be unstable, then we have the problem of how to replace it by a more stable one. In the early days of numerical computing the difficulties in analysing the numerical properties of algorithms led to a pessimistic view about achieving reliable computation in finite precision arithmetic. Pioneering work by Jim Wilkinson at NPL and others during the 1960s and 1970s in the development of new error analysis methods made the problem much more tractable. Today we have access to a large number of numerical software routines, particularly for numerical linear algebra, that implement algorithms with proven stability properties.

This Case Study is an abridgement of an article that first appeared in Counting on IT Issue 11.