Understand the need for Welford's online algorithm

Question

I am puzzled by the Wikipedia entry discussing many online algorithms for computing the sample variance, including the Welford's online algorithm.

In particular, the sample variance $s_n^2$ can be computed from the sample $\{X_1,..,X_n\}$ as: $$ x_n = \frac{n-1}n x_{n-1} + \frac 1n X_n $$ $$ s_n^2=\frac{n-2}{n-1}s_{n-1}^2+\frac 1n(X_n-x_{n-1})^2 $$ with the caveat that

These formulas suffer from numerical instability, as they repeatedly subtract a small number from a big number which scales with n.

I can see it, because for $n\to\infty$ the first term is of the order of the variance, while the second goes to 0. I can see the numerical issue arising from the finite precision of a float.

The entry goes on arguing that the Welford's online algorithm is designed to tackle this issue: $$ M^2_n = M^2_{n-1}+(X_n-x_n)(X_n-x_{n-1}),\qquad s^2_n=\frac{M^2_n}{n-1}. $$

My question is: how does it solve it? I don't see it happening. $M_{n-1}^2$ is a sum of terms looking like the second one, so not only is prone to overflow, but also to become too big with respect to the second term, that may eventually become smaller than the last digit of $M_{n-1}^2$. It seems to me that no improvement is achieved, apart of risking floating point overflow ..

Answer 1

Wikipedia is wrong. The first of two formulas you quote are not in fact numerically unstable; they just add a small number to a large number (which is what Welford's formula also does), and that's a reasonable thing to do as long as you have enough digits in your floating point number to not lose too much accuracy. With today's double precision floating point accuracy, you need to be in the range of $n=10^{12}$ samples so that you have the correct answer to less than 4 digits. For most cases, that's good enough.

The correct comparison the Wikipedia page should be making is with the first set of formulas on the page: $$ \sigma^2 = \overline{(x^2)} - \bar x^2 = \displaystyle\frac {\sum_{i=1}^N x_i^2 - (\sum_{i=1}^N x_i)^2/N}{N}. $$ This formula is unstable because it subtracts two numbers of (generally) equal size, and can therefore lead to catastrophic cancellation of digits.

As a side remark, I agree with the statement that not dividing every time by something of size $n$ is bad because it will lead to overflow. The first set of formulas in the Wikipedia page's section on Welford's algorithm do this right.

---

Wikipedia, like this page, is a volunteer project. You might want to take the opportunity and update the text there!