The formula sets the power, the sum sets the constant

Cancellation takes the answer, not a digit put a rule beside the subtraction it was about. If a computation forms a small quantity $s$ as the difference of operands of size $m$, the relative error in $s$ is about $u\,m/|s|$. For a single subtraction of correctly rounded operands the rule held, and was conservative by between two and a half and twenty-eight times. For the one-pass variance — the mean of the squares minus the square of the mean — it failed in the dangerous direction, under-predicting by 11, 11, 69 and 46 times on four spreads of a thousand values near $10^6$. The explanation given was that the operands of that subtraction were not stored numbers but the output of a thousand-term accumulation, so $u\,m$ was the wrong estimate of the error they carried.

That is true, and it leaves the factor unpriced. Eleven and sixty-nine are not the same number, the essay did not say where between them the next dataset would land, and nothing said which part of the failure belongs to the formula and which to the sums inside it.

The rule’s ratio is a condition number #

For the one-pass variance the operands are $\frac1n\sum x_i^2$ and $\bar x^{\,2}$, and the difference is the variance $v$. So the rule’s $m/|s|$ is the mean square over the variance:

$$\frac{m}{|s|} \;=\; \frac{\bar x^{\,2} + v}{v} \;=\; 1 + \frac{\bar x^{\,2}}{v} \;=\; \kappa^2 .$$

$\kappa$ is the condition number of the variance with respect to the data: how far a relative change in the values can move it. A sample whose mean is a thousand standard deviations from zero has $\kappa = 1000$, and its variance is determined only to about $\kappa\,u$ by the data whatever is done with it — but every formula below is applied to data stored exactly, so the question is purely what the arithmetic adds. The rule says the one-pass formula adds $\kappa^2 u$.

To measure that without a reference that carries its own error, every value is taken as the exact rational it is — a double is a 53-bit integer times a power of two — and $\sum x_i$ and $\sum x_i^2$ are formed exactly in integer arithmetic before one division to sixty bits — the move an answer that is known makes for a linear system, applied to two sums. The earlier essay compared against the two-pass formula, which is accurate to a few units of roundoff and, measured here, never more than 43 off on its data; that was good enough for its table, and it is not good enough to separate a constant of one from a constant of three.

Six formulas against κ #

Six ways of computing the same number, each in the working precision at every operation: the one-pass formula with running sums, as written in every textbook and many instruction sets; the same formula with both sums formed pairwise; the same with both sums compensated; the one-pass formula applied to every value minus the first; Welford’s running-mean update; and two passes, the mean first and then the squares of the deviations. Each is applied to a thousand Gaussian values of standard deviation one and mean $\sqrt{\kappa^2-1}$, eleven seeds, the median error drawn.

The picture is four slopes, not six curves. Everything built on the one-pass formula climbs with slope one in $\kappa^2$ — the square of the condition number — whatever is done to the sums inside it. At $\kappa^2 = 3.2\cdot 10^{14}$ the running-sum version is $7.4\,\kappa^2 u$ from the answer and the compensated one exactly $1.0\,\kappa^2 u$. Welford’s update climbs with slope one half: its error divided by $\kappa u$ sits between 0.06 and 1.2 across the range. And the two formulas that subtract the mean before squaring do not climb at all. Two passes end at $4\,u$ at the top of the range, where the one-pass formula is wrong by a quarter of the answer.

So the rule was right about the exponent. What it called “the operands’ size over the result’s” is a property of the formula: any formula that squares values first and subtracts the squares afterwards forms a difference whose operands are $\kappa^2$ times the result, and spends $\kappa^2$. Any formula that subtracts first squares numbers the size of the answer and spends nothing.

The dial does not change that. Above $\kappa^2 = 10$, at 16, 24 and 32 bits, the one-pass constant is 3 to 26 times $\kappa^2 u$, the compensated one 0.3 to 1.3, Welford’s 0.07 to 1.7 times $\kappa u$, and two passes 3 to 8 $u$ — the same bands double precision gives. What the precision moves is where the one-pass formula runs out — at roughly $\kappa^2 = 1/(10u)$, which is $1.7\cdot 10^6$ in single precision — and where the stored data stops having a spread worth measuring at all.

The constant in front belongs to the sum #

That leaves the constant: seven at the top of one sweep, eleven and sixty-nine in the earlier essay’s table. Hold $\kappa^2$ at $10^8$ and change the number of values.

The one-pass constant is $0.4\sqrt n$ to within the scatter of fifteen seeds: 2.1 at a hundred values, 14.5 at a thousand, 38 at ten thousand, 111 at a hundred thousand. That is the running sum’s own random walk. Each addition in $\sum x_i^2$ rounds a partial sum near $k\,\bar x^{\,2}$, the errors take both signs, and their total grows like the square root of the count — the same half-power three walks and one bound found under a left-to-right sum, a chain of rotations and a residual recurrence. The one-pass formula then multiplies that walk by $\kappa^2$, because the subtraction exposes it.

Sum pairwise and the walk is gone: a tree of depth $\log_2 n$ has nothing long enough to wander, and the constant is 0.5 to 1.7 at every length from a hundred to a hundred thousand. Compensated summation does the same thing by carrying each rounding forward. Neither does anything to $\kappa^2$.

Two passes grow like $\sqrt n$ as well — 1.6 $u$ at a hundred values, 73 at a hundred thousand — and that is not a contradiction. Their second pass sums the squared deviations with a running sum, and that walk is there to be seen. It is simply not multiplied by anything, because the quantities being summed are already the size of the answer.

So the factor the earlier essay left unexplained has two parts, and they come from different lines of the program. The $\kappa^2$ is the formula’s: it squares before it subtracts. The $0.4\sqrt n$ is the summation’s: it walks. A thousand values give $0.4\sqrt{1000} = 12.6$, against the 11.1 and 11.4 the earlier table found at its two narrowest spreads.

A better sum restores the rule and no more #

There is a piece of advice attached to the one-pass variance almost as often as the warning against it: if the formula must be used, sum it accurately. The measurement says exactly what that buys.

It buys the $\sqrt n$. With compensated sums the error is $0.5$ to $1.2\,\kappa^2 u$ at every length and $0.3$ to $1.3\,\kappa^2u$ at every precision — the digits-lost rule, met to its constant. It does not buy anything below the rule, and there is a reason that has nothing to do with how well the sums are formed. Even an exactly computed $\sum x_i^2$ has to be rounded once to be stored, and so does $(\sum x_i)^2/n$; each of those roundings is $u$ relative to a quantity $\kappa^2$ times the variance. Two numbers each wrong by $u\kappa^2 v$ cannot be subtracted to anything better. At $\kappa^2 = 10^{12}$ the compensated formula is still $4.7\cdot 10^{11}\,u$ from the exact variance, where two passes over the same data are six.

A summation that is perfect inside a formula that squares first is the rule, exactly. That is an improvement of $0.4\sqrt n$, which for a hundred thousand values is a factor of 130, and it is worth having when a single pass is forced. It is not a repair.

The sixty-nine was the data #

Thirty-six times $\kappa^2 u$ is the largest constant forty-one seeds of a well-mixed random generator give at a thousand values near $10^6$ with a spread of one. The earlier essay’s 69 at a spread of one, and 46 at a spread of ten, sit outside that.

Rerunning its exact recipe — its 31-bit linear congruential generator, values $10^6 + (r - \tfrac12)$, a thousand of them — against the exact variance reproduces its numbers to the digit: 11.1, 11.4, 68.9 and 45.5. Then running the same recipe from forty-one seeds, and the same mean and spread from the well-mixed generator:

The earlier essay’s data does not scatter around $0.4\sqrt{n}$. Its constant is 48 to 87 with a median of 62, and on all forty-one seeds the formula over-estimates. The well-mixed generator, asked for the same mean and the same spread, gives 0.7 to 36 with a median of 12, and over-estimates on sixteen seeds and under-estimates on twenty-five — the signs a random walk should have. With compensated sums the two generators are indistinguishable: both below one.

A constant five times the walk and one sign on forty-one seeds is not a walk. Splitting the error by sum says where it lives.

The lean is in $\sum x_i$. On the congruential data the running sum of the values falls short by 24 to 39 units of roundoff on every one of nine seeds, while on the well-mixed data it scatters from $-11$ to $+12$. $\sum x_i^2$ looks alike on both. Since $\sum x_i$ enters the formula squared and subtracted, a sum that is always short makes a variance that is always long, by about twice the shortfall times $\kappa^2$ — which is the 62.

That is the phenomenon the direction the error leans measured with the rounding mode as the cause. Under round-to-nearest the errors of ten thousand operations combined with an exponent of 0.47; under rounding toward infinity, 1.01, because every error had the same sign. Here the mode is round-to-nearest throughout and the errors lean anyway. The data, not the arithmetic, is supplying the bias. A coin flip that fixes the average is the same fact from the other side: errors that share a sign accumulate in proportion to their number, and rounding at random is what gave them back both signs there.

What about the data does it is not established here. The obvious candidate fails: the generator’s values sit on a lattice of $2^{-31}$, two bits coarser than the double grid at $10^6$ that what a float can hold draws, but the well-mixed values quantised to the same lattice give a median of 12.9 and no lean. Shuffling the congruential values with an unrelated generator leaves the median at 56, so the lean is in the multiset of values and not in their order. The low-order bits of a power-of-two congruential generator are known to be far from random — the lowest bit alternates — and the bits a running sum near $10^9$ discards are exactly those; that is the mechanism a reader would reach for, and it has not been tested.

The useful fact does not depend on it. A one-pass variance’s constant is $0.4\sqrt n$ only when the rounding errors of its sums are free to cancel, and nothing guarantees that for data with structure in its low bits. The earlier essay happened to draw its examples from such a source. A compensated sum does not care: it removes the walk and the lean together.

A shift is a smaller κ #

Two of the six formulas never climbed. The shifted one-pass formula is the instructive one, because it is the one-pass formula: it squares first and subtracts afterwards, exactly as the failing version does. It is applied to $x_i - x_1$ instead of $x_i$.

That changes $\kappa$, not the formula. The variance of $x_i - K$ is the variance of $x_i$, and the shifted data’s mean is $\bar x - K$, so its condition number is $1 + (\bar x - K)^2/v$. Shift by the first value and $\bar x - x_1$ is about one standard deviation, so $\kappa^2$ is about two whatever the original was.

Shift by a number $t$ standard deviations from the mean, with the data’s own $\kappa^2$ held at $10^8$, and the error lands on the unshifted curve at $\kappa^2 = 1 + t^2$. Past $t = 100$ it is 10 to 18 times $(1+t^2)\,u$, the same constant the unshifted formula carries at a thousand values, and it climbs with slope two in $t$ to within 0.3 at every step. A shift at the mean costs 6 $u$, one standard deviation away 14. A shift a million standard deviations off is no shift at all.

So the cheapest repair of the one-pass formula is not a better sum and not a second pass. It is any value near the data — the first one, a running guess, last month’s mean — subtracted before squaring. The formula is unchanged and still spends $\kappa^2$; the shift chooses which $\kappa$ it spends.

Welford spends one power #

Welford’s update keeps a running mean $m_k$ and a running sum of squared deviations, and adds $(x_k - m_{k-1})(x_k - m_k)$ at each step. It never forms $\sum x_i^2$. It still climbs, with slope one half in $\kappa^2$ — slope one in $\kappa$ — and its constant is 0.15 to 0.43 at every length from a hundred to a hundred thousand.

The mechanism is one rounding in the wrong place. The running mean is a number of size $\bar x$, so it carries an absolute error of about $u\,\bar x$, which is $u\kappa$ standard deviations. Every deviation $x_k - m_k$ is formed against that mean and inherits a relative error of $u\kappa$, and so does each product added to the sum. The errors do not compound, because each step starts from the values rather than from the previous deviation — which is why the constant does not grow with $n$ — but they do not shrink either.

That makes Welford the formula of choice where a single pass is forced and the data cannot be shifted: it costs a division and two subtractions per value more than the one-pass formula and spends a single power of $\kappa$ where the one-pass formula spends two. At $\kappa = 10^4$, a mean ten thousand standard deviations from zero, the one-pass formula is $1.2\cdot 10^9\,u$ from the answer and Welford’s $1.9\cdot 10^3\,u$: seven correct digits against twelve.

What each formula costs, and what it spends #

Read down the third column and the formulas sort themselves; read down the fourth and the sums do. The two columns answer different questions and the earlier essay’s rule answered only the first. The order they are added in is about the fourth column; the condition number is an amplifier is about what the third column multiplies.

The same split is the one the vector that hides it found in a parallel inner product. The summation condition number of the data decides how much of the reduction’s walk reaches the answer; the reduction decides how long the walk is. Positive values sat at the safe end there, with condition number one. The one-pass variance is the opposite: its two sums are of positive numbers and well conditioned on their own, and it is the subtraction afterwards that hands their error a factor of $\kappa^2$.

What a thousand Gaussian values do not show #

Every sweep here uses Gaussian data of unit variance, and every formula is evaluated with each operation rounded once to the working precision, the way the formats are specified. Data with heavy tails could make the running sum’s walk lopsided — a single large value dominates its partial sums — and the constant would no longer be a clean $0.4\sqrt n$. Hardware that accumulates in a wider register than it stores, as some vector units do, would be expected to shrink the one-pass constant toward the compensated one without changing the power; that is not measured here. And the congruential generator’s lean is measured on one mean and one spread; how it varies with either is untested, which is the first question below.

Still open: the low bits, a pairwise Welford, and a mean that moves #

What the congruential data supplies. The lean sits in $\sum x_i$, survives shuffling, and is absent from the well-mixed generator on the same lattice. The prediction with a sign is that replacing only the lowest five bits of each congruential value with bits from an unrelated generator removes the lean — the median constant falls from 62 to within a factor of 1.5 of the well-mixed 12, and the over-estimates fall from forty-one seeds to between ten and thirty — and that replacing only the top bits of the fraction leaves it. If the lean survives the low-bit replacement, it is in the distribution of the values rather than their bit patterns, and the explanation the low bits invite is wrong.

Welford combined pairwise. Welford’s update merges two partial results with a formula of Chan, Golub and LeVeque, and merging them as a tree is how a parallel reduction would compute a variance. The prediction is that the tree keeps Welford’s single power of $\kappa$ and its constant stays under 0.5 — that the update’s error is set by the running mean’s rounding and not by the order of combination — while a tree of one-pass partials keeps two powers.

A mean that drifts. The shift by $x_1$ works because $x_1$ is within a standard deviation or two of the mean. A stream whose mean drifts by $d$ standard deviations over its length makes that shift stale. The prediction is that the shifted formula’s error then follows $1 + d^2/3$ in place of $\kappa^2$ — the mean square distance of a linear drift from its start — and that re-shifting at every block of a hundred values, merging the blocks exactly, removes the drift’s cost for under one per cent more arithmetic.