Linear Algebra for ML

Why Linear Algebra for ML?

Every dataset is a matrix. Every prediction is a dot product. Every model update is a matrix operation. Linear algebra isn't a prerequisite you check off — it's the language ML is written in. This lesson covers what you'll actually use, with Python implementations throughout.

Vectors and Matrices in ML

Data matrix $\mathbf{X}$Shape $(m imes n)$: $m$ samples, $n$ features. Each row is one example.

Weight vector $\mathbf{w}$Shape $(n,)$: one weight per feature. The "learned parameters."

Prediction$\hat{\mathbf{y}} = \mathbf{X}\mathbf{w}$, shape $(m,)$ — one prediction per sample.

Python

import numpy as np

# Simulate a dataset: 5 samples, 3 features
X = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9],
              [2, 0, 1],
              [3, 3, 3]], dtype=float)

w = np.array([0.5, -0.2, 0.8])   # weights
b = 1.0                            # bias (intercept)

# Linear prediction — the core of every linear model
y_hat = X @ w + b                  # (5,)
print("Predictions:", y_hat.round(2))

Matrix Multiplication

Matrix Multiply$$c_{ij} = \sum_{k} a_{ik}\, b_{kj} \qquad \mathbf{A}_{m imes k}\,\mathbf{B}_{k imes n} = \mathbf{C}_{m imes n}$$Inner dimensions must match. The result shape is (outer rows × outer cols).

Python

A = np.array([[1, 2], [3, 4], [5, 6]])   # (3,2)
B = np.array([[7, 8, 9], [10, 11, 12]])  # (2,3)
C = A @ B                                 # (3,3)

# Properties
print(A.T @ A)            # Gram matrix — always square, symmetric, PSD
print(np.allclose(A.T.T, A))  # double transpose = original

# Neural network forward pass: 3 layers
W1 = np.random.randn(4, 3)   # hidden: 4 neurons, 3 inputs
W2 = np.random.randn(2, 4)   # output: 2 neurons, 4 inputs
x  = np.random.randn(3)      # single input
h  = np.maximum(0, W1 @ x)   # ReLU activation
y  = W2 @ h                   # final output

Transpose and Symmetric Matrices

$(A^T)_{ij} = A_{ji}$ — rows become columns. Key identities: $(AB)^T = B^T A^T$, $(A^T)^T = A$.

A matrix is symmetric if $A = A^T$. Covariance matrices, Gram matrices, and Hessians are always symmetric — a property exploited by eigensolver algorithms.

Matrix Inverse

$A^{-1}A = AA^{-1} = I$. Exists only for square, full-rank matrices. The analytical solution to linear regression (Normal Equation) requires it:

Normal Equation$$\mathbf{w}^* = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$$

Python

np.random.seed(0)
X = np.random.randn(100, 3)
X = np.c_[np.ones(100), X]       # add bias column
true_w = np.array([1, 2, -1, 0.5])
y = X @ true_w + np.random.randn(100) * 0.3

# Normal equation
w_star = np.linalg.solve(X.T @ X, X.T @ y)   # NEVER use inv() — use solve()
print("Recovered weights:", w_star.round(3))   # should be close to true_w

⚠️

Use np.linalg.solve(A, b), never np.linalg.inv(A) @ b. solve() uses LU decomposition internally — 3× faster and numerically stabler.

Norms and Distance

L^p$ Norm$$\|\mathbf{x}\|_p = \left(\sum_i |x_i|^p ight)^{1/p} \qquad ext{L1: Manhattan}\quad ext{L2: Euclidean}\quad ext{L}\infty ext{: max element}$

Python

x = np.array([3.0, -4.0, 0.0])
print(np.linalg.norm(x, 1))     # L1 = 7.0   (Lasso regulariser)
print(np.linalg.norm(x, 2))     # L2 = 5.0   (Ridge regulariser)
print(np.linalg.norm(x, np.inf))# L∞ = 4.0

# Euclidean distance matrix (efficient)
from scipy.spatial.distance import cdist
X = np.random.randn(50, 10)
D = cdist(X, X, metric="euclidean")    # (50,50) pairwise distances

Eigenvalues & PCA Preview

If $A\mathbf{v} = \lambda\mathbf{v}$, $\mathbf{v}$ is an eigenvector with eigenvalue $\lambda$. For the covariance matrix $\Sigma = X^TX/(m-1)$, eigenvectors are the principal components — directions of maximum variance.

Python

X = np.random.randn(200, 4)
X_c = X - X.mean(axis=0)                 # centre
Sigma = X_c.T @ X_c / (len(X_c) - 1)    # covariance matrix (4,4)

vals, vecs = np.linalg.eigh(Sigma)        # eigh for symmetric matrices
order = np.argsort(vals)[::-1]            # sort descending
print("Explained variance ratio:", (vals[order] / vals.sum()).round(3))

X_pca = X_c @ vecs[:, order[:2]]         # project to 2D

Summary

Datasets are matrices: $\mathbf{X}$ is $(m imes n)$. Predictions are $\hat{\mathbf{y}} = \mathbf{X}\mathbf{w}$.
Matrix multiply: inner dims must match. Use @ operator in Python.
Use np.linalg.solve() not inv() for solving linear systems.
L1 norm → Lasso; L2 norm → Ridge. L∞ = max absolute element.
Eigendecomposition of the covariance matrix gives PCA components.