Matrix Algebra in Data Analysis

Info sheet · Statistics for Psychology & Neuroscience

Author

Andrew Bell

Published

August 13, 2026

Info sheet 0.2 (draft) · Prerequisites: matrix multiplication, transpose & inverse; multiple regression · Give feedback ↗

Working notes for the author — not shown to students once collapsed; remove before publishing.

  • [ ]

What you’ll get from this sheet

The payoff. Matrix algebra is the engine behind tools you already use — knowing it lets you see what R, Python, or MATLAB are doing behind the scenes, speed things up, and spot problems. By the end you should be able to:

  1. Build a covariance (and correlation) matrix with one matrix product.
  2. Solve multiple regression with the normal equation.
  3. Recognise vectorisation, and say why it matters for speed.

Three everyday tools are matrix algebra: the covariance matrix \(S = \frac{1}{n-1}X^{\mathsf T}X\) (after centring), the regression coefficients \(\theta = (X^{\mathsf T}X)^{-1}X^{\mathsf T}y\) (the normal equation), and vectorisation — operating on whole arrays at once instead of looping.

The covariance matrix in one product

Covariance measures the joint variability of two variables — positive if they rise together, negative if one rises as the other falls, near zero if unrelated. Computing every pairwise covariance by hand is tedious; matrix algebra does the lot at once. Centre the columns of your data matrix \(X\) (subtract each column’s mean), and then

\[S = \frac{1}{n-1}\,X^{\mathsf T}X\]

is the whole covariance matrix — a \(p \times p\) grid of every predictor’s variance (diagonal) and covariances (off‑diagonal). Standardise instead of just centring and the same product gives the correlation matrix. (This is exactly the matrix PCA then decomposes.)

Regression is one matrix line

Multiple regression finds the coefficients \(\theta\) that minimise squared error. Written in matrix form, the whole least‑squares solution is a single expression — the normal equation:

\[\theta = (X^{\mathsf T}X)^{-1}X^{\mathsf T}y\]

where \(X\) is your predictors with a column of 1s bolted on for the intercept, and \(y\) is the outcome. That’s it — transpose, multiply, invert, multiply. Below, a slider adds noise to some data; watch the normal equation recompute the intercept and slope (identical to what lm() returns).

With no noise the normal equation recovers the true intercept (2) and slope (1.5) almost exactly; add noise and the coefficients wobble around them, but the machinery never changes — transpose \(X\), multiply, invert the little \(X^{\mathsf T}X\), multiply again. Every regression you’ve ever run is doing this behind the scenes.

Vectorisation: stop looping

One last idea, about speed. Non‑vectorised code loops through elements one at a time; vectorised code applies an operation to a whole array at once. Squaring a million numbers with a for loop is slow; writing data**2 (which hands the whole array to optimised, compiled routines) is often orders of magnitude faster. The rule of thumb: when you spot a loop over a big array, ask can this be vectorised? — the matrix formulations above are vectorisation in action.

Try it: vectorise this

The cell below computes, with a loop, how far each song sits from the average danceability. Your turn: replace the loop with a single whole‑array expression, and check the timing drops. (It runs live — edit and re‑run.)

start = time.time()
out_vec = (danceability - danceability.mean()) ** 2   # whole array at once
print(f"vectorised: {time.time() - start:.4f}s")
print("match:", np.allclose(out, out_vec))            # True — same numbers, a fraction of the time

The vectorised version hands the entire array to NumPy’s compiled routines instead of stepping through a million Python iterations — typically tens to hundreds of times faster, for identical output.

See it in code

x = (1:8)';  y = 2 + 1.5*x + randn(8,1);
X = [ones(8,1) x];
theta = (X' * X) \ (X' * y);   % normal equation ( \ is the stable solve )

S = cov(data);                 % covariance matrix
R = corrcoef(data);            % correlation matrix

The R and Python tabs run live; MATLAB is a static reference.

What does \(\theta = (X^{\mathsf T}X)^{-1}X^{\mathsf T}y\) give you — and what must be true of \(X^{\mathsf T}X\)?

It gives the least‑squares regression coefficients (the intercept and slopes), the same numbers lm() reports. For it to work, \(X^{\mathsf T}X\) must be invertible (non‑singular). That fails under perfect multicollinearity — if one predictor is an exact linear combination of the others, \(X^{\mathsf T}X\) is singular and the inverse doesn’t exist, which is why you can’t put perfectly redundant predictors in a regression.

The normal equation is perfect for understanding regression, but forming \((X^{\mathsf T}X)^{-1}\) explicitly is numerically unstable and blows up when predictors are highly collinear. Real software doesn’t invert — it uses QR decomposition or a direct solver (R’s lm, MATLAB’s \). Treat the formula as the idea, not the implementation. And remember covariance needs its columns centred first, or you’ll get nonsense.

Where this shows up next

This closes the matrix‑algebra appendix — and connects it back to the rest of the book: the covariance matrix is what PCA decomposes, and the normal equation is OLS in one line. See Appendix (Matrix Algebra) and Chapter (Regression) for the full treatments.