Algorithms for calculating variance

Algorithms for calculating variance are computational procedures used to determine the variance—the second central moment and a fundamental measure of statistical dispersion—of a set of numerical data. Although the mathematical definition of variance is elementary, its direct numerical evaluation on digital computers is surprisingly subtle. Floating-point arithmetic introduces rounding errors that can be catastrophically amplified by the most obvious formulations, while practical constraints such as very large datasets, streaming data, weighted observations, and parallel architectures impose further requirements on memory use and pass structure. As a result, statisticians and numerical analysts have developed a family of algorithms—ranging from the simple two-pass method to Welford's online update and Chan–Golub–LeVeque parallel combination formulas—each embodying different trade-offs among numerical accuracy, speed, memory consumption, and applicability to sequential or distributed computation.

Background

For a population of N values x₁, x₂, …, xN with mean μ = (1/N) Σ xᵢ, the population variance is defined as

σ² = (1/N) Σᵢ (xᵢ − μ)²,

while for a sample of n observations the usual unbiased estimator of the population variance is the sample variance

s² = (1/(n − 1)) Σᵢ (xᵢ − x̄)²,

where x̄ is the sample mean. In both cases the variance is an average (or weighted average) of squared deviations from a mean, and the standard deviation is its square root.

Translating this definition into a computer program raises two distinct classes of problems. The first is numerical: in binary floating-point arithmetic, subtraction of nearly equal quantities causes loss of significance (catastrophic cancellation), and naive formulations of the variance are particularly vulnerable because the variance is often small relative to the magnitude of the data itself. The second is architectural: data may arrive as a stream that cannot be stored in its entirety, may be partitioned across the nodes of a distributed system, or may carry non-uniform weights, and an algorithm suited to one setting may be unsuited to another. The literature on variance computation is therefore organized around algorithms that are optimized for accuracy (two-pass methods), for a single traversal of the data (online algorithms), and for parallel combination (merge formulas).

The Two-Pass Algorithm

The most accurate straightforward method is the two-pass algorithm. In the first pass, the mean x̄ is computed; in the second pass, the sum of squared deviations Σ(xᵢ − x̄)² is accumulated and divided by n (or n − 1). Because the quantities subtracted—xᵢ and x̄—are close in magnitude, the subtraction is well conditioned, and the accumulated deviations are relatively small numbers whose sum loses little precision.

Formal error analysis (notably by Chan, Golub, and LeVeque) shows that the two-pass algorithm has forward error that grows only moderately with n, making it the reference standard for accuracy. Its drawbacks are practical: it requires that all data be available twice (either stored or re-read), which may be impossible for streams or for datasets exceeding memory, and it doubles data-transfer costs, which can dominate runtime in modern I/O-bound systems.

A refinement, the corrected two-pass algorithm, compensates for the rounding error committed in computing the mean itself by subtracting a correction term based on the sum of the deviations:

s² = [ Σᵢ (xᵢ − x̄)² − (Σᵢ (xᵢ − x̄))² / n ] / (n − 1).

Since Σ(xᵢ − x̄) should be zero but is not exactly so in floating point, its square estimates the error contribution and improves accuracy at negligible cost. Accuracy can be improved further by accumulating the sums with compensated (Kahan) summation or by pairwise summation.

The Textbook One-Pass Formula

The algebraically equivalent identity

s² = ( Σᵢ xᵢ² − n·x̄² ) / (n − 1),

equivalently Var(X) = E[X²] − (E[X])², permits a one-pass computation: Σxᵢ and Σxᵢ² are accumulated simultaneously, and the variance is formed at the end. This formulation appears in many textbooks, handbooks, and early calculators because it is compact and traverses the data once.

Numerically, however, it is the least reliable of the standard methods. Both Σxᵢ² and nx̄² are approximately equal to n(μ² + σ²); their difference, nσ², is obtained by subtracting two large, nearly equal numbers. When the mean is large compared with the standard deviation—that is, when the coefficient of variation σ/|μ| is small—catastrophic cancellation destroys most or all of the significant digits, and the computed "variance" can even be negative. In single-precision arithmetic, computing the variance of values clustered near 10⁸ with a spread of a few units yields almost entirely noise. Consequently, numerical analysts recommend this formula only when the data are known to be centered near zero or when extended-precision accumulation is used. Its failures are a standard cautionary example in numerical analysis.

Welford's Online Algorithm

Welford's algorithm, published by B. P. Welford in 1962 and popularized by Donald Knuth in The Art of Computer Programming (Volume 2, Seminumerical Algorithms), computes the variance in a single pass with constant memory while retaining nearly the accuracy of the two-pass method. It maintains three running quantities—the count n, the current mean, and M₂, the sum of squared deviations Σ(xᵢ − x̄)²—which are updated for each new observation x as follows:

• n ← n + 1
• Δ ← x − mean
• mean ← mean + Δ / n
• Δ₂ ← x − mean
• M₂ ← M₂ + Δ · Δ₂

After the last observation, the sample variance is M₂ / (n − 1) and the population variance is M₂ / n. The key idea is that the update M₂ ← M₂ + Δ·(x − mean) is mathematically equivalent to accumulating squared deviations from the final mean, because the mean used in the correction shifts consistently with the data. The products Δ·Δ₂ are non-negative, so no cancellation between large terms occurs, and the running mean itself remains well conditioned. Empirically, Welford's algorithm is nearly as accurate as the two-pass method while being suitable for streams, sensors, and any setting in which observations arrive sequentially and cannot be stored.

Weighted Incremental Algorithms

Many applications require observations to carry weights wᵢ—counts in grouped data (frequency weights) or measures of precision or reliability. D. H. D. West's 1979 algorithm extends the incremental approach to this case, maintaining the total weight W, the weighted mean, and S, the weighted sum of squared deviations. For each new pair (x, w):

• W ← W + w
• mean_old ← mean
• mean ← mean + (w / W)(x − mean)
• S ← S + w·(x − mean_old)(x − mean)

The choice of divisor for the variance depends on the interpretation of the weights: S/(W − 1) gives the unbiased estimator for frequency weights, whereas reliability weights require the denominator W − (Σwᵢ²)/W. Weighted versions of the parallel combination formulas (below) have likewise been derived, allowing weighted statistics from independent samples to be merged exactly.

Parallel and Distributed Algorithms

In parallel and distributed settings, partial statistics computed on disjoint subsets of the data must be combined. The combination formula of Tony Chan, Gene Golub, and Randall LeVeque (1979, 1983) merges two partial results—(nₐ, meanₐ, M₂,ₐ) and (n_b, mean_b, M₂,b)—as follows:

• Δ = mean_b − meanₐ
• n = nₐ + n_b
• mean = meanₐ + Δ · n_b / n
• M₂ = M₂,ₐ + M₂,b + Δ² · (nₐ·n_b / n)

These updates are associative, so statistics can be combined in a balanced tree over any number of processors or data chunks, in MapReduce frameworks, or in database online-aggregation engines. Chan, Golub, and LeVeque's 1983 paper Algorithms for Computing the Sample Variance: Analysis and Recommendations analyzed the round-off behavior of all the principal methods and recommended such pairwise (tree) reduction schemes, which combine the parallelism of one-pass methods with accuracy approaching that of the two-pass algorithm. The approach has been adopted in stream-processing systems and statistical libraries that must aggregate counts, means, and variances across partitions.

Extensions: Covariance, Higher Moments, and Forgetting Factors

Welford-style updates extend naturally to multivariate statistics. For paired observations (xᵢ, yᵢ), maintaining running means and a co-moment accumulator C enables the online computation of covariance; each update shifts both means by their respective increments and adds Δx·(y − mean_y) to C, with the sample covariance given by C/(n − 1). Pébay's 2008 technical report systematized these "robust, one-pass parallel" formulas and generalized them to covariances and to arbitrary-order statistical moments (skewness, kurtosis, and beyond), including their weighted and mergeable forms. For data streams in which recent observations should count more than old ones, exponentially weighted moving variance algorithms apply an analogous recurrence with a forgetting factor, and are widely used in financial volatility estimation and process monitoring.

Numerical Stability and Accuracy

The central theme uniting these algorithms is the control of rounding error. Catastrophic cancellation arises when subtracting nearly equal floating-point numbers; the textbook one-pass formula suffers from it structurally, whereas the two-pass and Welford algorithms avoid it by never subtracting large comparable accumulations. Theoretical forward-error bounds confirm this ordering: the two-pass method and the corrected two-pass method are most accurate, Welford's algorithm and pairwise-combination methods are close behind, and the naive one-pass formula can be arbitrarily inaccurate as the ratio |μ|/σ grows. Practical remedial strategies include shifting all data by a known approximate mean (which changes the variance not at all but dramatically improves conditioning), accumulating sums in double or extended precision, and using compensated summation. Standard deviations computed by squaring the result inherit the same conditioning issues, and modern libraries expose options for choosing among these strategies.

History

Awareness that the obvious formula for the variance is numerically treacherous developed alongside the spread of digital computers in the 1950s and 1960s. Welford's brief 1962 note in Technometrics, "Note on a Method for Calculating Corrected Sums of Squares and Products," introduced the stable single-pass update later codified by Knuth as Algorithm for computing running mean and variance. West's 1979 note in Communications of the ACM extended the incremental scheme to weighted data. The systematic numerical analysis of all principal methods was accomplished by Chan, Golub, and LeVeque in the late 1970s and early 1983, providing both the merge formulas for parallel computation and the definitive recommendations among competing algorithms. Subsequent work by Pébay and others generalized the mergeable formulations to covariance matrices and arbitrary moments, supplying the theoretical foundation for modern distributed statistical aggregation.

Applications and Significance

Variance-calculation algorithms are embedded throughout scientific and statistical software. Numerical libraries such as NumPy, the GNU Scientific Library, R, and Boost.Accumulators implement two-pass, incremental, or combination variants chosen according to input format; stream-processing and big-data platforms use Chan-style merges to aggregate variance statistics across clusters; online aggregates in database systems and telemetry pipelines rely on Welford's constant-memory update; and machine-learning pipelines use running mean and variance estimators in data normalization schemes such as batch normalization. In finance and quality control, streaming and exponentially weighted variance estimators underpin volatility measurement and control-chart monitoring.

Beyond its direct utility, the problem of computing the variance has had significant pedagogical and methodological influence. It is a canonical case study in numerical analysis—illustrating how algebraically equivalent formulas can differ enormously in stability—and a standard demonstration that algorithm design must account jointly for mathematics, floating-point behavior, and computing architecture. The evolution from the naive one-pass formula to Welford's online update and to parallel combination formulas exemplifies how a seemingly trivial statistical quantity can motivate a rich and practically important body of computational technique.