Numerically most robust way to compute sum of products (standard deviation) in floating-point?

Question

I stumbled across a paper by Welford (1962), where he proclaims a method that should compute the standard deviation numerically more robust than the naive algorithms (http://www.jstor.org/stable/1266577). Here, "numerically robust" means that round-off errors are reduced.

He gives a recurrence for the sum of squares $S_n = \sum_{i=1}^n (x_i - \mu_n)^2 = S_{n-1} + \frac{n-1}{n} (x_n - \mu_{n-1})^2 $ , $\mu$ being, of course, the mean. That way, the standard deviation can be computed iteratively in a single pass.

As far as I understand, he claims that his iterative recurrence formula is numerically more robust, because all terms in it are of the same order (provided the input data are all of the same order).

This is what I don't understand. It seems to me that, as $n$ gets incremented, $S_n$ becomes larger and larger, so more and more significant digits from the second term are lost, aren't they?

Googling a bit further, I have found a paper by Youngs & Cramer, 1971, who looked at a number of methods, including Welford's, of computing the sum of products / standard deviation more robustly (http://www.jstor.org/stable/1267176).

Conducting a number of experiments, they found that Welford's method does not provide any benefits. So that seems to confirm my doubts about Welford's method.

Now, Youngs & Cramer propose another method, which computes $S_n = S_{n-1} + \frac{n-1}{n}(nx_n - s_n)^2 $, where $s_n = s_{n-1} + x_n$.
Empirically, they found their method to be superior.

Again, I don't understand why this should be the case: isn't there some catastrophic cancellation going on in $(nx_n - s_n)$ ? Don't the terms $S_n$ and the fraction differ in their magnitude more and more, so that more and more digits of the fraction term get rounded off?

I would by most grateful if somebody could shed some light on these questions.
In addition, I'd like to know which is the best method (in terms of roundoff errors) to compute these sums. Surprisingly, I haven't found anything about this in Numerical Recipes (or I overlooked it).

Thank you very much in advance. Gabriel.

Answer 1

I found a discussion of this exact problem in Higham, Accuracy and stability of numerical algorithms, Section 1.9.

The author suggests an alternative algorithm and claims that it is numerically stable (in the mixed backward-forward sense); proofs for the precise accuracy bounds of the formulas are left as exercises.

Answer 2

For the variance (square of std), personally I would use this simple one:

( n (sum x*x ) - (sum x)*(sum x) ) / (n*n) or: ( (sum x*x ) - (sum x)*(sum x) / n ) / n

[ For the case of relatively small variance ("catastrophic cancellation") check also: http://www.daheiser.info/excel/notes/NOTE%20P.pdf ]

Edit:

A viki page also reports a nice discussion: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#cite_note-1

In particular, it mentions: A numerically stable algorithm is given below. It also computes the mean. This algorithm is due to Knuth,[1] who cites Welford.[2] [...] This algorithm is much less prone to loss of precision due to massive cancellation, but might not be as efficient because of the division operation inside the loop.

About the 2-pass algo, it says: "This algorithm is often more numerically reliable than the naïve algorithm for large sets of data, although it can be worse if much of the data is very close to but not precisely equal to the mean and some are quite far away from it."

I'd say that if the variance is not small (as often happens) and ease and speed of computation are a priority, the simpler version (closer to the definition) could still be considered.