Why use Welford's Method over a more naive approach?

I came across this pair of blog posts explaining the problem with calculating variance on floating point data:

https://jvns.ca/blog/2023/01/13/examples-of-floating-point-problems/#example-3-a-variance-calculation-gone-wrong

https://jonisalonen.com/2013/deriving-welfords-method-for-computing-variance/

They show that a totally naive implementation of the variance calculation runs into the problem that the register in which you're accumulating the sum of squares will get so large that the floating point error becomes a fatal problem. In some cases this can even cause the variance to be negative at the end of the calculation.

To solve this, they propose using Welford's method, which (if I understand correctly) keeps a running estimate of the mean as you loop through the numbers, giving you an approximation of $$x_i - \mu$$. The blog author gives the following pseudocode:

variance(samples):
M := 0
S := 0
for k from 1 to N:
x := samples[k]
oldM := M
M := M + (x-M)/k
S := S + (x-M)*(x-oldM)
return S/(N-1)


But then I thought "why can't we just calculate the mean first, and use that when calculating the variance?

variance(samples):
M := 0
S := 0
for k from 1 to N:
M := M + samples[k]
M := M / N
for k from 1 to N:
S := S + (samples[k]-M)**2
return S/(N)


It looks like this should be more accurate and run in about the same amount of time. To test it, I wrote it up in Python and compared numpy to my implementation of Welford, my fairly naive implementation ("good"), and the blogger's example of an overly naive implementation ("bad").

Numpy variance:     0.00029272912371828137
Good variance:      0.0002927291237197718
Bad variance:       -1.0752
Welford variance:   0.0002927584024568408

Numpy time: 2.579808235168457e-05
Good time:  0.0026930201053619383
Bad time:   0.0018948330879211427
Welford time:   0.004232816934585572

Good method suffered 42.12% peformance hit
Welford method suffered 123.39% performance hit


My expectations seem to have been borne out, assuming that the numpy implementation does it "right". So, is Welford actually the best choice for this for some other reason? Did I make a mistake somewhere in my implementation? Did the blog author just think it was cool so he decided to write a blog post? All of the code I wrote is given below:

import numpy, time

def calculate_bad_variance(nums):
sum_of_squares = 0
sum_of_nums = 0
N = len(nums)
for num in nums:
sum_of_squares += num**2
sum_of_nums += num
mean = sum_of_nums / N
variance = (sum_of_squares - N * mean**2) / N

return variance

def calculate_good_variance(nums):
sum_of_squares = 0
sum_of_nums = 0
N = len(nums)
for num in nums:
sum_of_nums += num
mean = sum_of_nums / N
for num in nums:
sum_of_squares += (num - mean)**2
variance = (sum_of_squares) / N

return variance

def welford_variance(nums):
M = 0
S = 0
N = len(nums)
for k in range(N):
x = nums[k]
oldM = M
M = M + (x-M)/(k+1)
S = S + (x-M)*(x-oldM)
return S/(N-1)

nums = np.random.uniform(10000000., 10000000.06, 10000)
N = 1000

true_var = np.var(nums)
bad_var = calculate_bad_variance(nums)
good_var = calculate_good_variance(nums)
welford_var = welford_variance(nums)

start = time.time()
for _ in range(N):
np.var(nums)
np_t = (time.time() - start)/N

start = time.time()
for _ in range(N):
calculate_bad_variance(nums)
bad_t = (time.time() - start)/N

start = time.time()
for _ in range(N):
calculate_good_variance(nums)
good_t = (time.time() - start)/N

start = time.time()
for _ in range(N):
welford_variance(nums)
welford_t = (time.time() - start)/N

print(f"Numpy variance:\t\t{true_var}\nGood variance:\t\t{good_var}\nBad variance:\t\t{bad_var}\nWelford variance: \t{welford_var}\n")
print(f"Numpy time:\t{np_t}\nGood time:\t{good_t}\nBad time:\t{bad_t}\nWelford time:\t{welford_t}\n")
print(f"Good method suffered {((good_t/bad_t)-1):.2%} peformance hit")
print(f"Welford method suffered {((welford_t/bad_t)-1):.2%} performance hit")
$$$$


2 Answers

The fact that you first pass every value to calculate mean and then calculate variance means you need to store data and that's the problem ! For a lot of applications you want to be able to compute "on the fly"

• This actually uses less state than Welford, since I got rid of the oldM register. Welford still takes advantage of the M` register to calculate a running mean approximation. I agree that it's "interruptible", i.e. I can stop execution of Welford at any time, and the registers will have an approximate variance, whereas with mine you need to at least finish calculating the mean before-hand. I'm not really sure how that's useful, though.
– fpf3
Oct 19, 2023 at 19:12

It’s a good algorithm if you want a single pass through your data. If it is not precise enough you have two choices:

1. A more complicated algorithm. Calculate a rough estimate of the mean. Then get the precise mean: Add the differences between data points and the estimated mean, reducing your numbers. And add numbers > mean and numbers < mean separately: If the sum is > 0 add the next number < mean first, otherwise the next number > mean. So the sum stays around zero with small rounding errors.

Then you calculate the average of (data point - mean) squared. You use exactly the same method. So all in all four passes.

1. Use higher precision. If your computer has well-implemented floating-point arithmetic then you can calculate the exact rounding error of a sum a+b: calculate the sum s. If |a| >= |b| then the exact error is (a - s) + b, otherwise it is (b - s) + a. So you have one sum where you add up values, and one where you add up rounding errors, and finally you add both sums.