Accuracy and stability of numerical algorithms

What Hacker News says

1. cmrx64Jul 2025
   Displeased with this trend of LLM-assisted “research”. The central claims are false, the examples of floating point are false and off by a factor of 1e11, the latency numbers are complete WTF and disconnected from reality (‘cycles’ of what).

   Fixed point arithmetic with a dynamic scale is presented along the way motivating floating point in probably every computer architecture class. It’s a floating point.

   This guy needs to open a book. I recommend Nick Higham’s _Accuracy and Stability of Numerical Algorithms_.

   Show morethread
2. dahartNov 2024
   > More generally, relative error works only if your computation scales “multiplicatively” from zero; if there’s any additive component, it’s suspect.

   IEEE floating point inherently scales from zero, the absolute error in any computation is proportional to the magnitude of the input numbers, whether you’re adding or multiplying or doing something else. It’s the reason that subtracting two large numbers has a higher relative error than subtracting two small numbers, c.f. catastrophic cancellation.

   > How does one do this properly?

   There’s a little bit of an art to it, but you can start by noting the actual result of any accurate operation has a maximum error of 0.5 LSB (least significant bit) simply as a byproduct of having to store the result in 32 or 64 bits; essentially just think about every single math instruction being required to round the result so it can fit into a register. Now write an expression for your operation in terms of perfect math. If I’m adding two numbers it will look like a[+]b = (a + b)*(1 + e), where e is your 0.5 LSB epsilon value. For 32 bit float, e == +/- 1e^-24. In this case I differentiate between digital addition with a finite precision result, and perfect addition, using [+] for digital addition and + for perfect addition.

   This gets hairy and you need more tricks for anything complicated, but multiplying each and every operation by (1+e) is the first step. It quickly becomes apparent that the maximum error is bounded by |e| * (|a|+|b|) for addition or |e| * (|a| * |b|) for multiply… substitute whatever your operation is.

   When doing more complicated order-dependent error analysis, it’s helpful to use bounds and to allow error estimates to grow slightly in order to simplify expressions. This way you can prove the error is less than a certain expression, but the expression might be conservative.

   A 3d dot product is a good example to work though using (1+e). Typically it’s reasonable to drop e^2 terms, even though it will technically compromise your error bound proof by some minuscule amount.

       a[*]x [+] b[*]y [+] c[*]z = ((ax(1+e) + bx(1+e)) + cx(1+e))(1+e)
       = ((ax+axe + by+bye)(1+e) + cz+cze)(1+e)
       = (ax+by+(ax+by)e + axe+bye+(ax+by)ee + cz+cze)(1+e)
       = (ax + bx + 2e(ax+by) + e^2(ax+by) + cz+cze)(1+e)
       = ax + by + 2e(ax+by) + e^2(ax+by) + cz + cze + axe + bye + 2e^2(ax+by) + e^3(ax+by) + cze + cze^2
       = ax+by+cz + 3e(ax+by) + 2e(cz) + 3e^2(ax+by) + e^2cz + e^3(ax+by)
   

   Now drop all the higher order terms of e.

       = ax+by+cz + 3e(ax+by) + 2e(cz)
   

   Now also notice that 2e|cz| <= 3e|cz|, so we can say the total error bound:

       <= (ax + by + cz) + 3e( |a||x| + |b||y| + |c||z| )
   

   And despite the intermediate mess, this suddenly looks very conceptually simple and doesn’t depend on the order of operations. If the input values are all positive, then we can say the error is proportional to 3 times the magnitude of the dot product. And it’s logical too because we stacked 3 math operations, one multiply for each element of the sum and two adds.

   Sorry if that was way too much detail… I got carried away. :P I glossed over some topics, and there could be mistakes but that’s the gist. I’ve had to do this for my work on a much more complicated example, and it took a few tries. There is a good linear algebra book about this, I think called Accuracy and Stability of Numerical Algorithms (Nicholas Higham). The famous PBR graphics book by Pharr et. al. also talks about error estimation techniques.

   Show morethread
3. bee_riderOct 2022
   He also wrote "Accuracy and Stability of Numerical Algorithms" which is part of the canon of the field.
   thread
4. ChrisLomontJun 2022
   >run it twice, once rounding towards negative infinity, and once towards positive infinity, and see how different they are

   That's actually a good method that lots of numerical analysts use. I think C99 has several rounding modes - I often have code that runs the same computation under each mode and look at what your doing.

   >not all crappy code generates a wide window

   True, often you need to find the values that make it puke. I think there was a paper recently on doing this automatically.... If I can refind it I'll post - it was cool stuff.

   >students are only allowed to use single (32-bit) precision floats and operations

   For stuff this small you can often brute force all values or at least all edge case values. I sometimes do the trick of converting 32 bit integers (watch out for weird C aliasing rules) to hit lots of numbers near precision edges and run them through, denormals, other edge cases, and pseudo-fuzz the results. This will sometimes turn up more.

   If you teach, and have not worked through some numerical analysis books that provide good background for this stuff, read Goldberg's paper "What Every Computer Scientist Should Know About Floating-Point Arithmetic" [1]. The next good thing to dig through is Higham's book "Accuracy and Stability of Numerical Algorithms" [2]

   Those will give you lots of places to look for issues, and the ability to analyze algorithms more concretely for floating point issues.

   [1] https://docs.oracle.com/cd/E19957-01/800-7895/800-7895.pdf [2] http://ftp.demec.ufpr.br/CFD/bibliografia/Higham_2002_Accura...

   Show morethread
5. conistonwaterOct 2014
   If you are interested in writing accurate numerical algorithms, as further reading I highly recommend reading "Accuracy and Stability of Numerical Algorithms" by Nick Higham, which is a very good comprehensive book. Like this article says, floating-point arithmetic is often thought of as mysterious, but it really isn't: it just obeys its own precisely specified rules. I think it's very good to dispel the mystery, so I really like this article.

   Another way to do the same type of error analysis of numerical algorithms is to use a parameter sensitivity library. In the end, numerical stability can be thought of as the sensitivity of the final answer to separate independent changes in each intermediate value obtained in the computation of the answer. This is sometimes overlooked, I think, and it's quite easy to implement as well.

   I have one quibble. They say "Use (x-y)(x+y) instead of x^2-y^2". This is true, but (in my opinion) misleading: the function (x,y) -> x^2-y^2 is itself ill-conditioned when x is near y, and using a more numerically stable algorithm to evaluate an ill-conditioned function will not improve the function's condition. If x and y are close to each other, the inaccuracy introduced by using x^2-y^2 is close to the inaccuracy due to using approximate values for x and y. So it you have x^2-y^2, I think the first thing you should worry about is whether your computation is ill-conditioned, and only then worry whether it's numerically stable. I say this because I sometimes see people try to evaluate x^2-1, notice it is inaccurate for x near 1, and try to replace x^2-1 with (x-1)*(x+1), which is still inaccurate.

   Show morethread
6. delhantyMar 2017
   Floating point is Nick Higham's thing isn't? I have his book "Accuracy and Stability of Numerical Algorithms, Second Edition"

   http://www.ma.man.ac.uk/~higham/asna/

   but I haven't got very far with it yet.

   Maybe the article in the Princeton Companion is a summarized version of his book?

   thread