Case Study 2: Building a Variance You Can Trust

DataField.Dev

Case Study 2: Building a Variance You Can Trust

"The textbook formula for variance is one subtraction away from a disaster. The fix is not a faster computer — it is a better arrangement of the same arithmetic."

Executive Summary

Where Case Study 1 diagnosed an order-dependent sum, this one builds a numerically trustworthy replacement for a routine every scientific code needs: computing a mean and a variance. The famous one-pass "computational formula" for variance — the mean of the squares minus the square of the mean — is seductive because it needs a single loop, and it is a textbook example of catastrophic cancellation (§20.3): for data with a large mean it can return a badly wrong answer, and even a negative variance, which is mathematically impossible. We will make the failure concrete, build the stable two-pass algorithm, validate it against a quadruple-precision reference oracle, and add compensated summation so the routine is safe to drop into the heat solver's diagnostics.

Skills applied - Recognizing catastrophic cancellation in the one-pass variance formula (§20.3). - Designing a stable algorithm that subtracts the mean before squaring (§20.5, stability). - Using quadruple precision as a reference oracle, and knowing its cost (§20.6). - Choosing dp deliberately and reasoning about accumulated error (§20.2, §20.6).

Background

Given data $x_1, \dots, x_n$ with mean $\bar{x}$, the (population) variance is $\sigma^2 = \frac{1}{n}\sum_i (x_i - \bar{x})^2$. Algebra lets you expand that into the one-pass formula

$$ \sigma^2 = \frac{1}{n}\sum_i x_i^2 \;-\; \bar{x}^2, $$

which is appealing because you can accumulate $\sum x_i$ and $\sum x_i^2$ in a single sweep and never store the data. Mathematically the two forms are identical. Numerically they are worlds apart, and the difference is the entire subject of this chapter.

The danger is in that final subtraction. When the data has a large mean, both $\frac{1}{n}\sum x_i^2$ and $\bar{x}^2$ are large and nearly equal — and the true variance is the small difference between them. Subtracting two nearly equal large numbers is catastrophic cancellation: the leading digits, which agree, cancel away, and what remains is dominated by the rounding errors in the two big quantities. You can lose every significant digit of the answer.

Phase 1 — Watch the One-Pass Formula Fail

Take the smallest dataset that exhibits the problem: two values near $10^{8}$ that differ by $1$.

The true mean is $10^{8} + 0.5$, and the true population variance is $\sigma^2 = \frac{(-0.5)^2 + (0.5)^2}{2} = 0.25$. Now trace the one-pass formula in double precision, using the ULP facts from §20.2 (near $10^{16}$, one ULP is $2$):

$x_2^2 = (10^{8}+1)^2 = 10000000200000001$ — an odd integer just above $10^{16}$, so it rounds to the nearest even representable value, $10000000200000000$. One unit already gone.
$\frac{1}{2}\sum x_i^2 = \frac{10^{16} + 10000000200000000}{2} = 10000000100000000$.
$\bar{x}^2 = (10^{8}+0.5)^2 = 10000000100000000.25$, which rounds to $10000000100000000$.
$\sigma^2_{\text{naive}} = 10000000100000000 - 10000000100000000 = 0$.

The formula reports a variance of zero for two distinct numbers — a 100% error, produced entirely by cancellation. The correct $0.25$ was hiding in digits that both large quantities had already rounded away.

Phase 2 — Build the Stable Two-Pass Algorithm

The cure is to never form that subtraction. Compute the mean first, then subtract it from each datum before squaring, so the squared quantities are small and no cancellation occurs. This is the two-pass algorithm, and it is the stable one:

module kinds
  implicit none
  private
  public :: dp, qp
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: qp = selected_real_kind(33, 4931)   ! quad; platform-dependent
end module kinds

module stats
  use kinds
  implicit none
contains
  pure function var_naive(x) result(v)      ! one-pass: mean of squares - square of mean
    real(dp), intent(in) :: x(:)
    real(dp)             :: v, mean
    integer              :: n
    n    = size(x)
    mean = sum(x) / real(n, dp)
    v    = sum(x**2) / real(n, dp) - mean**2
  end function var_naive

  pure function var_twopass(x) result(v)    ! two-pass: subtract the mean, THEN square
    real(dp), intent(in) :: x(:)
    real(dp)             :: v, mean
    integer              :: n
    n    = size(x)
    mean = sum(x) / real(n, dp)
    v    = sum((x - mean)**2) / real(n, dp)
  end function var_twopass

  pure function var_quad(x) result(v)       ! reference oracle: one-pass, but in quad precision
    real(dp), intent(in) :: x(:)
    real(dp)             :: v
    real(qp)             :: mean
    integer              :: n
    n    = size(x)
    mean = sum(real(x, qp)) / real(n, qp)
    v    = real(sum(real(x, qp)**2) / real(n, qp) - mean**2, dp)
  end function var_quad
end module stats

program variance_trust
  use kinds, only: dp
  use stats, only: var_naive, var_twopass, var_quad
  implicit none
  real(dp) :: hard(2) = [ 1.0e8_dp, 1.0e8_dp + 1.0_dp ]
  real(dp) :: text(8) = [ 2.0_dp, 4.0_dp, 4.0_dp, 4.0_dp, 5.0_dp, 5.0_dp, 7.0_dp, 9.0_dp ]

  print '(a)',       '--- {1e8, 1e8+1}: true population variance = 0.25 ---'
  print '(a, f8.4)', 'naive one-pass  = ', var_naive(hard)
  print '(a, f8.4)', 'two-pass        = ', var_twopass(hard)
  print '(a, f8.4)', 'quad reference  = ', var_quad(hard)
  print '(a)',       '--- {2,4,4,4,5,5,7,9}: true population variance = 4.0 ---'
  print '(a, f8.4)', 'two-pass        = ', var_twopass(text)
end program variance_trust

$ gfortran -std=f2018 -Wall case2.f90 -o case2 && ./case2
--- {1e8, 1e8+1}: true population variance = 0.25 ---
naive one-pass  =   0.0000
two-pass        =   0.2500
quad reference  =   0.2500
--- {2,4,4,4,5,5,7,9}: true population variance = 4.0 ---
two-pass        =   4.0000

The two-pass algorithm recovers the correct 0.2500 where the one-pass gave 0.0000. In the two-pass form the deviations are exactly $-0.5$ and $+0.5$, their squares exactly $0.25$, and nothing cancels. On the well-behaved textbook dataset (whose deviations from the mean of $5$ are the small integers $-3,-1,-1,-1,0,0,2,4$, with squares summing to $32$), it returns the exact 4.0000.

Phase 3 — Validate Against a Quad-Precision Oracle

How do we know 0.25 is right and 0.0 is wrong, rather than the reverse? By computing the answer a third way, in a precision so much higher that its own rounding is negligible, and treating that as ground truth. That is what var_quad does: it runs the same one-pass formula — the unstable one — but in quadruple precision (§20.6), whose 113-bit significand carries about 34 decimal digits. At that width, $(10^{8}+1)^2$ and $(10^{8}+0.5)^2$ are represented exactly, the cancellation removes only correct digits, and the result is the true 0.25.

This is a technique worth internalizing: use a higher precision as an oracle to test a lower one. When a double-precision result is suspect, recompute it in quad; if they agree, cancellation was not a problem; if they diverge, you have found your instability. The quad oracle confirms our hand analysis: the two-pass double result (0.25) matches the oracle, and the naive double result (0.0) does not.

⚡ Performance Note — the oracle is expensive. Quadruple precision on mainstream CPUs is software emulated (via libquadmath in gfortran), commonly tens of times slower than hardware double. That is fine for a one-off validation of a small dataset, and unacceptable as the production path over a billion-cell grid. Use quad to check, then ship the stable double algorithm. Note too that selected_real_kind(33, 4931) returns a negative value on a platform with no quad type at all, which would make the code fail to compile — a portable code should guard against that or fall back to double.

Phase 4 — Going Further: Compensated Summation

The two-pass algorithm removes the cancellation, but both passes still contain a large summation, and Case Study 1 showed that long summations drift through absorption. For very large $n$, the remaining error lives in those sums, and the tool for it is compensated (Kahan) summation — the fix Case Study 1 pointed toward. The idea is to keep a small running correction c that captures the low-order bits lost at each addition, and feed them back in:

pure function kahan_sum(v) result(s)
  use kinds, only: dp
  real(dp), intent(in) :: v(:)
  real(dp)             :: s, c, y, t
  integer              :: i
  s = 0.0_dp
  c = 0.0_dp                  ! the running compensation
  do i = 1, size(v)
     y = v(i) - c             ! bring back what was lost last time
     t = s + y                ! this addition may lose the low bits of y ...
     c = (t - s) - y          ! ... and this recovers exactly what was lost
     s = t
  end do
end function kahan_sum

Substituting kahan_sum for the plain sum in the two passes gives an algorithm whose error barely grows with $n$ at all — the running sum behaves almost as if it had twice the precision. For the heat solver's energy diagnostic over a large grid, this is the difference between a conservation check that holds to a few ULPs and one that drifts visibly over a long run. The cost is three extra floating-point operations per element and one extra variable — a small, constant price for order-independent, near-exact sums.

Phase 5 — Package It for the Solver, and Sanity-Check

The solver gains a small diagnostics module exposing field_variance(u) built on var_twopass and, for long runs, kahan_sum. Three checks confirm it before it goes in:

Zero variance on a flat field. A uniform plate must report exactly 0.0. The two-pass form guarantees it (every deviation is exactly zero); the one-pass form does not.
Shift invariance. Adding a constant to every cell must not change the variance. The two-pass form is (nearly) shift-invariant because it subtracts the mean; the one-pass form fails this test badly, as Phase 1 showed with the $10^{8}$ offset.
Oracle agreement. On a representative field, the double two-pass result matches the quad oracle to full double precision.

With those three green, the diagnostic is trustworthy — and, crucially, its correctness rests on the algorithm's structure (subtract the mean first, compensate the sums), not on hoping double precision is "good enough." That is the §20.5 lesson made operational: the problem (computing a variance) is well-conditioned; the one-pass algorithm was unstable; we replaced it with a stable one.

Discussion Questions

The one-pass formula can return a negative variance. Construct a small dataset (by hand or by experiment) for which it does, and explain why the two-pass formula cannot.
var_quad runs the unstable one-pass formula yet gives the right answer. Does higher precision fix the instability, or merely hide it? What happens to the quad version as the mean grows toward $10^{16}$?
The two-pass algorithm reads the data twice, which is impossible for a true streaming computation (data arriving once, too large to store). Look up Welford's online algorithm; what does it trade to get stability in a single pass?
When would the one-pass formula actually be safe to use? Characterize the datasets for which it and the two-pass form agree to full precision.

Your Turn: Extensions

Option A (build). Implement var_kahan by substituting kahan_sum for sum in both passes of var_twopass, and design a dataset of $10^{6}$ values where it visibly beats the plain two-pass form.
Option B (design). Implement Welford's single-pass online variance and compare it with the two-pass algorithm on the $\{10^{8}, 10^{8}+1\}$ case and on a large random dataset. Which is more accurate? Which is faster?
Option C (optimize). Add field_variance to the heat solver as a per-diagnostic-step call, and use system_clock (Chapter 28) to measure the overhead of two-pass versus one-pass versus Kahan on your grid. Is the stable version's cost ever noticeable next to the solver's stencil update?

Key Takeaways

The one-pass variance formula ($\overline{x^2} - \bar{x}^2$) is a classic catastrophic cancellation trap: for data with a large mean it loses all accuracy and can return a negative variance.
The two-pass algorithm — subtract the mean, then square — is stable: it never subtracts nearly equal large numbers, and it returns the correct answer for the same well-conditioned problem.
Use quadruple precision as an oracle to test a double-precision result; it is authoritative but slow and platform-dependent, so validate with it and ship double.
For very large sums, add compensated (Kahan) summation to control the residual drift. Stability is a property you build into the algorithm, not one you buy with more precision.