# How to identify unstable floating point computations?

In numerics, it is very important to be able to identify unstable schemes and to improve their stability. How to identify unstable floating point computations?

I am working on a very complex simulation where many numerical schemes work together and I am looking for a methodic to identify its weak parts. I am working on a physical model involving differential equations. A bird's eye view of the overall process is:

1. (Preliminary step) Gather physical observations P.

2. Determine initial parameters of the simulation. This uses an optimisation algorithm, where we walk in a parameter space and look for parameters C such that some error function E(F(C), P) is minimised, where F is some derived quantity of the parameters.

3. Plug C in the simulation engine. This is an Euler scheme of the EDP, so that at each time step, we compute the terms driving the dynamic (each of them is a complex function, potentially subject to instability) and feed the Euler scheme with these dynamic terms to compute the next state. This goes on for thousands time points.

4. At the end of the simulation, we compute some function Proof(S) of the final state S and compare to some quantities Require(P) deduced from the observed quantities. This is not a formal proof of the result, more a plausibility check.

Also, I see a tower of complex operations (computation of dynamic terms, within the Euler scheme, within the Proof). And would like to recognise “bad parts” and fix them.

I speculate that using a software implementation of floating point numbers with reduced precision would magnify the unstability of numerical schemes, thus easing the comparison between different implementations. Is this a common technique to investigate this question? Is it possible to use a virtual machine, as Bochs, to achieve this without altering the program?

To deal appropriately with the stability question, it is sometimes acceptable to target the typical input of the numerical procedure, so that it can be tuned to do well on that input and maybe less well on other valid, yet unlikely, input. Given a sample of typical inputs, it is possible to snoop some intermediate results and prepare a statistical profile for them. Again, is this a common technique to study stability issues? Is a virtual machine useful for this?

• may be you will have more interesting answers on math.stackexchange.com – Simon Bergot Dec 11 '13 at 9:54
• @Simon You might be right, but this is definitely a cross-domain question. I guess persons able to answer are wether registered to both math and programmers or to none… Let's wait a bit an see if this question finds its answer here! – user40989 Dec 11 '13 at 11:33
• Interval arithmetics? – SK-logic Dec 11 '13 at 21:12
• – SK-logic Dec 12 '13 at 9:50
• Using Euler to propagate state isn't necessarily evil; neither is optimization, but you really have to split the problem into subtasks. Numerical instability may be the least of your woes - convergence to a false maximum and problems related to the stiffness of the ODEs/PDEs loom larger than that. And yes, never ever use single precision :) – Deer Hunter Dec 20 '13 at 11:21

The study of stability of floating point computation is part of numerical analysis and if you really want a sound result, you really want someone knowledgeable in that domain to do the analysis of the algorithms used.

There are some things which can help to experimentally identify unstable algorithms. Running with rounding set to different modes (up/down/random) or with different precision and checking that the result don't vary too much. Answering is this too much? isn't simple at all and even when the answer is no that doesn't mean the algorithm is stable, just that it wasn't detected unstable on the data set you used.

Interval arithmetic has been proposed in comments. When I looked at it, even the most rabid proponent of interval arithmetic admitted that it worked well with algorithms designed for interval arithmetic but that switching to it without analyzing the algorithm and ensuring that it hadn't patterns known not to work well wouldn't be useful (opponents seemed of the opinion that the pre-conditions for interval arithmetic to be useful where too restrictive to be of practical interest)

Designing stable floating point algorithms is highly nontrivial. Those more mathematically skilled than myself suggest using well regarded libraries where possible rather than attempting to roll your own. The standard reference in the area seems to be:

N. J. Higham. Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second edition, 2002. ISBN 0-89871-521-0

Not knowing more about what types of computations, languages, etc. makes it hard to give much of a concrete answer. There is a good lecture here: http://introcs.cs.princeton.edu/java/91float/ this may be a little basic, but it's a good introduction if you are using java.

How to identify unstable floating point computations? Is this a common technique to investigate this question?

I think unless you need to show some statistics on error, you don't really need to gather samples. What you need is Numerical Analysis, which also falls under the subjects of Numerical Methods, Numerical Linear Algebra, etc. And they're part of computer science, so you might get some answers too in cs.stackexchange.

Anyhow, in general programming, most problems are easy to spot given some basic understanding about how floating point works and basic numerical methods. But even more complex problem are "easier" to solve today with availability of 128-bit floats, even less reason to produce error samples. Here are some example problems to show my point:

1. using floating point for calculating monetary values.
2. using floating point for big numbers.
3. not doing divisions before other operations when it's possible to do so. (to make value closer to 0).
4. long calculation without special treatment for error propagation.

There's also example of a naive algorithm and error compensated algorithm here algorithm for calculating variance. In the example, looking at the naive version, you just kind of can smell that doing calculation in loops will carry some errors and is not being compensated.

• Thank you for your answer, I am however looking for more detailed information. I have a very large computation and want to identifiy its weak parts. I edited the question accordingly. – user40989 Dec 12 '13 at 9:54
• I'm not exactly sure what's your situation when you say you have large computation and want to identify weak parts. Numeric calculations inherently have errors in them, even one simple add operation. So, unless your large computation is error compensated, then they're as a whole need to be fixed. Improving weaker spots might not be good enough. If you now the "epsilon" of your floating point model, a simple analysis will show how big the error can be as they propagate through the long computation. – imel96 Dec 14 '13 at 22:33

You can avoid numerical errors by using appropriate data types (like, for example, continued fractions). If you need or want to use floating point arithmetic, you need to apply numerical know-how to know the errors.

• I do not want to avoid numerical errors, I want to find which parts of the computation are unstable. It is similar to localising speed bottlenecks when optimising for speed. So I want to optimise precision and therefore wants to find precision bottlenecks. (Continued fractions are not useful here.) – user40989 Dec 12 '13 at 11:32
• @user40989, then you definitely need interval arithmetics. – SK-logic Dec 12 '13 at 11:44