Running means and covariances

Welfords algorithmPermalink

The mean of a string of numbers is

μ_{n} = (x_{1} + x_{2} + \dots x_{n}) / n

and this gives a simple update equation for a running mean, when updating online.

$μ_{0} = 0$ and $μ_{n} = μ_{n - 1} \frac{n - 1}{n} + x_{n} \frac{1}{n}$

The similar formlua for variance is by accumulating the sum of squares $S_{n} = \sum_{i = 1}^{n} x_{i} x_{i}^{T} = S_{n - 1} + x_{n} x_{n}^{T}$ , $S_{0} = 0$ and then compute

Σ_{n} = \frac{1}{n} \sum_{i = 1}^{n} (x_{i} - μ_{n}) (x_{i} - μ_{n})^{T} = \frac{1}{n} (\sum_{i = 1}^{n} x_{i} x_{i}^{T}) - μ_{n}^{2} = \frac{S_{n}}{n} - μ_{n}^{2}

The algorithm is not great, since we compute the answer as a difference between two large quantities, and may lose numeric precision.

We want to find an accumulator that is “closer” to the answer, not requireing a subtraction at the end. Examine the difference of sum of central squares, $M_{n} = n Σ_{n}$ . We obtain

\begin{aligned} Σ_{n} & = M_{n} / n \\ M_{0} & = 0 \\ M_{n} & = M_{n - 1} + \frac{n - 1}{n} (x_{n} - μ_{n - 1}) (x_{n} - μ_{n - 1})^{T} \end{aligned}

The proof of this formula is quite simple. Some steps use the recursion $μ_{n} = μ_{n - 1} + \frac{x_{n} - μ_{n - 1}}{n}$ . Others use $\sum_{i = 1}^{n - 1} x_{i} = (n - 1) μ_{n - 1}$ . The rest uses straightforward algebra. the computation below is for scalar data, but the generalization to vector valued data is straightforward since I won’t use commutativity. Just replace multiplication with outer products.

\begin{aligned} M_{n} - M_{n - 1} \\ = \sum_{i = 1}^{n} (x_{i} - μ_{n})^{2} - \sum_{i = 1}^{n - 1} (x_{i} - μ_{n - 1})^{2} \\ = (x_{n} - μ_{n})^{2} + \sum_{i = 1}^{n - 1} [(x_{i} - μ_{n})^{2} - (x_{i} - μ_{n - 1})^{2}] \\ = (x_{n} - μ_{n})^{2} + \sum_{i = 1}^{n - 1} [(μ_{n} - x_{i}) μ_{n} + (μ_{n - 1} - μ_{n}) x_{i} + (x_{i} - μ_{n - 1}) μ_{n - 1}] \\ = (x_{n} - μ_{n})^{2} + (n - 1) (μ_{n} - μ_{n - 1})^{2} \\ = \frac{(n - 1)^{2}}{n^{2}} (x_{n} - μ_{n - 1})^{2} + \frac{n - 1}{n^{2}} (x_{n} - μ_{n - 1})^{2} \\ = \frac{n - 1}{n} (x_{n} - μ_{n - 1})^{2} \\ = (x_{n} - μ_{n - 1})^{2} - (x_{n} - μ_{n - 1})^{2} / n \end{aligned}

The paper Comparison of Several Algorithms for Computing Sample Means and Variances by Robert F. Ling (1974) doi attributes the final equation to Welford (1962), but also posits some alternative algorithms and compare them. While writing this article, I found that others have written this or very similar blog post before me, e.g. Joni Salonen, Jakob Ludewig and John D Cook. The formulas reported in their posts may differ slightly from mine, and might be better. Do your own work if you want to be sure which one is the best for your data. You can also read about this kind of algorithms on Wikipedia.

Ludvig Hult

Welfords algorithmPermalink