Floating point addition: loss-of-precision issues

Question

In short: how can I execute a+b such that any loss-of-precision due to truncation is away from zero rather than toward zero?

The Long Story

I'm computing the sum of a long series of floating point values for the purpose of computing the sample mean and variance of the set. Since Var(X) = E(X²) - E(X)², it suffices to maintain running count of all numbers, the sum of all numbers so far, and the sum of the squares of all numbers so far.

So far so good.

However, it's absolutely required that E(X²) > E(X)², which due to floating point accuracy isn't always the case. In pseudo-code, the problem is this:

int count;
double sum, sumOfSquares;
...
double value = <current-value>;
double sqrVal = value*value; 

count++;
sum += value; //slightly rounded down since value is truncated to fit into sum
sumOfSquares += sqrVal; //rounded down MORE since the order-of-magnitude 
//difference between sqrVal and sumOfSquares is twice that between value and sum;

For variable sequences, this isn't a big issue - you end up slightly under-estimating the variance, but it's often not a big issue. However, for constant or almost-constant sets with a non-zero mean, it can mean that E(X²) < E(X)², resulting in a negative computed variance, which violates expectations of consuming code.

Now, I know about Kahan Summation, which isn't an attractive solution. Firstly, it makes the code susceptible to optimization vagaries (depending on optimization flags, code may or may not exhibit this problem), and secondly, the problem isn't really due to the precision - which is good enough - it's because addition introduces systematic error towards zero. If I could execute the line

sumOfSquares += sqrVal;

in such a way as to ensure that sqrVal is rounded up, not down, into the precision of sumOfSquares, I'd have a numerically reasonable solution. But how can I achieve that?

_{Edit: Finished question - why does pressing enter in the drop-down-list in the tag field submit the question anyhow?}

Answer 1

IEEE provides four rounding modes, (toward -inf, toward +inf, toward 0, tonearest). Toward +inf is what you seem to want. There is no standard control in C90 or C++. C99 added the header <fenv.h> which is also present as an extension in some C90 and C++ implementation. To respect the C99 standard, you'd have to write something like:

#include <fenv.h>
#pragma STDC FENV_ACCESS ON

int old_round_mode = fegetround();
int set_round_ok = fesetround(FE_UPWARD);
assert(set_round_ok == 0);
...
int set_round_ok = fesetround(old_round_mode);
assert(set_round_ok == 0);

It is well known that the algorithm you use is numerically unstable and has precision problem. It is better for precision to do two passes on the data.

Answer 2

There's another single-pass algorithm which rearranges the calculation a bit. In pseudocode:

n = 0
mean = 0
M2 = 0

for x in data:
    n = n + 1
    delta = x - mean
    mean = mean + delta/n
    M2 = M2 + delta*(x - mean)  # This expression uses the new value of mean

variance_n = M2/n         # Sample variance
variance = M2/(n - 1)     # Unbiased estimate of population variance

(Source: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance )

This seems better behaved with respect to the issues you pointed out with the usual algorithm.

Answer 3

If you don't worry about the precision, but just about a negative variance, why don't you simply do V(x) = Max(0, E(X^2) - E(X)^2)