Covariance, Correlation and the Multivariate Gaussian Explained

Do two things move together? Students who study more tend to score higher. Taller people tend to weigh more. In machine learning, every pair of features either moves together or it does not, and that single fact shapes almost every model you will build. Today we learn to measure it exactly.

• Hours studied vs exam score
• Can we measure togetherness?

Last time, in lesson 37, we studied convexity and Lagrange multipliers. Now we turn to probability. We will build covariance by hand, normalize it into correlation, assemble the covariance matrix, and finally meet the multivariate Gaussian, whose shape is an ellipse. Along the way we will see exactly where correlation lies to you.

• Covariance and correlation
• Covariance matrix
• Multivariate Gaussian

This is not a side topic. Principal component analysis diagonalizes a covariance matrix. Gaussian mixture models fit one ellipse per cluster. Gaussian processes are built from a covariance function. If covariance feels shaky, all three feel like magic. If it feels solid, all three feel like natural consequences.

• PCA uses it
• GMMs use it
• Gaussian processes use it

Covariance starts with deviations. A deviation is how far a value sits from its own mean. Take hours studied: one, two, three, four, five. The mean is three, so the deviations are minus two, minus one, zero, plus one, plus two. Positive means above average, negative means below.

• Subtract the mean
• Deviation = value minus mean

Now pair each student with a score: two, four, five, four, five. The mean score is four. The score deviations are minus two, zero, plus one, zero, plus one. Notice that student one is below average on both measures, and student five is above average on both. That shared direction is what we want to capture.

• x: 1, 2, 3, 4, 5
• y: 2, 4, 5, 4, 5
• Mean of y is 4

Multiply each pair of deviations. Minus two times minus two is four. Zero times anything is zero, and so on, ending with plus one times plus one, plus two times plus one is two. The products are four, zero, zero, zero, and two. Their sum is six. Positive, because the points mostly lie where both deviations share a sign.

• Products: 4, 0, 0, 0, 2
• Sum = 6

Sample covariance divides that sum by n minus one, here four, giving one point five. Why n minus one? Because we estimated the means from the same data, which uses up one degree of freedom, and dividing by n would systematically underestimate. The formula is: covariance equals one over n minus one, times the sum of x deviation times y deviation.

• Cov = sum / (n − 1)
• 6 / 4 = 1.5

The sign tells the story: positive means they rise together, negative means one rises as the other falls, and near zero means no linear trend. But the size is slippery. Measure study time in minutes instead of hours and the covariance jumps by a factor of sixty. A number that changes with units is hard to interpret. We need to normalize.

• Positive: move together
• Negative: move opposite
• Size depends on units

To remove units, divide by both standard deviations. This is the Pearson correlation coefficient, called r. For our data, the variance of x is ten over four, which is two point five, so its standard deviation is about one point five eight. The variance of y is six over four, one point five, so its standard deviation is about one point two two five.

• Divide by std of x and y
• r = Cov / (sx · sy)

Now plug in. One point five divided by one point five eight one times one point two two five is one point five over one point nine three six, about zero point seven seven five. Check with the shortcut: six divided by the square root of ten times six, which is six over seven point seven four six. Same answer. A correlation of about zero point seven seven is a fairly strong positive link.

• 1.5 / (1.581 × 1.225)
• r ≈ 0.775

Correlation always lies between minus one and plus one. Plus one means the points fall exactly on an upward line, minus one on a downward line, and zero means no linear relationship. And because the units cancel, hours versus minutes gives the identical r. Notice the word linear, because it will matter a lot in a few minutes.

• r = 1: perfect line up
• r = −1: perfect line down
• Units cancel

Here is the geometric reading, using lesson 3. Treat the centered x values as a vector, and the centered y values as another. Their dot product is six, their lengths are root ten and root six, and the cosine of the angle between them is six over root sixty, which is zero point seven seven five. So the angle is about thirty nine degrees. Correlation is a cosine.

• Centered vectors in n dimensions
• r = cos θ

That picture explains the bounds. A cosine can never leave the range minus one to one. Vectors pointing the same way give one, perpendicular vectors give zero, opposite vectors give minus one. So uncorrelated literally means the centered data vectors are orthogonal. Centering is essential here: forget it, and you are just computing cosine similarity of raw numbers.

• Angle 0°: r = 1
• Angle 90°: r = 0
• Angle 180°: r = −1

Real data has many features. Put it in a matrix X with n rows, one per observation, and d columns, one per feature. First, center it: subtract each column's mean from that column, giving a centered matrix. After this step, every column has mean zero, and the data cloud is centered at the origin.

• Data matrix X: n rows, d columns
• Center each column

The covariance matrix is capital sigma, equal to X centered transpose, times X centered, all divided by n minus one. It is d by d. Entry in row i, column j is the dot product of centered feature i with centered feature j, divided by n minus one. That is exactly the covariance of those two features, all computed at once by one matrix multiplication.

• Σ = XcᵀXc / (n − 1)
• Shape d × d

For our study data, the matrix is two point five, one point five on top, one point five, one point five underneath. The diagonal holds the variances: two point five for hours, one point five for score. The off-diagonals both hold the covariance, one point five. It is symmetric, because covariance of x with y equals covariance of y with x.

• [[2.5, 1.5],[1.5, 1.5]]
• Diagonal = variances

Now a key property. Pick any direction v. The quantity v transpose sigma v equals the variance of the data projected onto v. A variance can never be negative, so this quadratic form is always at least zero. That is exactly what positive semi-definite means, from lesson 5. It is zero only if some direction carries no spread at all.

• vᵀΣv = variance of projection
• Never negative

Check it on our matrix. The trace is four and the determinant is two point five times one point five minus one point five squared, which is one point five. The eigenvalues solve lambda squared minus four lambda plus one point five, giving about three point five eight and zero point four two. Both are positive, as promised. These are the variances along the two natural axes of the data.

• λ ≈ 3.58 and 0.42
• Trace 4, determinant 1.5
• Both positive

In NumPy, call np dot cov on your data with rowvar equal to False, so that columns are features. You get our two point five, one point five matrix. A common mistake is forgetting that argument, which silently treats rows as variables and returns a wrong-shaped matrix. And np dot corrcoef gives the correlation matrix, with ones on the diagonal and zero point seven seven five off it.

• np.cov(X, rowvar=False)
• np.corrcoef(X, rowvar=False)

Pearson only sees straight lines. Take x equal to one through five and y equal to x cubed: one, eight, twenty seven, sixty four, one hundred twenty five. The relationship is perfectly increasing, yet Pearson gives only about zero point nine four, because the points curve away from any line. A single outlier can distort r even more.

• Pearson measures linear
• Outliers can dominate

Spearman correlation fixes this. Replace each value by its rank, smallest is one, next is two, and so on, then compute Pearson on the ranks. Here the ranks of x and y are both one through five, so the ranked points lie on a perfect line and Spearman equals one. It measures how monotonic a relationship is, and it shrugs off outliers.

• Replace values by ranks
• Pearson on the ranks

So which to use? Pearson when you care about linear structure and the data is roughly well behaved. Spearman when the relationship might be curved but consistently increasing, when there are outliers, or when the data is ordinal. Looking at both is a cheap diagnostic. If they disagree strongly, plot the data before trusting either number.

• Pearson: linear, sensitive to outliers
• Spearman: monotonic, robust

Now the famous trap. Let x be minus two, minus one, zero, one, two, and let y equal x squared: four, one, zero, one, four. The mean of y is two, so the y deviations are two, minus one, minus two, minus one, two. Multiply with the x deviations: minus four, one, zero, minus one, four. They sum to zero.

• x: −2, −1, 0, 1, 2
• y = x²

Covariance is zero, so correlation is zero, yet y is completely determined by x. These variables are as dependent as it gets. The lesson: independence implies zero correlation, but never the other way around. Correlation only detects linear association. In machine learning, this is why a feature with zero correlation with the target may still be hugely useful to a nonlinear model.

• Cov = 0, r = 0
• Knowing x fixes y exactly

There is one important exception. If two variables are jointly Gaussian, meaning together they follow a multivariate normal distribution, then zero correlation does imply independence. This is a special gift of the Gaussian, not a general law. Remember it, because it is one reason Gaussians are so convenient, and also why people wrongly generalize the rule.

• Jointly Gaussian: uncorrelated ⇒ independent
• Special to the Gaussian

Even a large correlation does not mean one variable causes the other. A hidden third factor can drive both, like hot weather raising both ice cream sales and swimming accidents. Or the causation may run backwards. Or it is pure coincidence in a small sample. Correlation is evidence worth investigating, but experiments or careful causal reasoning are needed to claim a cause.

• Hidden confounders
• Reverse causation
• Coincidence

Simpson's paradox shows how a hidden grouping can flip a trend. Here is a made-up example. Among easy cases, drug A cured nine of ten, ninety percent, and drug B cured eighty of one hundred, eighty percent. Among hard cases, A cured thirty of one hundred, thirty percent, and B cured two of ten, twenty percent. Drug A wins in both groups.

• Easy cases: A 9/10, B 80/100
• Hard cases: A 30/100, B 2/10

Now combine the groups. Drug A cured thirty nine out of one hundred ten, about thirty five percent. Drug B cured eighty two out of one hundred ten, about seventy five percent. The overall comparison reverses. Why? Drug B was mostly given to easy cases, and drug A mostly to hard ones. The aggregate correlation is confounded by case difficulty, so the right comparison conditions on the group.

• A overall: 39/110 ≈ 35%
• B overall: 82/110 ≈ 75%

Now let us turn the covariance matrix into a distribution. In one dimension the Gaussian is a bell curve with mean and variance. In two dimensions it becomes a smooth hill, and a top-down view shows nested elliptical contours. The mean mu says where the hill is centered, and the covariance matrix sigma says how it is stretched and tilted.

• 1D: bell curve
• 2D: hill with elliptical contours

The density is one over the square root of two pi to the d, times the determinant of sigma, multiplied by e to the power of minus one half, times x minus mu, transposed, times sigma inverse, times x minus mu. The exponent decides the shape. The front factor is just a constant that makes the total probability equal one, and the determinant accounts for how stretched the hill is.

• p(x) = exp(−½ (x−μ)ᵀΣ⁻¹(x−μ)) / √((2π)ᵈ|Σ|)
• d = number of dimensions

Contours of constant density are ellipses. Take unit variances and covariance zero point eight. The eigenvalues are one point eight and zero point two, with eigenvectors along the diagonal directions. So the ellipse points along one, one, with one standard deviation half length about one point three four, and it is only about zero point four five wide the other way. More correlation means a thinner, more tilted ellipse.

• Axes = eigenvectors of Σ
• Lengths ∝ √eigenvalues

Euclidean distance treats every direction equally, but data does not. Take a distribution with variance four in x and one in y, no correlation. The points two, zero and zero, two are both Euclidean distance two from the center. Yet getting two units away along x is routine, since the standard deviation is two, while two units along y is two standard deviations out.

• Euclidean ignores spread
• Mahalanobis respects covariance

Mahalanobis distance fixes this. Squared distance is x minus mu, transposed, times sigma inverse, times x minus mu. With variances four and one, the point two, zero gives four over four, which is one. The point zero, two gives four over one, which is four, so distance two. Mahalanobis distance counts how many standard deviations away you are, in the data's own geometry.

• d² = (x−μ)ᵀΣ⁻¹(x−μ)
• Distance in standard deviations

Now include correlation. With covariance zero point eight, the determinant is one minus zero point six four, which is zero point three six. So sigma inverse is one over zero point three six times one, minus zero point eight, minus zero point eight, one. For the point one, one the squared distance is about one point one one. For one, minus one it is exactly ten. Same Euclidean length, wildly different surprise.

• Σ = [[1, 0.8],[0.8, 1]]
• (1,1): d² ≈ 1.11
• (1,−1): d² = 10

Mahalanobis distance is just Euclidean distance after whitening, rotating along the eigenvectors and rescaling by the square roots of the eigenvalues so the ellipse becomes a circle. This makes it a classic tool for anomaly detection: points with a large squared distance are unlikely under the fitted Gaussian. One caution: the formula needs sigma to be invertible, so highly redundant features break it.

• Whitening: Σ becomes identity
• Large d² flags anomalies

Here is where it all leads. PCA, which we already met in lesson 23, simply finds the eigenvectors of this covariance matrix, and we will not re-derive it here. A Gaussian mixture, coming in lesson 44, fits a mean and covariance matrix to every cluster. And a Gaussian process in lesson 61 uses a covariance function to say how strongly any two inputs move together.

• PCA: eigenvectors of Σ
• GMM: one Gaussian per cluster
• GP: Σ from a kernel

A quick checklist of mistakes. Forgetting to center before using the matrix formula. Dividing by n when you meant n minus one. Trusting a single correlation number without plotting the data, since Anscombe-style datasets can hide anything. Reading zero correlation as independence. And reading any correlation as causation. Each of these appears in real projects, so keep this list handy.

• Forgetting to center
• Trusting r without a plot
• Causation from correlation

Let us recap. Covariance averages products of paired deviations, here one point five. Dividing by both standard deviations gives correlation, zero point seven seven five, which is the cosine of the angle between centered vectors. The covariance matrix is X centered transpose times X centered over n minus one, symmetric and positive semi-definite. The Gaussian turns it into ellipses, and Mahalanobis distance measures surprise inside them.

• Cov: paired deviations
• r = cosine of centered vectors
• Σ = XcᵀXc/(n−1), PSD
• Gaussian ellipse, Mahalanobis

Everything here used a sample to estimate a covariance, but how much can we trust such an estimate? That is exactly the question of lesson 39: sampling, estimation, confidence intervals and hypothesis tests.

• Lesson 39: estimation
• Confidence intervals, hypothesis tests
• Subscribe for the rest