Welford's streaming variance — numerically stable computation

Variance looks simple. Sum of squares minus square of sum, divided by N. But try computing that with values around 10^9 and a variance of 1.

The naive formula computes E[X²] - E[X]². For large values, both terms are enormous — on the order of 10^18 — and their difference is tiny. You lose every significant digit to floating point cancellation. The result can be negative. It can be NaN. It’s never correct.

Welford’s algorithm avoids this by updating the mean incrementally. Instead of recomputing everything from scratch, it adjusts three running values per sample: count, mean, and M2 (sum of squared differences from the mean).

The accumulator

WelfordAccumulator is the core. It’s a struct with three fields and a push method that does exactly three updates:

use tflo_ops::primitives::WelfordAccumulator;

let mut acc = WelfordAccumulator::new();

acc.push(2.0);
acc.push(4.0);
acc.push(4.0);
acc.push(4.0);
acc.push(5.0);
acc.push(5.0);
acc.push(7.0);
acc.push(9.0);

assert!((acc.mean() - 5.0).abs() < 1e-10);
assert!((acc.population_variance() - 4.0).abs() < 1e-10);

Each push does:

1. Increment count: n += 1
2. Update mean: mean += (x - mean) / n
3. Update M2: M2 += (x - mean_old) * (x - mean_new)

No squaring of large numbers. No catastrophic cancellation. The algorithm is stable for millions of samples.

Querying the results

let count = acc.count();              // u64 — number of values pushed
let mean = acc.mean();                // f64 — arithmetic mean
let pop_var = acc.population_variance(); // f64 — variance (divide by N)
let sample_var = acc.sample_variance();  // f64 — variance (divide by N-1)
let pop_std = acc.population_std();      // f64 — standard deviation
let sample_std = acc.sample_std();       // f64 — sample standard deviation

All variance and std methods return f64::NAN if there are fewer than 2 values. Mean returns f64::NAN if the accumulator is empty. You won’t accidentally use uninitialized statistics.

Sliding windows

For streaming use, you usually don’t want an infinite accumulator — you want the variance over the last N seconds or last N samples. That’s what WelfordWindow handles:

use tflo_ops::primitives::WelfordWindow;
use std::time::Duration;

let mut window = WelfordWindow::new(Duration::from_secs(300)); // 5 minute window

window.push(timestamp_ms, value);
window.push(timestamp_ms, value);
// ...

let mean = window.mean();
let variance = window.variance();
let std = window.std();

The window evicts old values when they fall outside the time range. It uses Welford’s remove() operation to subtract old values from the running statistics. And it periodically recomputes from scratch to prevent drift from accumulating across thousands of remove operations.

Higher moments

Variance is the second moment. But distributions have shape — skewness and kurtosis tell you about asymmetry and tail weight. MomentsCountWindow extends the same incremental approach to higher moments:

use tflo_ops::primitives::MomentsCountWindow;

let mut window = MomentsCountWindow::new(100);

for value in samples {
    window.push(value);
}

let skew = window.skewness();   // Third moment — positive means right-tailed
let kurt = window.kurtosis();   // Fourth moment — positive means heavy tails

Skewness needs at least 3 values. Kurtosis needs at least 4. Both return f64::NAN when there’s insufficient data.

Correlation

Pearson correlation measures whether two series move together. CorrelationCountWindow uses a Welford-style incremental algorithm that tracks joint statistics — sum of x, sum of y, sum of x², sum of y², sum of xy:

use tflo_ops::primitives::CorrelationCountWindow;

let mut window = CorrelationCountWindow::new(20);

for (x, y) in paired_samples {
    window.push(x, y);
}

let r = window.correlation();   // Pearson's r, from -1 to 1
let cov = window.covariance();  // Covariance
let beta = window.beta();       // Regression slope
let alpha = window.alpha();     // Regression intercept

The window maintains these running sums and recomputes correlation from them on demand. Same numerical stability as Welford, but for the joint distribution of two variables.

Available as graph nodes

You don’t have to use the raw structs. These are all available as graph nodes via Comp<R, f64>:

use tflo_core::prelude::*;
use tflo_ops::prelude::*;

let results: Vec<f64> = events.into_iter()
    .tflo(|t| {
        t.timestamp(|x| x.ts);
        let value = t.prop(|x| x.value);
        let other = t.prop(|x| x.other);

        // Rolling variance over 20 samples
        let var = value.variance(20usize);

        // Rolling skewness and kurtosis
        let skew = value.skewness(20usize);
        let kurt = value.kurtosis(20usize);

        // Rolling correlation between the two signals
        let corr = value.correlation(&other, 20usize);

        (var, skew, kurt, corr)
    })
    .collect();

Time-based versions also exist: value.variance(20_u64.secs()), value.skewness(30_u64.secs()), value.kurtosis(30_u64.secs()), value.correlation(&other, 5_u64.mins()).

Welford’s algorithm is a textbook solution to a textbook problem. It’s been known since 1962. Most libraries just don’t bother implementing it. tflo does — because when your pipeline runs for days against millions of samples, the difference between stable variance and catastrophic cancellation is the difference between a working system and a silent failure.

Source code

The full working example is available on GitHub: tflo-examples/examples/welford-streaming-variance

— Matt