ChungMG

Welford’s Algorithm

Welford’s algorithm is a numerically stable online algorithm for computing the running mean, variance, and covariance of a data stream in a single pass. Unlike the naïve two-pass or textbook one-pass formulas, Welford’s method avoids catastrophic cancellation by maintaining a running sum of squared deviations rather than a sum of squares — making it the standard choice for variance computation in production machine learning systems, including [[ResNet|Batch Normalization]], reinforcement learning, and distributed training.

---

1. Core Concept

1.1 The Problem: Numerically Unstable Variance

The textbook variance formula is:

$${\sigma }^2=\frac{1}{n}\sum _{i=1}^n{x}_i^2-{\left(\frac{1}{n}\sum _{i=1}^n{x}_i\right)}^2=\overline{x^2}-{\overline{x}}^2$$

Problem: When data values are large but the variance is small (e.g., $x_i\approx {10}^8$ with ${\sigma }^2\approx 1$ ), the two terms $\overline{x^2}$ and ${\overline{x}}^2$ are nearly equal large numbers. Their subtraction causes catastrophic cancellation, destroying significant digits.

[!NOTE] Catastrophic Cancellation Example
Consider $x=[{10}^8+1,{10}^8+2,{10}^8+3]$ . The true variance is $1.0$ . Using the naïve formula in float32, $\overline{x^2}\approx {10}^{16}$ and ${\overline{x}}^2\approx {10}^{16}$ , and their difference can yield zero or even negative variance due to floating-point rounding.

1.2 Welford’s Solution

Instead of accumulating $\sum x_i^2$ , Welford maintains a running sum of squared deviations from the current mean:

$$M_{2,n}=\sum _{i=1}^n(x_i-{\overline{x}}_n{)}^2$$

Since deviations $(x_i-{\overline{x}}_n)$ are small even when $x_i$ is large, the sum $M_{2,n}$ never suffers from cancellation.

1.3 Key Properties

• Single-pass: Each data point is processed exactly once — no need to store the full dataset
• Online: New data can be incorporated incrementally without recomputing from scratch
• Numerically stable: Avoids catastrophic cancellation in all regimes
• Memory efficient: Only $\mathcal{O}(1)$ state variables regardless of data size

---

2. Mathematical Derivation

2.1 Running Mean Update

Definition: The arithmetic mean of the first $n$ observations is:

$${\overline{x}}_n=\frac{1}{n}\sum _{i=1}^n{x}_i$$

Derivation of the incremental form:

We want to express ${\overline{x}}_n$ in terms of ${\overline{x}}_{n-1}$ without re-accessing previous data points. Start by separating the last observation:

$${\overline{x}}_n=\frac{1}{n}\sum _{i=1}^n{x}_i=\frac{1}{n}(\sum _{i=1}^{n-1}{x}_i+x_n)$$

Recognize that $\sum _{i=1}^{n-1}{x}_i=(n-1){\overline{x}}_{n-1}$ , so:

$${\overline{x}}_n=\frac{1}{n}((n-1){\overline{x}}_{n-1}+x_n)$$

Rewrite $(n-1){\overline{x}}_{n-1}$ as $n{\overline{x}}_{n-1}-{\overline{x}}_{n-1}$ :

$${\overline{x}}_n=\frac{1}{n}(n{\overline{x}}_{n-1}-{\overline{x}}_{n-1}+x_n)={\overline{x}}_{n-1}+\frac{x_n-{\overline{x}}_{n-1}}{n}$$ $$\boxed{{\displaystyle {\overline{x}}_n={\overline{x}}_{n-1}+\frac{x_n-{\overline{x}}_{n-1}}{n}}}$$

Interpretation: The new mean equals the old mean plus a correction term proportional to the prediction error $(x_n-{\overline{x}}_{n-1})$ , scaled by $1/n$ . As more data is observed, each new point has diminishing influence on the mean.

2.2 Running Variance Update (Core Result)

Definition: Define the sum of squared deviations from the current mean:

$$M_{2,n}=\sum _{i=1}^n(x_i-{\overline{x}}_n{)}^2$$

Theorem (Welford, 1962): The quantity $M_{2,n}$ satisfies the recurrence:

$$\boxed{{\displaystyle M_{2,n}=M_{2,n-1}+(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)}}$$

Complete Proof:

Step 1: Split the sum into the first $n-1$ terms and the $n$ -th term:

$$M_{2,n}=\sum _{i=1}^{n-1}(x_i-{\overline{x}}_n{)}^2+(x_n-{\overline{x}}_n{)}^2$$

Step 2: Express $(x_i-{\overline{x}}_n)$ in terms of $(x_i-{\overline{x}}_{n-1})$ . Define the shift $\delta ={\overline{x}}_n-{\overline{x}}_{n-1}=\frac{x_n-{\overline{x}}_{n-1}}{n}$ , so:

$$x_i-{\overline{x}}_n=(x_i-{\overline{x}}_{n-1})-\delta$$

Step 3: Expand the sum of squares of the first $n-1$ terms:

$$\sum _{i=1}^{n-1}(x_i-{\overline{x}}_n{)}^2=\sum _{i=1}^{n-1}((x_i-{\overline{x}}_{n-1})-\delta {)}^2$$ $$=\sum _{i=1}^{n-1}(x_i-{\overline{x}}_{n-1}{)}^2-2\delta \sum _{i=1}^{n-1}(x_i-{\overline{x}}_{n-1})+(n-1){\delta }^2$$

Step 4: The cross-term vanishes. By definition of ${\overline{x}}_{n-1}$ :

$$\sum _{i=1}^{n-1}(x_i-{\overline{x}}_{n-1})=\sum _{i=1}^{n-1}{x}_i-(n-1){\overline{x}}_{n-1}=0$$

Therefore:

$$\sum _{i=1}^{n-1}(x_i-{\overline{x}}_n{)}^2=M_{2,n-1}+(n-1){\delta }^2$$

Step 5: Substitute $\delta =\frac{x_n-{\overline{x}}_{n-1}}{n}$ :

$$(n-1){\delta }^2=(n-1)\cdot \frac{(x_n-{\overline{x}}_{n-1}{)}^2}{n^2}=\frac{n-1}{n^2}(x_n-{\overline{x}}_{n-1}{)}^2$$

Step 6: Handle the $n$ -th term. Since ${\overline{x}}_n={\overline{x}}_{n-1}+\delta$ :

$$x_n-{\overline{x}}_n=x_n-{\overline{x}}_{n-1}-\delta =(x_n-{\overline{x}}_{n-1})-\frac{x_n-{\overline{x}}_{n-1}}{n}=\frac{n-1}{n}(x_n-{\overline{x}}_{n-1})$$

So:

$$(x_n-{\overline{x}}_n{)}^2=\frac{(n-1{)}^2}{n^2}(x_n-{\overline{x}}_{n-1}{)}^2$$

Step 7: Combine all terms:

$$M_{2,n}=M_{2,n-1}+\frac{n-1}{n^2}(x_n-{\overline{x}}_{n-1}{)}^2+\frac{(n-1{)}^2}{n^2}(x_n-{\overline{x}}_{n-1}{)}^2$$

Factor out $(x_n-{\overline{x}}_{n-1}{)}^2$ :

$$M_{2,n}=M_{2,n-1}+(x_n-{\overline{x}}_{n-1}{)}^2\cdot \frac{(n-1)+(n-1{)}^2}{n^2}$$ $$=M_{2,n-1}+(x_n-{\overline{x}}_{n-1}{)}^2\cdot \frac{(n-1)(1+n-1)}{n^2}=M_{2,n-1}+(x_n-{\overline{x}}_{n-1}{)}^2\cdot \frac{(n-1)n}{n^2}$$ $$=M_{2,n-1}+\frac{n-1}{n}(x_n-{\overline{x}}_{n-1}{)}^2$$

Step 8: Rewrite into the symmetric Welford form. From Step 6, $(x_n-{\overline{x}}_n)=\frac{n-1}{n}(x_n-{\overline{x}}_{n-1})$ , therefore:

$$(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)=(x_n-{\overline{x}}_{n-1})\cdot \frac{n-1}{n}(x_n-{\overline{x}}_{n-1})=\frac{n-1}{n}(x_n-{\overline{x}}_{n-1}{)}^2$$

This is exactly the correction term from Step 7. Hence:

$$M_{2,n}=M_{2,n-1}+(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)◼$$

[!NOTE] Why This Is Stable
The product $(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)$ involves deviations from the mean, not raw values. These deviations are small even when $x_i$ is large, so no cancellation occurs when accumulating into $M_{2,n}$ . In floating-point arithmetic, both factors are computed via subtraction of quantities of similar magnitude, preserving relative precision.

[!NOTE] Non-Negativity Guarantee
Since ${\overline{x}}_n$ always lies between ${\overline{x}}_{n-1}$ and $x_n$ , both factors $(x_n-{\overline{x}}_{n-1})$ and $(x_n-{\overline{x}}_n)$ have the same sign (or one is zero). Therefore the correction term $(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)\ge 0$ always, ensuring $M_{2,n}\ge 0$ — the computed variance can never become negative, unlike the naïve formula.

2.3 Equivalence of Two Forms

The recurrence has two algebraically equivalent forms, both useful in practice:

Form A (product of deviations — numerically preferred):

$$M_{2,n}=M_{2,n-1}+(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)$$

Form B (scaled squared deviation — analytically convenient):

$$M_{2,n}=M_{2,n-1}+\frac{n-1}{n}(x_n-{\overline{x}}_{n-1}{)}^2$$

Form A requires two subtractions and one multiplication. Form B requires one subtraction, one squaring, and one division. In practice, Form A is preferred because it avoids the explicit division by $n$ in the correction term, reducing floating-point rounding.

2.4 Final Variance and Standard Deviation

From $M_{2,n}$ , the variance estimates are obtained by a single division:

$${\sigma }_{\text{population}}^2=\frac{M_{2,n}}{n},s_{\text{sample}}^2=\frac{M_{2,n}}{n-1}$$

where $s^2$ uses Bessel’s correction ( $n-1$ instead of $n$ ) to yield an unbiased estimator of the population variance.

$$\sigma =\sqrt{\frac{M_{2,n}}{n}},s=\sqrt{\frac{M_{2,n}}{n-1}}$$

---

3. Algorithm

3.1 Basic Welford Algorithm

class WelfordOnline:
    """Welford's online algorithm for running mean and variance."""
    
    def __init__(self):
        self.n = 0
        self.mean = 0.0
        self.M2 = 0.0  # Sum of squared deviations
    
    def update(self, x):
        """Process a new data point."""
        self.n += 1
        delta = x - self.mean
        self.mean += delta / self.n
        delta2 = x - self.mean
        self.M2 += delta * delta2
    
    @property
    def variance(self):
        """Population variance."""
        return self.M2 / self.n if self.n > 0 else 0.0
    
    @property
    def sample_variance(self):
        """Sample variance (Bessel's correction)."""
        return self.M2 / (self.n - 1) if self.n > 1 else 0.0
    
    @property
    def std(self):
        """Population standard deviation."""
        return math.sqrt(self.variance)

3.2 Batched Update (Vectorized)

For processing mini-batches in deep learning:

class WelfordBatched:
    """Welford's algorithm with batched (vectorized) updates."""
    
    def __init__(self, shape):
        self.n = 0
        self.mean = torch.zeros(shape)
        self.M2 = torch.zeros(shape)
    
    def update(self, x_batch):
        """
        Update statistics with a batch of data.
        
        Args:
            x_batch: Tensor of shape (batch_size, *shape)
        """
        batch_size = x_batch.shape[0]
        
        # Batch mean and variance
        batch_mean = x_batch.mean(dim=0)
        batch_M2 = ((x_batch - batch_mean) ** 2).sum(dim=0)
        
        # Parallel Welford merge
        delta = batch_mean - self.mean
        total_n = self.n + batch_size
        
        self.mean = self.mean + delta * (batch_size / total_n)
        self.M2 = (self.M2 + batch_M2 
                   + delta ** 2 * (self.n * batch_size / total_n))
        self.n = total_n

---

4. Parallel and Distributed Welford

4.1 Chan’s Parallel Algorithm

When data is split across multiple workers (distributed training), Welford statistics can be merged using Chan et al.'s formula (1983).

Setup: Given two disjoint subsets $A=\{a_1,\dots ,a_{n_A}\}$ and $B=\{b_1,\dots ,b_{n_B}\}$ with precomputed statistics $(n_A,{\overline{x}}_A,M_{2,A})$ and $(n_B,{\overline{x}}_B,M_{2,B})$ .

Goal: Compute $(n_{AB},{\overline{x}}_{AB},M_{2,AB})$ for $A\cup B$ without re-accessing the raw data.

Merged count and mean:

$$n_{AB}=n_A+n_B$$ $${\overline{x}}_{AB}=\frac{n_A{\overline{x}}_A+n_B{\overline{x}}_B}{n_A+n_B}$$

Derivation of merged $M_2$ :

Start from the definition:

$$M_{2,AB}=\sum _{i\in A}(a_i-{\overline{x}}_{AB}{)}^2+\sum _{j\in B}(b_j-{\overline{x}}_{AB}{)}^2$$

For the $A$ -part, write $a_i-{\overline{x}}_{AB}=(a_i-{\overline{x}}_A)+({\overline{x}}_A-{\overline{x}}_{AB})$ and expand:

$$\sum _{i\in A}(a_i-{\overline{x}}_{AB}{)}^2=\sum _{i\in A}(a_i-{\overline{x}}_A{)}^2+2({\overline{x}}_A-{\overline{x}}_{AB})\underset{=0}{\underbrace{\sum _{i\in A}(a_i-{\overline{x}}_A)}}+n_A({\overline{x}}_A-{\overline{x}}_{AB}{)}^2$$ $$=M_{2,A}+n_A({\overline{x}}_A-{\overline{x}}_{AB}{)}^2$$

By the same argument for $B$ :

$$\sum _{j\in B}(b_j-{\overline{x}}_{AB}{)}^2=M_{2,B}+n_B({\overline{x}}_B-{\overline{x}}_{AB}{)}^2$$

Therefore:

$$M_{2,AB}=M_{2,A}+M_{2,B}+n_A({\overline{x}}_A-{\overline{x}}_{AB}{)}^2+n_B({\overline{x}}_B-{\overline{x}}_{AB}{)}^2$$

Now simplify the correction terms. Define $\mathrm{\Delta }={\overline{x}}_A-{\overline{x}}_B$ . From the merged mean formula:

$${\overline{x}}_A-{\overline{x}}_{AB}={\overline{x}}_A-\frac{n_A{\overline{x}}_A+n_B{\overline{x}}_B}{n_A+n_B}=\frac{n_B({\overline{x}}_A-{\overline{x}}_B)}{n_A+n_B}=\frac{n_B}{n_{AB}}\mathrm{\Delta }$$ $${\overline{x}}_B-{\overline{x}}_{AB}=-\frac{n_A}{n_{AB}}\mathrm{\Delta }$$

Substitute:

$$n_A({\overline{x}}_A-{\overline{x}}_{AB}{)}^2+n_B({\overline{x}}_B-{\overline{x}}_{AB}{)}^2=n_A{\left(\frac{n_B\mathrm{\Delta }}{n_{AB}}\right)}^2+n_B{\left(\frac{n_A\mathrm{\Delta }}{n_{AB}}\right)}^2$$ $$=\frac{n_A{n}_B^2+n_B{n}_A^2}{n_{AB}^2}{\mathrm{\Delta }}^2=\frac{n_A{n}_B(n_B+n_A)}{n_{AB}^2}{\mathrm{\Delta }}^2=\frac{n_A{n}_B}{n_{AB}}{\mathrm{\Delta }}^2$$

Final result:

$$\boxed{{\displaystyle M_{2,AB}=M_{2,A}+M_{2,B}+\frac{n_A\cdot n_B}{n_A+n_B}({\overline{x}}_A-{\overline{x}}_B{)}^2}}$$

[!NOTE] Interpretation of the Correction Term
The extra term $\frac{n_A{n}_B}{n_A+n_B}{\mathrm{\Delta }}^2$ is the between-group sum of squares — it accounts for the variance introduced by the difference in group means. This is the same decomposition used in ANOVA: ${\text{SS}}_{\text{total}}={\text{SS}}_{\text{within}}+{\text{SS}}_{\text{between}}$ .

def merge_welford(stats_a, stats_b):
    """
    Merge two Welford statistics (Chan's parallel algorithm).
    
    Args:
        stats_a, stats_b: Tuples of (n, mean, M2)
    
    Returns:
        Merged (n, mean, M2)
    """
    n_a, mean_a, m2_a = stats_a
    n_b, mean_b, m2_b = stats_b
    
    n = n_a + n_b
    delta = mean_a - mean_b
    mean = (n_a * mean_a + n_b * mean_b) / n
    m2 = m2_a + m2_b + delta ** 2 * (n_a * n_b / n)
    
    return n, mean, m2

4.2 Distributed Training Application

In Distributed Data Parallel (DDP) training, each GPU computes local batch statistics, which are then merged via Chan’s algorithm:

GPU 0: (n₀, μ₀, M₂₀) ──┐
GPU 1: (n₁, μ₁, M₂₁) ──┤── Merge (Chan) ──→ Global (n, μ, M₂)
GPU 2: (n₂, μ₂, M₂₂) ──┤
GPU 3: (n₃, μ₃, M₂₃) ──┘

Advantage over all-reduce of sums:

• All-reduce of $\sum x$ and $\sum x^2$ : suffers from catastrophic cancellation
• Chan’s merge: numerically stable, identical results regardless of GPU count

---

5. Welford for Covariance and Correlation

5.1 Online Covariance — Derivation

Definition: The co-deviation sum for two streams $\{x_i\}$ and $\{y_i\}$ is:

$$C_n=\sum _{i=1}^n(x_i-{\overline{x}}_n)(y_i-{\overline{y}}_n)$$

Derivation of the recurrence: Split the sum as before:

$$C_n=\sum _{i=1}^{n-1}(x_i-{\overline{x}}_n)(y_i-{\overline{y}}_n)+(x_n-{\overline{x}}_n)(y_n-{\overline{y}}_n)$$

Using the identities from Section 2:

$$x_i-{\overline{x}}_n=(x_i-{\overline{x}}_{n-1})-{\delta }_x,{\delta }_x=\frac{x_n-{\overline{x}}_{n-1}}{n}$$ $$y_i-{\overline{y}}_n=(y_i-{\overline{y}}_{n-1})-{\delta }_y,{\delta }_y=\frac{y_n-{\overline{y}}_{n-1}}{n}$$

Expanding the product:

$$\sum _{i=1}^{n-1}(x_i-{\overline{x}}_n)(y_i-{\overline{y}}_n)=C_{n-1}-{\delta }_x\underset{=0}{\underbrace{\sum _{i=1}^{n-1}(y_i-{\overline{y}}_{n-1})}}-{\delta }_y\underset{=0}{\underbrace{\sum _{i=1}^{n-1}(x_i-{\overline{x}}_{n-1})}}+(n-1){\delta }_x{\delta }_y$$ $$=C_{n-1}+(n-1){\delta }_x{\delta }_y$$

For the $n$ -th term, from Section 2.2 Step 6:

$$x_n-{\overline{x}}_n=\frac{n-1}{n}(x_n-{\overline{x}}_{n-1}),y_n-{\overline{y}}_n=\frac{n-1}{n}(y_n-{\overline{y}}_{n-1})$$

Combining:

$$C_n=C_{n-1}+(n-1){\delta }_x{\delta }_y+\frac{(n-1{)}^2}{n^2}(x_n-{\overline{x}}_{n-1})(y_n-{\overline{y}}_{n-1})$$

Substituting ${\delta }_x{\delta }_y=\frac{(x_n-{\overline{x}}_{n-1})(y_n-{\overline{y}}_{n-1})}{n^2}$ :

$$C_n=C_{n-1}+\frac{n-1+(n-1{)}^2}{n^2}(x_n-{\overline{x}}_{n-1})(y_n-{\overline{y}}_{n-1})$$ $$=C_{n-1}+\frac{n-1}{n}(x_n-{\overline{x}}_{n-1})(y_n-{\overline{y}}_{n-1})$$

Rewriting in the symmetric Welford form (analogous to Section 2.3):

$$\boxed{{\displaystyle C_n=C_{n-1}+(x_n-{\overline{x}}_{n-1})(y_n-{\overline{y}}_n)}}$$

Note the asymmetry: the first factor uses the old mean ${\overline{x}}_{n-1}$ while the second uses the new mean ${\overline{y}}_n$ . This is not a typo — it is the algebraically correct form.

The covariance estimates are:

$$\text{Cov}(x,y)=\frac{C_n}{n}\text{(population)},\frac{C_n}{n-1}\text{(sample)}$$

class WelfordCovariance:
    """Online covariance via Welford's algorithm."""
    
    def __init__(self):
        self.n = 0
        self.mean_x = 0.0
        self.mean_y = 0.0
        self.C = 0.0  # Co-deviation sum
    
    def update(self, x, y):
        self.n += 1
        dx = x - self.mean_x
        self.mean_x += dx / self.n
        dy = y - self.mean_y  # Note: uses OLD mean_y
        self.mean_y += (y - self.mean_y) / self.n
        self.C += dx * (y - self.mean_y)  # NEW mean_y
    
    @property
    def covariance(self):
        return self.C / self.n if self.n > 0 else 0.0
    
    @property
    def correlation(self):
        # Requires also tracking M2_x and M2_y
        return self.C / math.sqrt(self.M2_x * self.M2_y) if self.n > 0 else 0.0

5.2 Online Covariance Matrix

For $d$ -dimensional data ${\mathbf{x}}_i\in {\mathbb{R}}^d$ , the full covariance matrix generalizes the scalar case via an outer product:

$${\mathbf{C}}_n={\mathbf{C}}_{n-1}+({\mathbf{x}}_n-{\overline{\mathbf{x}}}_{n-1})({\mathbf{x}}_n-{\overline{\mathbf{x}}}_n{)}^{\mathrm{\top }}$$

Structure: At each step, this is a rank-1 update — the outer product of two $d$ -dimensional vectors. The computational cost is $\mathcal{O}(d^2)$ per observation, and the storage is $\mathcal{O}(d^2)$ for the full covariance matrix.

The correlation matrix is then:

$${\mathbf{R}}_n=\text{diag}({\mathbf{C}}_n{)}^{-1/2}{\mathbf{C}}_n\text{diag}({\mathbf{C}}_n{)}^{-1/2}$$

---

6. Weighted and Exponentially-Weighted Variants

6.1 Weighted Welford

For weighted observations $(x_i,w_i)$ with positive weights, define:

$$W_n=\sum _{i=1}^n{w}_i,{\overline{x}}_n=\frac{1}{W_n}\sum _{i=1}^n{w}_i{x}_i,M_{2,n}=\sum _{i=1}^n{w}_i(x_i-{\overline{x}}_n{)}^2$$

Derivation of the weighted mean update:

$${\overline{x}}_n=\frac{W_{n-1}{\overline{x}}_{n-1}+w_n{x}_n}{W_n}=\frac{W_{n-1}{\overline{x}}_{n-1}+w_n{x}_n}{W_{n-1}+w_n}$$

Rewrite as:

$${\overline{x}}_n={\overline{x}}_{n-1}+\frac{w_n}{W_n}(x_n-{\overline{x}}_{n-1})$$

Derivation of the weighted $M_2$ update: Following the same algebraic strategy as Section 2.2, define ${\delta }_w={\overline{x}}_n-{\overline{x}}_{n-1}=\frac{w_n}{W_n}(x_n-{\overline{x}}_{n-1})$ :

$$M_{2,n}=\sum _{i=1}^{n-1}{w}_i(x_i-{\overline{x}}_n{)}^2+w_n(x_n-{\overline{x}}_n{)}^2$$

Expanding $(x_i-{\overline{x}}_n)=(x_i-{\overline{x}}_{n-1})-{\delta }_w$ and using $\sum _{i=1}^{n-1}{w}_i(x_i-{\overline{x}}_{n-1})=0$ :

$$\sum _{i=1}^{n-1}{w}_i(x_i-{\overline{x}}_n{)}^2=M_{2,n-1}+W_{n-1}{\delta }_w^2$$

For the $n$ -th term: $x_n-{\overline{x}}_n=(x_n-{\overline{x}}_{n-1})-{\delta }_w=(1-\frac{w_n}{W_n})(x_n-{\overline{x}}_{n-1})=\frac{W_{n-1}}{W_n}(x_n-{\overline{x}}_{n-1})$ , so:

$$w_n(x_n-{\overline{x}}_n{)}^2=w_n\cdot \frac{W_{n-1}^2}{W_n^2}(x_n-{\overline{x}}_{n-1}{)}^2$$

Combining and simplifying:

$$W_{n-1}{\delta }_w^2+w_n(x_n-{\overline{x}}_n{)}^2=\frac{W_{n-1}{w}_n^2}{W_n^2}(x_n-{\overline{x}}_{n-1}{)}^2+\frac{w_n{W}_{n-1}^2}{W_n^2}(x_n-{\overline{x}}_{n-1}{)}^2$$ $$=\frac{w_n{W}_{n-1}(w_n+W_{n-1})}{W_n^2}(x_n-{\overline{x}}_{n-1}{)}^2=\frac{w_n{W}_{n-1}}{W_n}(x_n-{\overline{x}}_{n-1}{)}^2$$

Rewriting in the symmetric form:

$$\boxed{{\displaystyle M_{2,n}=M_{2,n-1}+w_n(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)}}$$

This reduces to the standard (unweighted) Welford recurrence when $w_i=1$ for all $i$ .

6.2 Exponentially Weighted Moving Variance (EWMA)

For non-stationary data (e.g., loss tracking, reward monitoring), we use exponential forgetting with decay rate $\alpha \in (0,1)$ :

$${\overline{x}}_t=(1-\alpha ){\overline{x}}_{t-1}+\alpha x_t$$

The exponentially weighted variance update (West, 1979):

$$v_t=(1-\alpha )(v_{t-1}+\alpha (x_t-{\overline{x}}_{t-1}{)}^2)$$

Derivation: This follows from the weighted Welford formula with exponentially decaying weights $w_t=(1-\alpha {)}^{T-t}$ , which yields effective $W_t\to \frac{1}{\alpha }$ as $T\to \mathrm{\infty }$ . Substituting into the weighted recurrence and normalizing gives the formula above.

[!NOTE] Connection to Adam Optimizer
The Adam optimizer maintains an exponentially weighted second moment $v_t={\beta }_2{v}_{t-1}+(1-{\beta }_2)g_t^2$ . While structurally similar, Adam tracks $\mathbb{E}[g^2]$ (raw second moment), not $\text{Var}(g)$ (centered variance). The EWMA variance formula tracks the centered second moment, which is a strictly more informative quantity.

---

7. Applications in Deep Learning

7.1 [[ResNet|Batch Normalization]]

BatchNorm computes per-channel mean and variance during training:

$$\hat{x}=\frac{x-{\mu }_B}{\sqrt{{\sigma }_B^2+ϵ}}$$

While standard BatchNorm uses the naïve formula (batch sizes are typically small enough), Welford’s algorithm is used in production implementations (PyTorch, TensorFlow) for numerical robustness, especially with:

• Small batch sizes (e.g., batch size 1-4 in detection/segmentation)
• Float16 / BFloat16 training (limited mantissa bits)
• Running statistics accumulation across many batches

# PyTorch internally uses Welford for running stats
# Conceptual equivalent:
class RunningBatchNorm:
    def __init__(self, num_features, momentum=0.1):
        self.mean = torch.zeros(num_features)
        self.var = torch.ones(num_features)
        self.welford = WelfordBatched(shape=(num_features,))
        self.momentum = momentum
    
    def forward(self, x):
        if self.training:
            # Batch statistics (naïve is fine for single batch)
            batch_mean = x.mean(dim=[0, 2, 3])
            batch_var = x.var(dim=[0, 2, 3], unbiased=False)
            
            # Update running statistics (Welford-style)
            self.mean = (1 - self.momentum) * self.mean + self.momentum * batch_mean
            self.var = (1 - self.momentum) * self.var + self.momentum * batch_var
            
            return (x - batch_mean[None, :, None, None]) / \
                   torch.sqrt(batch_var[None, :, None, None] + 1e-5)
        else:
            return (x - self.mean[None, :, None, None]) / \
                   torch.sqrt(self.var[None, :, None, None] + 1e-5)

7.2 Layer Normalization & RMSNorm

In [[Vision Transformer (ViT)|Transformers]], LayerNorm requires computing variance over the feature dimension:

$$\text{LayerNorm}(x)=\frac{x-\mu }{\sqrt{{\sigma }^2+ϵ}}\cdot \gamma +\beta$$

For large hidden dimensions ( $d=4096+$ ), Welford ensures stable normalization even in mixed precision.

7.3 Reinforcement Learning: Advantage Estimation

In PPO and other policy gradient methods, advantage normalization requires running statistics of returns:

$${\hat{A}}_t=\frac{A_t-{\mu }_A}{{\sigma }_A+ϵ}$$

Since RL agents interact with the environment over millions of timesteps, Welford’s online algorithm is the natural choice — storing all returns would be prohibitively expensive.

7.4 Gradient Monitoring and Clipping

Tracking gradient statistics for adaptive clipping:

$$\text{grad\_norm\_clip}=\text{clip}(\frac{\parallel g\parallel }{{\sigma }_g},\text{max\_val})$$

Welford tracks ${\sigma }_g$ online without storing gradient history.

7.5 Loss Logging and Early Stopping

Monitoring training loss statistics:

• Running mean of loss for smoothed curves
• Running variance for detecting training instability (loss spikes)
• Running covariance between loss and learning rate for adaptive scheduling

---

8. Comparison with Other Methods

8.1 Variance Computation Methods

Method	Passes	Online	Numerically Stable	Memory	Parallel
Two-pass (subtract mean, then square)	2	❌	✅ Stable	$\mathcal{O}(n)$	❌
Naïve one-pass ( $\overline{x^2}-{\overline{x}}^2$ )	1	✅	❌ Catastrophic cancellation	$\mathcal{O}(1)$	✅
Welford	1	✅	✅ Stable	$\mathcal{O}(1)$	✅ (Chan)
Youngs & Cramer	1	✅	✅ Stable	$\mathcal{O}(1)$	✅
Kahan summation + naïve	1	✅	⚠️ Better but not guaranteed	$\mathcal{O}(1)$	❌

8.2 When to Use What

Scenario	Recommended Method
Small dataset, fits in memory	Two-pass (simplest, stable)
Streaming data, single pass	Welford
Distributed / multi-GPU	Chan’s parallel Welford
Mixed precision (float16)	Welford (essential)
High-dimensional covariance	Welford with rank-1 updates

---

9. Numerical Analysis

9.1 Error Bounds

For $n$ observations with machine precision ${ϵ}_{\text{mach}}$ :

Method	Relative Error Bound
Naïve	$\mathcal{O}(n{\kappa }^2{ϵ}_{\text{mach}})$ where $\kappa =\frac{max|x_i|}{\sigma }$
Welford	$\mathcal{O}(n{ϵ}_{\text{mach}})$ — independent of data scale

Key insight: Welford’s error does not depend on the condition number $\kappa$ , making it robust regardless of data magnitude.

9.2 Catastrophic Cancellation Demonstration

import numpy as np

# Data: large mean, small variance
data = np.array([1e8 + 1, 1e8 + 2, 1e8 + 3], dtype=np.float32)

# Naïve formula
mean_sq = np.mean(data ** 2)
sq_mean = np.mean(data) ** 2
naive_var = mean_sq - sq_mean
print(f"Naïve variance:   {naive_var:.6f}")  # Often 0.0 or negative!

# Welford
welford = WelfordOnline()
for x in data:
    welford.update(float(x))
print(f"Welford variance: {welford.variance:.6f}")  # Correct: 0.666667

9.3 Float16 Considerations

In mixed-precision training (AMP), variance computation is especially fragile:

Precision	Mantissa Bits	Cancellation Threshold
float32	23 bits	$\kappa >{10}^3$ problematic
float16	10 bits	$\kappa >{10}^1$ problematic
bfloat16	7 bits	$\kappa >{10}^0$ problematic

Welford’s algorithm is essential for float16/bfloat16 training — without it, normalization layers produce NaN or Inf values.

---

10. Core Formula Cards

[!QUOTE] Welford Running Mean

$${\overline{x}}_n={\overline{x}}_{n-1}+\frac{x_n-{\overline{x}}_{n-1}}{n}$$

[!QUOTE] Welford Running Variance (Core Recurrence)

$$M_{2,n}=M_{2,n-1}+(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)$$

[!QUOTE] Variance from $M_2$

$${\sigma }^2=\frac{M_{2,n}}{n},s^2=\frac{M_{2,n}}{n-1}$$

[!QUOTE] Chan’s Parallel Merge

$$M_{2,A\cup B}=M_{2,A}+M_{2,B}+\frac{n_A{n}_B}{n_A+n_B}({\overline{x}}_A-{\overline{x}}_B{)}^2$$

[!QUOTE] Online Covariance

$$C_n=C_{n-1}+(x_n-{\overline{x}}_{n-1})(y_n-{\overline{y}}_n)$$

[!QUOTE] Weighted Welford Update

$$M_{2,n}=M_{2,n-1}+w_n(x_n-{\overline{x}}_{n-1})(x_n-{\overline{x}}_n)$$

[!QUOTE] Exponentially Weighted Variance

$$v_t=(1-\alpha )(v_{t-1}+\alpha (x_t-{\overline{x}}_{t-1}{)}^2)$$

---

11. Summary

Aspect	Description
Core idea	Accumulate squared deviations from running mean, not raw squares
Key advantage	Numerically stable — immune to catastrophic cancellation
Complexity	$\mathcal{O}(1)$ time and space per update
Parallel extension	Chan’s merge formula for distributed computation
Generalization	Covariance, correlation, weighted, exponentially-weighted
Role in DL	BatchNorm running stats, LayerNorm, RL advantage normalization, gradient monitoring
When essential	Mixed precision (float16/bfloat16), large-scale data, distributed training

Welford’s algorithm exemplifies a fundamental principle in numerical computing: reformulate to avoid subtracting large, nearly-equal quantities. By tracking deviations rather than raw sums, it achieves stability that the naïve formula cannot — a lesson that recurs throughout machine learning, from [[ResNet|BatchNorm]] to the [[Score Function|score function]] estimation in diffusion models.

---

Related Concepts

• [[ResNet]]
• [[U-Net]]
• [[Vision Transformer (ViT)]]
• [[Diffusion Model]]
• [[Score Function]]

---

Dataview Query

LIST
FROM #welford OR #online_algorithm OR #numerical_stability
SORT file.ctime DESC

---

References

• Paper: Note on a Method for Calculating Corrected Sums of Squares and Products (Welford, Technometrics 1962)
• Paper: Updating Mean and Variance Estimates: An Improved Method (West, Communications of the ACM 1979)
• Paper: Algorithms for Computing the Sample Variance: Analysis and Recommendations (Chan, Golub, LeVeque, The American Statistician 1983)
• Paper: Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift (Ioffe & Szegedy, ICML 2015)
• Blog: Accurately computing running variance — John D. Cook
• Docs: PyTorch torch.nn.BatchNorm2d implementation notes
• Code: cpython/Modules/_statisticsmodule.c — Python’s statistics.variance uses Welford internally