Numerical Stability of Welford's Algorithm

Hello! Welcome back to our deep dive into Welford's algorithm.

In the last two lessons, we completed the full mathematical derivation of the algorithm. We established the recursive update rules for both the mean ( $μ_{k}$ ) and the sum of squared differences from the mean ( $M_{k}$ ):

$μ_{k} = μ_{k - 1} + \frac{x _{k} - μ _{k - 1}}{k}$
$M_{k} = M_{k - 1} + (x_{k} - μ_{k - 1}) (x_{k} - μ_{k})$

Now we arrive at the crucial "why" behind this course. The purpose of this entire derivation is to create an algorithm that is not only efficient but also numerically stable. Today's lesson addresses this directly. Our learning outcome is to explain why Welford's algorithm avoids catastrophic cancellation and maintains numerical stability. We will connect the mathematical formulas we just derived to the practical limitations of computer arithmetic that we first explored in Module 1.

1. A Refresher on Catastrophic Cancellation

Let's begin by briefly revisiting the core problem: catastrophic cancellation. This phenomenon occurs in floating-point arithmetic when you subtract two very large, nearly equal numbers. Because computers store numbers with finite precision, the leading, most significant digits of the two numbers cancel each other out, leaving a result dominated by the less-significant, and potentially inaccurate, trailing digits.

This isn't just a theoretical curiosity; it can cause real-world algorithms to fail spectacularly.

Numerical Instability in Function Evaluation — This plot shows the function $f(x) = \frac{1 - \cos(x)}{x^2}$ calculated using standard double-precision floating-point numbers for very small values of $x$. Mathematically, as $x \to 0$, the function should approach 0.5. However, because $\cos(x)$ gets extremely close to 1, the subtraction in the numerator `1 - cos(x)` suffers from catastrophic cancellation, leading to the noisy and incorrect results you see.

The key takeaway is that any algorithm that relies on subtracting two large, nearly equal numbers is at risk of producing highly inaccurate results.

2. Why the Naive Algorithm Fails

In Module 1, you analyzed why the naive single-pass variance formula is unstable. Let's formalize that analysis now. The naive algorithm computes the variance using the formula:
$s^{2} = \frac{1}{n - 1} (\sum_{i = 1}^{n} x_{i}^{2} - \frac{( \sum _{i = 1}^{n} x _{i} ) ^{2}}{n})$

Consider a dataset where the standard deviation is much smaller than the mean (e.g., measuring the diameter of a precision-engineered part in millimeters, where values might be 10.001, 10.002, 9.999, ...).

The term $\sum x_{i}^{2}$ will be a very large number.
The term $\frac{( \sum x _{i} ) ^{2}}{n}$ , which is equivalent to $n μ_{n}^{2}$ , will also be a very large number.
If the variance is small, these two large numbers will be almost identical.

The subtraction $\sum x_{i}^{2} - n μ_{n}^{2}$ is a textbook case of catastrophic cancellation waiting to happen. The small, true difference between these values gets lost in the floating-point rounding errors of the large values themselves. This can lead to a final variance that is wildly inaccurate, or even negative—which is mathematically impossible.

To reinforce this, it's helpful to see a practical demonstration.

Ten Little Algorithms, Part 3: Welford's Method (and Friends)

The article 'Ten Little Algorithms, Part 3: Welford's Method (and Friends)' by Jason Sachs, which you've seen before, provides a clear demonstration and explanation of this failure.

Please quickly review the sections showing the Python code for the naive_stats function and its output. Focus on the explanation that starts with 'The problem here is that when we do this calculation: (S2 - (S*S)/n)... we are subtracting two very large numbers'. This directly connects the formula to the problem of precision loss.

3. How Welford's Algorithm Succeeds

Now, let's analyze the Welford update rules to understand why they are immune to this problem.

$M_{k} = M_{k - 1} + (x_{k} - μ_{k - 1}) (x_{k} - μ_{k})$

The stability of this formula comes from what it calculates.
Instead of tracking two large, unrelated sums ( $\sum x_{i}$ and $\sum x_{i}^{2}$ ), Welford's algorithm directly tracks the sum of squared differences from the mean, $M_{k}$ .

Let's break down the update step:

The quantity being updated, $M_{k}$ : This value is directly proportional to the variance. If the variance is small, $M_{k}$ will be a relatively small number throughout the calculation.
The correction term, $(x_{k} - μ_{k - 1}) (x_{k} - μ_{k})$ : This term is the product of two deviations from the mean. If the data has a small standard deviation, these deviations will be small numbers.

The algorithm works by taking the existing sum of squared differences ( $M_{k - 1}$ ) and adding a small correction term to it. Crucially, at no point does it subtract two large, nearly equal numbers. The calculation is always centered around the mean, ensuring that the intermediate values being handled relate to the data's dispersion, not its absolute magnitude.

This insight is the core of the lesson. The following resource explains it concisely.

Welford algorithm for updating variance

The article 'Welford algorithm for updating variance' by Chang-Yao Chen has an excellent section on numerical stability that summarizes this point perfectly.

Please read the section titled 'Numerical stability'. It directly explains how catastrophic cancellation affects the naive variance calculation and states why Welford's algorithm avoids it: 'the two terms on either side of the subtraction are of the same order by design'. (Note: While the Welford update rule is an addition, the principle of working with numbers of similar, appropriate magnitude holds).

Test your understanding!

Let's illustrate the difference in magnitudes with a simple, extreme example. Consider the dataset [1,000,001, 1,000,003]. The mean is $1, 000, 002$ and the sample variance is 2.

For the naive formula:
The two terms you would subtract are $\sum x_{i}^{2}$ and $n μ_{n}^{2}$ .

$\sum x_{i}^{2} = (1 0^{6} + 1)^{2} + (1 0^{6} + 3)^{2} \approx 2 \cdot 1 0^{12}$
$n μ_{n}^{2} = 2 \cdot (1 0^{6} + 2)^{2} \approx 2 \cdot 1 0^{12}$
You are subtracting two numbers on the order of $1 0^{12}$ to find a result related to the variance of 2.

For Welford's algorithm:
Calculate the final sum of squared differences, $M_{2}$ , step-by-step. What are the magnitudes of the numbers involved in the second update step (calculating $M_{2}$ )?

Show answer

Let's trace the Welford calculation.
Initialize: $μ_{0} = 0, M_{0} = 0$ .

Step 1 (for $x_{1} = 1, 000, 001$ ):

$μ_{1} = 1, 000, 001$
$M_{1} = M_{0} + (x_{1} - μ_{0}) (x_{1} - μ_{1}) = 0 + (1000001) (0) = 0$

Step 2 (for $x_{2} = 1, 000, 003$ ):

$μ_{2} = μ_{1} + \frac{x _{2} - μ _{1}}{2} = 1000001 + \frac{1000003 - 1000001}{2} = 1000001 + \frac{2}{2} = 1, 000, 002$
$M_{2} = M_{1} + (x_{2} - μ_{1}) (x_{2} - μ_{2})$
- $x_{2} - μ_{1} = 1000003 - 1000001 = 2$
- $x_{2} - μ_{2} = 1000003 - 1000002 = 1$
- $M_{2} = 0 + (2) (1) = 2$

The final sum of squared differences is 2. To compute this, the crucial update step involved multiplying the small numbers 2 and 1. The algorithm completely avoids calculations with numbers in the $1 0^{12}$ range, thereby preserving numerical precision. This is the power of Welford's algorithm.

Conclusion

In this lesson, we connected the abstract mathematics of Welford's algorithm to the concrete problem of floating-point instability. By analyzing the structure of the update formulas, we've seen precisely how it sidesteps the fatal flaw of the naive single-pass approach.

Key Takeaways:

The naive variance formula is unstable because it involves the subtraction of two very large, nearly equal numbers ( $\sum x_{i}^{2}$ and $n μ_{n}^{2}$ ), a classic recipe for catastrophic cancellation.
Welford's algorithm is stable because it computes the sum of squared deviations from the mean incrementally.
The update step involves adding a small correction term—derived from deviations from the mean—to the running sum. This ensures all calculations are performed on numbers of a manageable magnitude, preserving precision.

Preview of the Next Lesson:

We have now covered the full theory: the problem (numerical instability), the mathematical derivation of the solution (Welford's recursive formulas), and the explanation for why the solution works. In the final module, we will move from theory to practice. You will implement Welford's algorithm, the naive algorithm, and the stable two-pass algorithm in code. We will then design a dataset specifically to trigger catastrophic cancellation and empirically compare the accuracy of the three methods, bringing all the concepts from this course together.

Recommended resources

1. A Refresher on Catastrophic Cancellation

2. Why the Naive Algorithm Fails

3. How Welford's Algorithm Succeeds

Conclusion