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

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.

-
For certain values of A and B (where B is the average), it might be helpful to compute the result as (A+B)(A-B). Perhaps a numerical analyst can say how bright/stupid this would be in terms of minimizing roundoff error? Gerhard "Email me About System Design" Paseman. 2011.07.14 – Gerhard Paseman Jul 14 '11 at 21:51

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.

-
The book is on Google books, books.google.com/… , although for some odd reason the pages have been scanned upside down. – Federico Poloni Sep 28 '11 at 11:38
It would help to summarize the algorithm at the link. I don't have time to, but here is a direct link to Section 1.9, which starts on page 11: Higman, p. 11 – nealmcb Jul 7 '15 at 14:41

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.

-
I am not an expert in numerical computation, but your simple formula has low precision when x_i’s are large and the variance is small. – Tsuyoshi Ito Jul 14 '11 at 18:45
I have done some problems that required iterative computation of mean/variance, and I have used this simple formula. It can have catastrophic cancellation. I used a multiprecision library (GMP) to store intermediate results at very high precision (128 bits or so) and that worked fine. – David Harris Jul 14 '11 at 19:47
Is the above formula from Youngs & Cramer written correctly? From a quick check i am not getting the sum of squares. – Luna Jul 14 '11 at 20:00
For the cases rightly considered by Tsuyoshi Ito, check out also: daheiser.info/excel/notes/NOTE%20P.pdf – Luna Jul 14 '11 at 20:14
David, as you are at it, what about if we use something like: ( sum (x^2 / n) - (sum (x / n) )^2 ) where the divisions are carried out before the sum (to avoid the difference between possibly similar large numbers) ? Do you still see the same problem with small variances ? – Luna Jul 15 '11 at 0:48