Stable Covariance and Cross-Product Computation

Why Centering Is a Numerical Problem

A covariance matrix is built from deviations around a mean. In exact arithmetic, this is straightforward; in floating-point arithmetic, the order of operations can determine whether small variation survives. This matters when measurements share a large common offset, when observations arrive in a stream, or when a computed covariance matrix is passed to a later matrix algorithm.

The central choice is between computing deviations first and accumulating their cross-products, or subtracting two large raw-moment quantities. The first approach is usually safer. An online update offers a way to maintain the same centered quantity as observations arrive, without storing the entire data set.

Definition: Let $x_1,\ldots,x_n$ be column vectors in $\mathbb{R}^p$, with $n\geq 1$, and let $\bar{x}=n^{-1}\sum_{i=1}^n x_i$. The centered cross-product matrix is $$ C_n=\sum_{i=1}^n (x_i-\bar{x})(x_i-\bar{x})^T. $$ For $n>1$, the usual sample covariance matrix is $S_n=C_n/(n-1)$. If $X$ is the $n\times p$ data matrix with rows $x_i^T$, then $C_n=X_c^T X_c$, where $X_c$ has rows $(x_i-\bar{x})^T$.

The cross-product matrix is unscaled; the divisor distinguishes it from a covariance estimate. For two paired vector measurements $x_i\in\mathbb{R}^p$ and $y_i\in\mathbb{R}^q$, the centered cross-product is $C_{xy}=\sum_i(x_i-\bar{x})(y_i-\bar{y})^T$, a $p\times q$ matrix. Its diagonal-block special case $x_i=y_i$ is $C_n$.

Two-Pass Computation

The direct two-pass method first estimates the mean and then forms deviations from it. In exact arithmetic, it is simply the definition evaluated in two stages. In floating-point work, accumulate the mean with a numerically sensible summation strategy, then accumulate the outer products of the deviations. The final result can be symmetrized to remove tiny asymmetry from rounding.

Compute the mean.
Set $\bar{x}=n^{-1}\sum_{i=1}^n x_i$, using pairwise or compensated summation when accumulation error may matter.

Center the observations.
For each observation, calculate $d_i=x_i-\bar{x}$.

Accumulate cross-products.
Compute $C_n=\sum_i d_i d_i^T$, then $S_n=C_n/(n-1)$ if $n>1$.

Each outer product $d_i d_i^T$ is symmetric and positive semidefinite. Their sum has those properties in exact arithmetic as well. The two-pass method does not eliminate all roundoff: the mean and deviations are still computed approximately, and the sum of outer products still incurs error. It does, however, avoid the most severe cancellation in the raw-moment formula.

Worked Example: A Two-Dimensional Covariance Matrix

Take $x_1=(1,2)^T$, $x_2=(2,0)^T$, and $x_3=(3,4)^T$. Compute the centered cross-product and sample covariance.

The mean is $\bar{x}=((1+2+3)/3,(2+0+4)/3)^T=(2,2)^T$. The deviations are $(-1,0)^T$, $(0,-2)^T$, and $(1,2)^T$. Their outer products are \[ \begin{pmatrix}1&0\\0&0\end{pmatrix},\qquad \begin{pmatrix}0&0\\0&4\end{pmatrix},\qquad \begin{pmatrix}1&2\\2&4\end{pmatrix}. \]

C_3= \begin{pmatrix}1&0\\0&0\end{pmatrix} + \begin{pmatrix}0&0\\0&4\end{pmatrix} + \begin{pmatrix}1&2\\2&4\end{pmatrix} = \begin{pmatrix}2&2\\2&8\end{pmatrix}, \qquad S_3=\frac{C_3}{2}= \begin{pmatrix}1&1\\1&4\end{pmatrix}.

As a check, $C_3$ is symmetric, and its determinant is $2(8)-2(2)=12>0$; its diagonal entries are positive. The computation uses the sample divisor $n-1=2$, not $n$.

Online Updates for Centered Cross-Products

When data arrive one observation at a time, recomputing all deviations after each arrival is unnecessary. Maintain the current count $k$, mean $\bar{x}_k$, and centered cross-product $C_k$. For a new observation $x_{k+1}$, let $\delta=x_{k+1}-\bar{x}_k$. Update the mean and cross-product as follows.

Online update: For $k\geq 1$, initialize $\bar{x}_1=x_1$ and $C_1=0$. Given $x_{k+1}$, use $$ \bar{x}_{k+1}=\bar{x}_k+\frac{\delta}{k+1}, \qquad C_{k+1}=C_k+\frac{k}{k+1}\delta\delta^T, \qquad \delta=x_{k+1}-\bar{x}_k. $$

The formula is often called the multivariate Welford update. Its coefficient and the fact that the outer product uses the old-mean difference are important; replacing either casually changes the quantity being accumulated.

Derivation: Write $m=\bar{x}_k$, $x=x_{k+1}$, and $\delta=x-m$. The new mean is $m+\delta/(k+1)$. The old deviations, when recentered at the new mean, each decrease by $\delta/(k+1)$. Since their sum about $m$ is zero, their new outer products sum to $C_k+k\delta\delta^T/(k+1)^2$. The new observation has deviation $k\delta/(k+1)$, contributing $k^2\delta\delta^T/(k+1)^2$. Adding these terms gives $C_k+k\delta\delta^T/(k+1)$, as stated.

For paired vectors $x_{k+1}$ and $y_{k+1}$, maintain their respective means and define $\delta_x=x_{k+1}-\bar{x}_k$, $\delta_y=y_{k+1}-\bar{y}_k$. The centered cross-product update is $C_{xy,k+1}=C_{xy,k}+\frac{k}{k+1}\delta_x\delta_y^T$. The same recentering argument applies, with a rectangular outer product.

Worked Example: Updating the Cross-Product Sequentially

Use the three observations from the previous calculation, in order: $(1,2)^T$, $(2,0)^T$, $(3,4)^T$. Initially $\bar{x}_1=(1,2)^T$ and $C_1=0$.

For the second observation, $k=1$ and $\delta=(1,-2)^T$. Thus $\bar{x}_2=(1,2)^T+\frac12(1,-2)^T=(1.5,1)^T$, and $C_2=\frac12\begin{pmatrix}1&-2\\-2&4\end{pmatrix} =\begin{pmatrix}0.5&-1\\-1&2\end{pmatrix}$.

For the third observation, $k=2$ and $\delta=(3,4)^T-(1.5,1)^T=(1.5,3)^T$. The mean becomes $(1.5,1)^T+\frac13(1.5,3)^T=(2,2)^T$. The increment is $\frac23\begin{pmatrix}2.25&4.5\\4.5&9\end{pmatrix} =\begin{pmatrix}1.5&3\\3&6\end{pmatrix}$. Therefore $C_3=\begin{pmatrix}0.5&-1\\-1&2\end{pmatrix} +\begin{pmatrix}1.5&3\\3&6\end{pmatrix} =\begin{pmatrix}2&2\\2&8\end{pmatrix}$, matching the two-pass result.

Why Raw-Moment Subtraction Can Fail

An algebraically equivalent identity is $C_n=\sum_i x_i x_i^T-n\bar{x}\bar{x}^T$. It can be attractive because it appears to require only sums and sums of squares. But when all observations are near a large vector while their deviations are small, both terms on the right are large and nearly equal. Their subtraction can erase the variation of interest.

Worked Example: Large Offset, Small Variation

Consider the scalar observations $M-1,M,M+1$, where $M=10^8$. Their mean is $M$, so the centered deviations are $-1,0,1$. The centered sum of squares is $1+0+1=2$, and the sample variance is $2/(3-1)=1$.

The raw-moment calculation instead subtracts $\sum_i x_i^2=3M^2+2$ and $3\bar{x}^2=3M^2$. In exact arithmetic their difference is indeed $2$. But the terms are about $3\times10^{16}$, while the desired answer is only $2$. In double precision, rounding errors in quantities of that scale can be comparable to or larger than the residual. The subtraction can therefore return a poor value, including a small negative value, even though the exact sum of squares cannot be negative.

Centering the data before squaring keeps the accumulated terms at the scale of the actual variation. If the offset itself makes mean accumulation problematic, choose a convenient reference $a$ and compute the mean as $a+n^{-1}\sum_i(x_i-a)$; the deviations from the resulting mean are unchanged in exact arithmetic. This shift reduces the size of intermediate sums without changing the target covariance.

Exact Properties and Floating-Point Diagnostics

Two properties are useful checks on a computed centered cross-product. First, $C_n$ is symmetric. Second, it is positive semidefinite: for every $v\in\mathbb{R}^p$,

v^T C_n v =\sum_{i=1}^n \left(v^T(x_i-\bar{x})\right)^2 \geq 0.

This identity proves positive semidefiniteness directly. It also shows why a substantially negative quadratic form or eigenvalue indicates a computation problem or an incorrect input to the calculation: such a result is impossible for the exact centered cross-product. Tiny negative eigenvalues can arise from roundoff, especially when the matrix is nearly singular.

In floating-point arithmetic, the usual model for a basic operation is $\operatorname{fl}(a\mathbin{\circ}b)=(a\mathbin{\circ}b)(1+\varepsilon)$, with $|\varepsilon|\leq u$, when the operation is performed without overflow or underflow; $u$ is the unit roundoff. This model explains why cancellation is dangerous: errors in two large quantities can dominate their small difference. It does not give a universal error bound for every matrix computation, since the accumulation order, dimension, data scale, and conditioning all matter.

Practical diagnostic: Check whether the computed matrix is close to its transpose, for example by comparing $\|C-C^T\|$ with the scale of $\|C\|$. If the discrepancy is consistent with rounding, replace $C$ by $(C+C^T)/2$ before a symmetric-matrix calculation. Also inspect the smallest eigenvalue relative to the matrix scale. Symmetrizing removes antisymmetric roundoff; it does not correct a materially negative eigenvalue or recover variation already lost through cancellation.

For a two-pass implementation, accumulating each matrix entry with pairwise or compensated summation can reduce summation error. In an online implementation, each update avoids subtracting two large raw moments, but repeated updates still accumulate rounding error; periodically comparing with a two-pass recomputation can be useful when observations are retained. In either case, scaling variables to comparable magnitudes can improve numerical behavior, provided the scaling is recorded and reversed correctly when interpreting the covariance.

Common Mistake

A common but incorrect argument is: “The raw-moment identity is algebraically equal to the centered sum of squares, so both methods are equally stable.” Algebraic equality in exact arithmetic does not imply equal behavior after rounding. When the mean is large relative to the spread, the raw-moment formula subtracts nearly equal large quantities; the two-pass and online methods work more directly with deviations. Conversely, a tiny negative eigenvalue is not automatic proof that the data are invalid: it may be roundoff, so assess its size relative to the matrix scale before drawing a conclusion.

Key takeaway: Form centered cross-products from deviations, either in two passes or with the online update. Treat symmetry and positive semidefiniteness as exact mathematical properties and numerical diagnostics, while recognizing that finite precision can perturb them.

Check Your Understanding

Use these questions to check both the formulas and the numerical reasoning.

For $n>1$, how do the centered cross-product matrix and the sample covariance matrix differ?
State the online updates for $\bar{x}_{k+1}$ and $C_{k+1}$, and identify which mean is used to define $\delta$.
Why can $\sum_i x_i^2-n\bar{x}^2$ be inaccurate when the observations have a large common offset?
Show from the definition that $v^T C_n v\geq 0$ for every vector $v$.
What does symmetrizing a computed matrix address, and what numerical problem does it not solve?