welford's algorithm and why i never store the samples

you want to know if the error rate just spiked. that’s a z-score: how many standard deviations is right now away from what normal looked like. simple formula, (x - mean) / stddev, you learned it before you could legally drink.

except a standard deviation needs a mean, and a mean needs all the numbers, and the numbers are a log file the size of a small moon. so where, exactly, are you keeping them?

last time in this series i ported the drain algorithm to rust to collapse forty thousand log lines into a dozen templates, and i closed by promising the sequel would be about flagging error-rate spikes without storing a single sample. here it is.

quick refresher on what this thing is: logmole ↗ is a little cli that points at a log file, auto-detects the format, and flags anomalies. it streams. it never loads the whole file into memory — mmap in, records out, one pass. that constraint is the whole personality of the tool, and it’s why the anomaly detector can’t do the obvious thing.

the z-score, and the two obvious ways to get it wrong

the naive approach is to collect your samples into a Vec<f64>, and at the end compute the mean, then loop again for the variance. two passes, and O(n) memory where n is every error rate you ever measured. for a cli triaging a bounded file you could maybe get away with it. for something that’s supposed to chew through gigabytes without breathing hard, keeping every sample around is exactly the thing i’m trying not to do.

so people reach for the “streaming” trick they half-remember: keep a running sum and a running sum-of-squares, and recover the variance from E[x²] - (E[x])². constant memory, one pass, done.

it’s also a numerical landmine. when your values are large and close together — which error rates near some steady baseline absolutely are — E[x²] and (E[x])² are two big, nearly-equal floating-point numbers, and subtracting them is catastrophic cancellation ↗ : the significant digits annihilate each other and you’re left holding rounding noise. you can compute a negative variance this way. i have seen it happen. it is deeply stupid to take the square root of a negative number because your arithmetic betrayed you.

there’s a better way, and it’s from 1962.

welford’s online variance

b. p. welford published the fix in a two-page note in technometrics — “note on a method for calculating corrected sums of squares and products” (1962, vol. 4, pp. 419–420). knuth reprints it in taocp vol. 2, and the running-mean / M2 recurrence ↗ is the canonical write-up. it keeps three numbers and updates them one sample at a time:

logmole-core/src/analysis/anomaly.rs

/// Rolling statistics for z-score computation using Welford's online algorithm.
/// O(1) per push and per z-score query — no stored values, constant memory.
struct RollingStats {
    count: u64,
    mean: f64,
    m2: f64,
}

that’s the entire state. a count, the running mean, and m2 — the running sum of squared deviations from the current mean. no samples. the update per new value:

logmole-core/src/analysis/anomaly.rs

fn push(&mut self, value: f64) {
    self.count += 1;
    let delta = value - self.mean;
    let count_f = self.count as f64;
    self.mean += delta / count_f;
    let delta2 = value - self.mean;
    self.m2 += delta * delta2;
}

the trick is those two deltas. delta is the gap before we update the mean; delta2 is the gap after. m2 accumulates their product. crucially you never subtract two big numbers — every term is a deviation, already small, so there’s nothing to catastrophically cancel. that’s the numerical stability, for free, as a side effect of the update order. variance is just m2 / count whenever you want it, and the samples are long gone.

O(1) memory, O(1) per push, single pass, and it doesn’t lie to you. this is the boring correct answer, and i love it.

the near-zero-stddev guard

there’s one edge case welford doesn’t cover because it’s not a numerical problem, it’s a semantic one. if every value in your baseline is identical — a service that logged exactly 0.00 errors for the whole warmup window — then the standard deviation is zero, and (x - mean) / stddev is a division by zero. even a whisper of an error afterwards is infinitely many standard deviations away.

so z_score guards it:

logmole-core/src/analysis/anomaly.rs

fn z_score(&self, value: f64) -> Option<f64> {
    if self.count < 2 {
        return None;
    }

    let variance = self.m2 / self.count as f64;
    let stddev = variance.sqrt();

    if stddev < 1e-10 {
        // All baseline values nearly identical.
        // Any meaningful deviation is anomalous.
        if (value - self.mean).abs() > 1e-10 {
            return Some((value - self.mean).signum() * 10.0);
        }
        return Some(0.0);
    }

    Some((value - self.mean) / stddev)
}

fewer than two samples and there’s no variance to speak of, so it returns None — the caller treats that as “not enough data, don’t flag.” a flat baseline with a real deviation returns a hard-coded ±10.0, comfortably past any sane threshold, instead of inf. and a flat baseline that stays flat returns 0.0. that 10.0 is a magic number and i’ll cop to it in the honest section.

baseline vs detection: the 80/20 split

a z-score is only meaningful against a baseline of what normal looked like. logmole doesn’t have a separate training run — it has one file — so it splits that file in place: the first ~80% builds the baseline, the rest gets scored against it.

logmole-core/src/analysis/anomaly.rs

// Check if we're still in baseline phase
if !self.baseline_complete {
    let baseline_threshold = (self.estimated_total as f64 * self.baseline_ratio) as u64;
    if self.total_lines <= baseline_threshold || self.total_lines < self.min_lines {
        // Still in baseline — collect statistics
        self.baseline_lines += 1;
        if let Some(rate) = error_rate {
            self.error_rate_stats.push(rate);
        }
        if template_id > 0 {
            *self.template_baseline.entry(template_id).or_insert(0) += 1;
        }
        return AnomalyResult::normal();
    }
    // Transition to detection phase
    self.baseline_complete = true;
    self.baseline_max_template_id =
        Some(self.template_baseline.keys().max().copied().unwrap_or(0));
}

during baseline, every line just feeds the stats — pushes its error rate into the RollingStats, bumps its template’s count — and returns normal() unconditionally. nothing is ever flagged in the first 80%. the moment we cross the threshold, baseline_complete flips, and we snapshot the highest template id we saw. that snapshot is the whole trick behind novel-pattern detection: any template id higher than it is, by construction, a shape that didn’t exist while things were healthy. (that keys().max() is the drain max_id() idea from the last post, cashed in.)

baseline_ratio defaults to 0.8 and min_lines to 100, so tiny files skip detection entirely rather than build a “baseline” out of nine lines and start crying wolf.

three detectors, combined by max

once we’re in the detection phase, three independent detectors run, each producing a score in [0.0, 1.0]:

logmole-core/src/analysis/anomaly.rs

//! Anomaly detection: error rate spikes, novel patterns, frequency deviations.
//!
//! Three detection methods, each producing an independent score (0.0-1.0).
//! Final anomaly score = max of individual scores.

detector one, the z-score spike — this is what all the welford machinery was for:

logmole-core/src/analysis/anomaly.rs

// 1. Error rate spike detection
if let Some(rate) = error_rate {
    if let Some(z) = self.error_rate_stats.z_score(rate)
        && z > self.z_threshold
    {
        let score = (z / (self.z_threshold * 2.0)).min(1.0);
        reasons.push(AnomalyReason {
            kind: AnomalyKind::ErrorRateSpike,
            score,
            description: format!("error rate spike (z={z:.1})"),
        });
    }
    // Keep updating rolling stats
    self.error_rate_stats.push(rate);
}

past z_threshold (default 3.0) it fires, and the raw z gets squashed into a score: twice the threshold — a z of 6 — saturates at 1.0. note it keeps push-ing during detection too, so the baseline gently drifts. and let ... && ... chained let-conditions, stable since edition 2024, one of the small quality-of-life wins from why i left go for rust .

detector two, the novel pattern — with an anti-noise rule baked in:

logmole-core/src/analysis/anomaly.rs

// 2. Novel pattern detection
if template_id > 0
    && let Some(baseline_max) = self.baseline_max_template_id
{
    if template_id > baseline_max {
        let novel_count = self.template_detection.entry(template_id).or_insert(0);
        *novel_count += 1;

        // Only flag if seen >= 3 times (avoid one-off noise)
        if *novel_count >= 3 {
            reasons.push(AnomalyReason {
                kind: AnomalyKind::NovelPattern,
                score: 0.8,
                description: format!(
                    "novel pattern (not in baseline, seen {} times)",
                    *novel_count
                ),
            });
        }
    } else {
        // ...
    }
}

a brand-new template id means a log shape we never saw while healthy. but a single weird line is usually just noise — a one-off stack trace, a truncated write. so a novel pattern has to show up at least three times before it earns a flag. that threshold is the difference between “useful signal” and “screams at every unusual line,” and it cost me exactly one >= 3.

detector three, frequency deviation, i’ll spare you the full block — it compares each known template’s rate during detection against its baseline rate, and flags when the ratio falls outside [0.33, 3.0] (a template firing 3× more or less than it used to). a login line that suddenly goes quiet is as interesting as an error line that goes loud.

and then the combine step, which is almost insultingly simple:

logmole-core/src/analysis/anomaly.rs

// Final score = max of individual scores
let score = reasons.iter().map(|r| r.score).fold(0.0_f64, f64::max);

AnomalyResult { score, reasons }

no weighting, no averaging, no bayesian anything. the line’s anomaly score is just the loudest of the three detectors. if the z-score is calm but a novel pattern showed up, you get 0.8 and the reason attached. fold(0.0, f64::max) because f64 isn’t Ord and i’m not going to pretend nan doesn’t exist.

an aside: percentiles never store the samples either

same philosophy shows up one file over. “what’s my p99 latency” sounds like it needs every latency sorted in memory. it doesn’t — hdrhistogram ↗ buckets values into a fixed set of ranges at a configured precision and reads quantiles straight off the buckets:

logmole-core/src/analysis/percentile.rs

/// Create a new percentile tracker.
/// Range: 1 to 3,600,000 ms (1ms to 1 hour) with 3 significant digits.
pub fn new() -> Self {
    Self {
        histogram: Histogram::new_with_bounds(1, 3_600_000, 3)
            .expect("valid histogram parameters"),
        count: 0,
        field_name: None,
    }
}

bounded memory, three significant digits from 1ms to an hour, and value_at_quantile(0.99) is O(1). i didn’t write the clever part — gil tene did, i just linked the crate — but it’s the same commitment as welford: decide up front that you will never keep the raw data, and design backward from there.

the honest part: what i cut

the rule in this journal is that i don’t get to pretend. so:

• the z-score assumes error rates are roughly normal. they aren’t, really — error rates are bounded at zero and often bursty and skewed. a z-score over a non-normal distribution still works as a “this is unusually far from typical” heuristic, but the 3.0 threshold doesn’t map to the clean “0.1% of samples” you’d get from an actual gaussian. it’s a smoke alarm, not a p-value.
• max-combining is crude. a line that trips all three detectors scores identically to a line that trips only the worst one. that’s clearly wrong if you think three weak signals should compound — but the moment you start weighting and summing, you’re inventing coefficients you can’t defend, and i’d rather ship an honest max than a made-up formula.
• the 80/20 split is a heuristic wearing a lab coat. it assumes the first 80% of the file is the healthy part. if your log starts mid-incident, logmole faithfully learns the incident as “normal” and flags the recovery. for post-hoc triage of a file you already suspect, that’s usually fine. it is not a monitoring system.
• it’s population variance, not sample. m2 / count, not m2 / (count - 1) — no bessel’s correction. over a baseline of thousands of lines the difference is a rounding error, and i’d rather not divide by zero at count == 1.
• that 10.0 is a magic number. the flat-baseline sentinel is “big enough to always flag” and nothing more principled than that.

none of these are bugs. they’re scope. the job was to surface “this file has a spike, a new shape, or a rate that moved” without holding the file in ram, and it does that.

lessons learned

• the boring 1962 answer is usually the right one. welford’s online variance is older than most languages, needs three fields, and quietly dodges an entire class of floating-point bug that the “obvious” streaming formula walks straight into. reach for the known-stable recurrence before you reach for cleverness.
• your memory budget is a design input, not an afterthought. “i will never store the samples” isn’t a constraint that limited the anomaly detector — it’s the decision that shaped it, straight through to picking hdrhistogram for percentiles. commit to it first and the rest follows.
• write down the seams. the >= 3 anti-noise rule, the 10.0 sentinel, the 80/20 split — each is a judgment call, and a one-line comment next to each turns “why is this here?” into “ah, that’s the trade-off.”

next in the logmole series: the query dsl — a hand-rolled lexer and recursive-descent parser for status >= 500 AND latency_ms > 2s, and why i didn’t reach for a parser generator.

the code is on github ↗ if you want to poke at it. it’s a work in progress, like everything here.