Variance, covariance & the covariance matrix

The last lesson left a general question hanging: a Markov chain linked the state now to the state next, so how do any two random quantities relate? Here is the answer, and it has a number. You already know expected value — the centre of a distribution; the next two questions are how data spreads around that centre, and how two features move together. Spread and co-movement: that is the whole content of the covariance matrix, and it powers half of classical ML.

Variance: spread of one feature

Variance is the average squared distance from the mean:

Var(X) = E[(X − μ)²]

Squaring keeps everything positive and punishes far-out points. Its square root, the standard deviation σ, is back in the original units and is what you usually report.

Covariance: do two features move together?

Covariance asks whether two features deviate from their means in sync:

Cov(X, Y) = E[(X − μₓ)(Y − μᵧ)]

When X is above its mean and Y is above its mean, the product is positive. Consistently positive → they rise together. Consistently negative → one rises as the other falls. Near zero → no linear relationship.

The catch: covariance’s units are units of X × units of Y — meaningless to compare across pairs. So we normalize.

Correlation: covariance, made comparable

corr(X, Y) = Cov(X, Y) / (σₓ · σᵧ)     ∈ [−1, 1]

+1 is a perfect upward line, −1 a perfect downward line, 0 no linear relationship. Now you can compare relationships across totally different scales.

The covariance matrix

For d features, collect every pairwise covariance into a d × d matrix Σ: variances on the diagonal, covariances off it. It’s symmetric (Cov(X,Y) = Cov(Y,X)) and positive semi-definite. The ellipse you just rotated is Σ — its eigenvectors are the axes, its eigenvalues the spreads along them.

import numpy as np
rng = np.random.default_rng(0)

height = rng.normal(170, 8, 500)
weight = 0.9 * height - 100 + rng.normal(0, 5, 500)   # correlated with height
noise  = rng.normal(0, 10, 500)                        # unrelated

X = np.vstack([height, weight, noise])   # 3 features x 500 samples
print("covariance matrix:\n", np.cov(X).round(1))
print("\ncorrelation matrix:\n", np.corrcoef(X).round(2))
# height-weight correlation is high; everything with 'noise' is ~0

covariance matrix:
 [[ 65.9  59.5   2.1]
 [ 59.5  75.8   0.7]
 [  2.1   0.7 102.2]]

correlation matrix:
 [[1.   0.84 0.03]
 [0.84 1.   0.01]
 [0.03 0.01 1.  ]]

The correlation matrix reads the relationships off at a glance: height and weight sit at 0.84 (tightly linked, as built), while everything paired with noise is ~0.03 — no linear relationship at all. Note the raw covariance between height and weight (59.5) is impossible to interpret on its own; only after dividing by the two standard deviations does it become the clean, comparable 0.84.

Where this lives in ML

• PCA eigendecomposes Σ — the directions of maximum spread are its top eigenvectors (you’ve seen this twice now).
• The multivariate Gaussian is defined by a mean vector and Σ.
• Mahalanobis distance uses Σ⁻¹ to measure distance “in standard deviations,” accounting for correlated, differently-scaled features.
• Multicollinearity is just high off-diagonal correlation — the rank problem from a statistical angle.

In one breath

Variance E[(X−μ)²] measures how one feature spreads (its root, the standard deviation σ, is back in original units). Covariance E[(X−μₓ)(Y−μᵧ)] measures whether two features deviate from their means in sync — positive means they rise together, negative means one rises as the other falls, near-zero means no linear link. Because covariance carries the product of both units, you normalise it to correlation Cov/(σₓσᵧ) ∈ [−1, 1], a unit-free strength of linear association. Stack all pairwise covariances into the symmetric, positive-semidefinite covariance matrix Σ (variances on the diagonal, covariances off it) — the same Σ that PCA eigendecomposes, that defines the multivariate Gaussian, and that Mahalanobis distance inverts. Two traps: correlation isn’t causation, and zero correlation rules out only linear dependence (a perfect parabola y = x² has correlation 0).

Practice

A question to carry forward

Look at what the covariance matrix handed us: an ellipse. Σ said the cloud of points is wide here, narrow there, tilted just so — its eigenvectors the axes, its eigenvalues the spreads. But Σ is only a description of the shape; it is not yet a distribution. It cannot tell you how likely any particular point is, only how the cloud is stretched.

So here is the thread onward: what is the full probability distribution that wears this shape — centred at a mean vector μ, stretched and tilted by exactly this Σ? It is the multivariate Gaussian, the bell curve lifted into many dimensions and the one rich distribution we can actually compute with. What is its density, why does that same elliptical “Mahalanobis distance” sit in the exponent, and what closure magic — marginals, conditionals, and linear maps all staying Gaussian — makes it the backbone of Gaussian processes, Kalman filters, VAEs, and anomaly detection?