Numerically stable calculations

The textbook formula computes variance as

E[x²] − (E[x])²

Each term is computed with a small relative error. When the variance is small compared with the mean squared (prices near 4000 moving by a few ticks, for example), the two terms are large and nearly equal. Subtracting them cancels the leading digits that agree and leaves mostly the rounding errors: the result can have no correct digits at all, and can be negative.

That is catastrophic cancellation: subtraction of nearly equal quantities magnifies their relative errors. The fix is to rearrange the computation so that the subtraction happens before the squaring, on quantities of the size of the answer.

Welford's method updates the mean and the sum of squared deviations from the current mean, one observation at a time:

struct RunningStats {
    std::int64_t n = 0;
    double mean = 0, m2 = 0;            // m2 = sum of squared deviations from the mean
    void add(double x) {
        ++n;
        double d = x - mean;            // deviation from the old mean
        mean += d / n;
        m2 += d * (x - mean);           // deviation from the old mean times deviation from the new
    }
    double variance() const { return n > 1 ? m2 / (n - 1) : 0.0; }
};

Every subtraction is between numbers close to each other only when the deviation is small, so the result keeps its precision. m2 can never become negative. The same idea extends to covariance, and to combining statistics computed on separate chunks or threads.

Adding many numbers into one running total loses the low-order bits of each small term once the total is large: adding 1 to 2⁵³ has no effect at all in double precision. Over millions of additions, the error grows. Two fixes:

Kahan (compensated) summation keeps a running correction for the bits lost at each step:

double sum = 0, c = 0;                 // c accumulates the lost low-order bits
for (double x : xs) {
    double y = x - c;
    double t = sum + y;
    c = (t - sum) - y;                 // what was lost in sum + y
    sum = t;
}

The error becomes independent of the number of terms, at the cost of four operations per addition instead of one. Note that -ffast-math may algebraically simplify c to zero and silently remove the correction.

Pairwise or blocked summation: add in a tree, or sum blocks and then sum the block totals. Error grows with the logarithm of the count instead of linearly, and it vectorises well.

Better still, where the quantities are really integers, such as cash in ticks times quantities, keep them as integers and the problem disappears.

The common update ewma = alpha * x + (1 - alpha) * ewma is stable. Two practical issues:

• Precision of alpha: for very small alpha the update ewma += alpha * (x - ewma) is better conditioned and cheaper.
• Time decay: for irregularly spaced events, alpha should depend on the time since the previous update, alpha = 1 - exp(-dt / tau), computed with -std::expm1(-dt / tau) for accuracy when dt is small. And a value decaying for a long time without input can drift into denormals (previous lesson).

Floating-point addition is not associative: (a + b) + c and a + (b + c) can differ. Consequences:

• A parallel or vectorised sum gives slightly different results from a sequential one, and different thread counts give different results. For reproducibility (Module 9's deterministic replay), fix the order of reduction.
• Comparing results between two implementations needs tolerances.
• Compilers do not reorder floating-point arithmetic unless told to, which is why -ffast-math changes results.