12

Can anyone recommend any C++ libraries/routines/packages that contain strategies for maintaining the stability of various floating point operations?

Example: suppose you would like to sum across a vector/array of one million long double in the unit interval (0,1), and that each number is of about the same order of magnitude. Naively summing for (int i=0;i<1000000;++i) sum += array[i]; is unreliable - for large enough i, sum will be of a much larger order of magnitude than array[i], and so sum += array[i] would be equivalent to sum += 0.00. (Note: the solution to this example is a binary summing strategy.)

I deal with sums and products of thousands/millions of miniscule probabilities. I am using the arbitrary-precision library MPFRC++ with a 2048 bit significand, but the same concerns still apply.

I am chiefly concerned with:

  1. Strategies for accurately summing many numbers (e.g. above Example).
  2. When is multiplication and division potentially unstable? (If I want to normalize a large array of numbers, what should my normalization constant be? The smallest value? The largest? A median?)
cmo
  • 3,762
  • 4
  • 36
  • 64

1 Answers1

6

Binary summation doesn't guarantee accurate result. The most reliable (albeit slower) method is to use Kahan summation. Boost.Accumulators has an implementation of the above and much more.

Multiplication and division stability: unless you get to denormalized floats they don't suffer from the same problems as summation and substraction. In fact the error of multiplication is at most 0.5 ulp (units last place).

... what should my normalization constant be?

What do you mean by 'normalize'? It depends on the norm you use. Possible candidates: use the maximum absolute value in the array, or any other generalized mean. (Other choices you listed do not work since they may be zero even for non-zero array.)

Yakov Galka
  • 70,775
  • 16
  • 139
  • 220
  • By "normalize", I mean: At step X, I have a set of miniscule-but-non-zero probabilities. At step X+1, I will multiply these numbers by even more miniscule probabilities. According to the details of my mathematical model, I am allowed to multiply the numbers at step X by any constant (i.e. normalize) before moving to step X+1. This has the practical benefit of helping avoid underflow. – cmo Jun 29 '12 at 16:26
  • @CycoMatto: then you should multiply by a power of two once in a while. The only way to avoid underflow/overflow is to guarantee that the ratio between the largest and the smallest numbers is representable with the floating points at hand (if not, you can't do anything). – Yakov Galka Jun 29 '12 at 16:36
  • 1
    Alternatively store their logarithms and use summation. Can't say about the stability of this approach though. – Yakov Galka Jun 29 '12 at 16:37
  • Why a power of 2? Currently, I am multiplying by the inverse of the smallest number at step X - this way, the numbers are all >= 1. Multiplying by probabilities at step X+1 will then "drag the numbers back toward 0". Right? – cmo Jun 29 '12 at 16:47
  • I was also thinking of multiplying by some nice power of 10. But I suppose a power of 2 makes more sense given that numbers are stored in binary. – cmo Jun 29 '12 at 16:48
  • You use a kind of "inverse transformation". This rule may help: "The specific transformation used depends on the extent of the deviation from normality. If the distribution differs moderately from normality, a square root transformation is often the best. A log transformation is usually best if the data are more substantially non-normal. An inverse transformation should be tried for severely non-normal data." http://dss.princeton.edu/online_help/analysis/regression_intro.htm – SChepurin Jun 29 '12 at 17:22
  • @SChepurin thank you, but I can't use any transformation - only constants. (constants will "cancel" in any division and can traverse sums, and this will not affect the mathematics). – cmo Jun 29 '12 at 18:05
  • 1
    @CycoMatto: multiplying or dividing by a power of two does not degrade accuracy *at all*. I just want to note that if you have the concerns of underflow, and you normalize it as you described, then you also have to think about possible overflows during normalization. – Yakov Galka Jun 29 '12 at 18:25