18

I am writing a physics simulator in C++ and am concerned about robustness. I've read that catastrophic cancellation can occur in floating point arithmetic when the difference of two numbers of almost equal magnitude is calculated. It occurred to me that this may happen in the simulator when the dot product of two almost orthogonal vectors is calculated. However, the references I have looked at only discuss solving the problem by rewriting the equation concerned (eg the quadratic formula can be rewritten to eliminate the problem) - but this doesn't seem to apply when calculating a dot product? I guess I'd be interested to know if this is typically an issue in physics engines and how it is addressed.

Adam Zalcman
  • 26,643
  • 4
  • 71
  • 92
Nick Kovac
  • 487
  • 3
  • 10
  • 1
    maybe something for [math.stackexchange](http://math.stackexchange.com/)? – Tom Knapen Jan 15 '12 at 11:11
  • 1
    @TomKnapen - While the question is probably in scope for math.stackexchange, the [faq at that site](http://math.stackexchange.com/faq) suggests using SO for questions about algorithm implementation. This question belongs here. – Ted Hopp Oct 30 '12 at 05:59
  • For a lot of applications, it doesn't matter. If you get catastrophic cancellation, the ratio of the correct dot product and the calculated one may be significantly different from one ... but both values will be very close to zero, and compared to other values you are calculating, will be insignificant. – Martin Bonner supports Monica Feb 02 '17 at 21:24

2 Answers2

13

One common trick is to make the accumulator variable be a type with higher precision than the vectors itself.

Alternatively, one can use Kahan summation when summing the terms.

Another approach is to use various blocked dot product algorithms instead of the canonical algorithm.

One can of course combine both the above approaches.

Note that the above is wrt general error behavior for dot products, not specifically catastrophic cancellation.

janneb
  • 36,249
  • 2
  • 81
  • 97
  • 1
    Thanks for your answer... it seems that those solutions address round-off errors rather than cancellation though. In my example, I am calculating a dot product of 2d vectors, x1*x2 + y1*y2 (where the variables are floats). I don't think converting them to doubles would solve the issue of cancellation (though it might make it less frequent). If x1*x2 is almost the same magnitude as y1*y2, then I also don't see how the Kahan algorithm would address the cancellation problem that would occur... it seems more to do with minimising error accumulation in multiple sums. – Nick Kovac Jan 15 '12 at 11:52
  • https://stackoverflow.com/a/63665011/768469 has a good overview of the general problem, including references and a different algorithm from Kahan (assuming your hardware provides FMA). – Nemo Jun 04 '23 at 23:32
5

You say in a comment that you have to calculate x1*x2 + y1*y2, where all variables are floats. So if you do the calculation in double-precision, you lose no accuracy at all, because double-precision has more than twice as many bits of precision as float (assuming your target uses IEEE-754).

Specifically: let xx, yy be the real numbers represented by the float variables x, y. Let xxyy be their product, and let xy be the result of the double-precision multiplication x * y. Then in all cases, xxyy is the real number represented by xy.

TonyK
  • 16,761
  • 4
  • 37
  • 72
  • 1
    +1: This sounds correct. If you have already lost the precision before you get to the dot-product, then there is nothing you can do to get it back. By doing the calculation in double-precision, you introduce no further error, at least for the case where you only have two terms. – Oliver Charlesworth Jan 15 '12 at 13:37
  • 1
    I suppose the problem is that the original float values are only rounded approximations of real numbers to start with, and so although performing the dot product in double precision doesn't introduce any more rounding errors (until you convert the value back to float anyway), the cancellation problem would still exist, because the dot product is operating on approximate values to start with (since catastrophic cancellation occurs when the error in these original values becomes dominant by subtracting them) - converting them to double doesn't make the original values more accurate. – Nick Kovac Jan 16 '12 at 02:15
  • 1
    I guess using doubles globally would help, since the original values would then be more accurate, but still it seems to me that as long as the two vectors were close enough to orthogonal then this problem would still occur. So basically converting to doubles globally might make the problem less frequent, but it doesn't seem to avoid it completely? – Nick Kovac Jan 16 '12 at 02:16
  • There is a good paper about floating point arithmetic: http://www.validlab.com/goldberg/paper.pdf If you search there for 'inner product' the suggestion is also to do the sum of multiplications in double precision. – brita_ Feb 02 '18 at 15:18
  • There is also this hands-on series of posts on floating points for more information on common pitfalls and on the topic in general: https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ – brita_ Feb 08 '18 at 18:28