Welford’s streaming variance

Paper

Note on a Method for Calculating Corrected Sums of Squares and Products, 1962. Technometrics 4(3), pp. 419–420. Read Finch’s later derivation, §§1–3.

The result

Starting from an empty state, the recurrence computes the exact arithmetic mean and the exact sum of squared deviations for any finite list of real observations.

δ = x − μ; μ′ = μ + δ/(n + 1) M₂′ = M₂ + δ(x − μ′) M₂ = Σᵢ(xᵢ − μ)²

Scope & assumptions

8 checked theorems. Update identities, nonnegative scatter increments, count and moment invariants, and full-list mean/scatter correctness. Products and covariances from the original paper are not covered; neither is floating-point stability.

Exact real arithmetic. The mean theorem requires a nonempty list; the scatter theorem also handles the empty list. The executable structure of compute processes the tail first, then adds the head. The update operation itself takes one new observation. M₂ is an unnormalized sum: divide by n for population variance, or n − 1 for sample variance when n > 1.

Prior work

Related Lean work was found in lutar-lean PR #180: an integer running-total invariant labeled as a weighted mean. Its inspected source does not implement the real mean/scatter recurrence or prove the centered sum-of-squares result here. This artifact was written independently. Priority for this stronger scope remains unconfirmed; see the exact source comparison in the search record.

Lean proof

Lean + Mathlib 4.33.1 · 8 checked theorems

import Mathlib.Algebra.BigOperators.Ring.List
import Mathlib.Data.Real.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Positivity

/-! Welford (1962): exact correctness of incremental mean and squared scatter.
The list computation is a fold from the tail; update can process new samples
online. Exact real arithmetic only, with no roundoff guarantee. -/
namespace PaperToLean.Welford
noncomputable section
structure State where
  count : ℕ
  mean : ℝ
  m2 : ℝ

def update (s : State) (x : ℝ) : State :=
  let μ := s.mean + (x - s.mean) / (s.count + 1 : ℝ)
  ⟨s.count + 1, μ, s.m2 + (x - s.mean) * (x - μ)⟩
def compute : List ℝ → State
  | [] => ⟨0, 0, 0⟩
  | x :: xs => update (compute xs) x

theorem mean_accumulator (s : State) (x : ℝ) :
    ((update s x).count : ℝ) * (update s x).mean = (s.count : ℝ) * s.mean + x := by
  have hn : (s.count : ℝ) + 1 ≠ 0 := by positivity
  simp only [update, Nat.cast_add, Nat.cast_one]
  field_simp
  ring

theorem square_accumulator (s : State) (x : ℝ) :
    (update s x).m2 + ((update s x).count : ℝ) * (update s x).mean ^ 2 =
      s.m2 + (s.count : ℝ) * s.mean ^ 2 + x ^ 2 := by
  have hn : (s.count : ℝ) + 1 ≠ 0 := by positivity
  simp only [update, Nat.cast_add, Nat.cast_one]
  field_simp
  ring

theorem scatter_increment (s : State) (x : ℝ) :
    (update s x).m2 = s.m2 + (s.count : ℝ) / (s.count + 1 : ℝ) * (x - s.mean) ^ 2 := by
  have hn : (s.count : ℝ) + 1 ≠ 0 := by positivity
  simp only [update]
  field_simp
  ring

theorem update_nonneg (s : State) (x : ℝ) (hs : 0 ≤ s.m2) : 0 ≤ (update s x).m2 := by
  rw [scatter_increment]
  positivity

/-- Both moment invariants and the count are derived from the computation. -/
theorem moments (xs : List ℝ) :
    (compute xs).count = xs.length ∧
    ((compute xs).count : ℝ) * (compute xs).mean = xs.sum ∧
    (compute xs).m2 + ((compute xs).count : ℝ) * (compute xs).mean ^ 2 =
      (xs.map (· ^ 2)).sum := by
  induction xs with
  | nil => simp [compute]
  | cons x xs ih =>
    refine ⟨by simp [compute, update, ih.1], ?_, ?_⟩
    · rw [compute, mean_accumulator, ih.2.1, List.sum_cons]
      ring
    · rw [compute, square_accumulator, ih.2.2, List.map_cons, List.sum_cons]
      ring

theorem mean_correct (xs : List ℝ) (hn : xs ≠ []) :
    (compute xs).mean = xs.sum / xs.length := by
  have hp : 0 < xs.length := List.length_pos_iff.mpr hn
  have hne : (xs.length : ℝ) ≠ 0 := by exact_mod_cast (ne_of_gt hp)
  apply (eq_div_iff hne).mpr
  have h := (moments xs).2.1
  rw [(moments xs).1] at h
  nlinarith

theorem centered_sum (xs : List ℝ) (μ : ℝ) :
    (xs.map (fun x => (x - μ) ^ 2)).sum =
      (xs.map (· ^ 2)).sum - 2 * μ * xs.sum + (xs.length : ℝ) * μ ^ 2 := by
  induction xs with
  | nil => simp
  | cons x xs ih =>
    simp only [List.map_cons, List.sum_cons, List.length_cons, Nat.cast_add, Nat.cast_one, ih]
    ring

theorem scatter_correct (xs : List ℝ) :
    (compute xs).m2 = (xs.map (fun x => (x - (compute xs).mean) ^ 2)).sum := by
  obtain ⟨hc, hμ, hM⟩ := moments xs
  rw [centered_sum, ← hM, ← hμ, ← hc]
  ring

#print axioms mean_accumulator
#print axioms square_accumulator
#print axioms scatter_increment
#print axioms update_nonneg
#print axioms moments
#print axioms mean_correct
#print axioms centered_sum
#print axioms scatter_correct
end
end PaperToLean.Welford

Lean 4.33.1
Dependencies: Mathlib 4.33.1 (see bundled lake-manifest.json)
File: Welford.lean
SHA-256: a40c05813a0849ee670b1eba00e08ade67868b3cb664011ffda2fc7f32f457eb
Result: passed

'PaperToLean.Welford.mean_accumulator' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.square_accumulator' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.scatter_increment' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.update_nonneg' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.moments' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.mean_correct' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.centered_sum' depends on axioms: [propext, Classical.choice, Quot.sound]
'PaperToLean.Welford.scatter_correct' depends on axioms: [propext, Classical.choice, Quot.sound]

Run it

Install Lean with elan and Git. Download the project ZIP, extract it, and run these commands in that folder. The first command downloads the pinned Mathlib dependencies.

python3 prepare.py
lake env lean -DwarningAsError=true Welford.lean

A successful check prints the axiom dependencies and exits without errors.