Negative variance: why a correct formula can fail on a computer

When an impossible number reveals a real problem

An AI pipeline measures variation before normalising data or checking a sensor. Values are large, but differences between them are small. The program returns negative variance. Since variance averages squares, that sign cannot describe the intended mathematical quantity. Should we investigate the data, the formula, or the execution of arithmetic? We isolate the third possibility using three synthetic numbers.

The question is precise: are expressions equivalent over real numbers interchangeable at finite precision? We show why not, compare a centred two-pass calculation with an incremental update, and distinguish arithmetic error from information already lost in the input. We measure neither a production line nor an EL-AI model’s accuracy. Units are generic: volt inputs would produce variance in squared volts, but no data comes from a real voltmeter.

Define deviations before choosing the algorithm

The arithmetic mean is the balance point of the values. A deviation is a value minus that mean. Variance describes the mean squared deviation; squaring prevents positive and negative differences from cancelling. Let n denote the number of observations, x_i the ith value, and mu their mean. We divide by n, using the population variance of the finite population comprising our three values.

Expand each square to obtain the last equality. Summing −2 mu x_i gives −2 n mu², while mu² appears n times. Dividing by n leaves the mean of squares minus the square of the mean. Exact arithmetic yields the same nonnegative result. The final expression appears convenient because it accumulates only sums and squared sums. Algebraic convenience, however, does not guarantee numerical stability.

The n−1 correction answers a different statistical question: estimating a broader population’s variance from a sample under suitable assumptions. It does not repair numerical cancellation. For these three values, changing the denominator multiplies variance by 3/2 without recovering lost digits. NumPy calls explicitly use ddof=0 so that all methods target the same mathematical quantity.

The worked example, even without a computer

Take b=1,000,000,000 and values b+1, b+2, b+4. Their mean is b+7/3. Deviations are −4/3, −1/3, and 5/3; their squares sum to 42/9, and dividing by three gives 14/9. The large baseline disappears because a common translation leaves distances from the mean unchanged. Exact variance is approximately 1.5555555556, neither zero nor negative.

Executed with Python 3.14.0, 53-bit float significands, and sequential accumulation, the raw formula instead returns −128. The deterministic search uses baselines 10^k for k=0 through 16 with fixed increments 1, 2, 4. The first negative result in that order occurs at k=9. Stating the search rule prevents presenting a selected example as an error frequency. We do not claim every triple or computer produces −128.

Where the digits disappear

A floating-point number retains a finite number of significant digits, not a fixed number of decimal places. As magnitude increases, representable numbers become farther apart. Near b, math.ulp(b) is about 1.19×10⁻⁷, so the integer increments remain clearly distinguishable. At k=9 the problem is therefore not that 1, 2, and 4 disappeared during input. It arises in later operations, when large values are squared.

The raw formula subtracts two quantities near 10¹⁸ to obtain a residual near 1.56. Rounding in products, sums, and the mean can be large relative to that residual. Subtraction exposes the uncertainty: it cancels common leading digits without recreating unavailable trailing digits. It is not correct to say every subtraction of nearby numbers is intrinsically inaccurate; existing operand errors dominate the small result we seek.

Replacing accumulation with math.fsum is useful but not a universal cure. fsum improves the sum of its inputs; it does not make already rounded x*x products or the squared mean exact. In the worked example, the fsum variant also returns −128. Clamping with max(v,0) removes the impossible sign by returning zero, while the correct variance remains 14/9. Clamping may be a containment policy in a wider system, never evidence of recovered information.

Centre before squaring

The two-pass method first estimates the mean, then sums squared deviations from it. The squared quantities are small; no final subtraction of enormous terms is required. Our code uses fsum for both the mean and deviations. If the computed centre is mu+e, exact arithmetic gives mean squared distance v+e² because true deviations sum to zero. A small centre error therefore does not automatically become a large variance error. Other arithmetic operations still round, however.

For a stream we do not wish to reread, the update commonly called Welford maintains a count, current mean, and M2, the centred sum of squares. When x_k arrives, delta is its distance from the old mean; the mean shifts by delta/k. M2 is updated using delta times the distance from the new mean. Final variance is M2/n, not M2/(n−1) for our chosen comparison.

The recurrence can be checked without invoking a library. Shifting the k−1 old values’ centre by delta/k increases their contribution by (k−1)delta²/k², since their old deviations sum to zero. The new value lies delta(k−1)/k from the new mean, adding delta²(k−1)²/k². Together these give delta²(k−1)/k, equal to delta(x_k−mu_k). This derivation is exact over reals; the implementation does not thereby inherit infinite precision.

Results and the meaning of an exact reference

The two-pass method and NumPy agree in this case; Welford is very close but not identical to the correct value. This is not a universal library contest. NumPy receives float64 with ddof=0; the other methods are short implementations in our experiment. Fraction.from_float converts each actual floating-point input to an exact rational, then computes variance without rounding. At k=9 the integer inputs are exactly represented, so this reference also matches the intended generated values.

At larger baselines an important distinction appears. Near 10¹⁶ adjacent floats are separated by 2, so an increment of 1 is no longer always preserved as intended. A variance algorithm cannot recover it after conversion. The plot therefore measures error against the exact variance of the actually represented inputs, rather than assuming 14/9 throughout the sweep. An inappropriate oracle would conflate input representation with arithmetic stability.

The program also checks all six permutations of the worked triple. In this environment the raw formula stays at −128, whereas Welford ranges approximately from 1.5555555224 to 1.5555555820. Even more stable algorithms can depend on ordering at finite precision. The baseline sweep and permutations are deterministic checks, not a formal stability proof for every sequence. Complete results are in the JSON, rather than only favourable cases.

Choose a method without confusing the trade-offs

Both two-pass and Welford methods require work proportional to n. Welford retains constant-size state; two passes require rereading or storing data. This describes the algorithms, not our script’s total memory, which retains inputs and results for reproducibility. Fraction supplies an exact reference here, but arithmetic on growing integers does not have constant float cost. No timing or memory benchmark was performed.

Centring does not solve every scale: non-finite inputs, product overflow, underflow of tiny deviations, or means outside representable range need further handling. Information lost before calculation may require a different representation or acquiring differences from a smaller reference. If statistics normalise model inputs, small variance often enters a denominator, so near-zero handling must be explicit too. These three numbers establish no predictive-quality improvement.

The formula is an identity; the algorithm is a choice

The answer is no: algebraic equivalence does not imply numerical equivalence. In the worked case, exact variance 14/9 survives in the centred calculation and nearly entirely in the incremental one, but is overwhelmed by the raw formula’s subtraction of large terms. The negative sign is a useful symptom, not a value to correct blindly. Before blaming an AI model, identify the intended quantity, the represented inputs, and the executed arithmetic sequence.

Sources, code, and scope

We read numpy.var technical notes on its formula, ddof, and precision, plus Python definitions of fsum and ulp. The derivations and synthetic triple are this article’s teaching analysis, not a reproduction of a published benchmark. NumPy 2.5 and Python 3.14 documentation were accessed on 21 October 2026; the executed environment was Python 3.14.0 and NumPy 2.5.3, recorded in the results. No sensor, GPU, or embedded-platform execution and no claimed human review.

NumPy 2.5 — numpy.var.

Python 3.14 — math.fsum, math.ulp.

from fractions import Fraction
xs=[1000000001.0,1000000002.0,1000000004.0]
total=squares=0.0
for x in xs:
    total+=x
    squares+=x*x
raw=squares/len(xs)-(total/len(xs))**2
ys=[Fraction.from_float(x) for x in xs]
mean=sum(ys)/len(ys)
reference=sum((x-mean)**2 for x in ys)/len(ys)
assert raw==-128.0 and reference==Fraction(14,9)
print(raw,reference)

Code, data, and instructions · JSON. Educational calculations executed with Python 3.14.0; figures with Matplotlib 3.11.2. AI-assisted analysis, without claiming peer review or human review. Original illustrative ImageGen cover: it does not document EL-AI people, premises, or installations. Sources accessed 21 October 2026.