Linear Algebra for AI

Every modern AI model is, at the level of the hardware, a long sequence of multiplications and additions on arrays of numbers. Linear algebra is the language for those arrays. When a search engine finds documents "similar" to your query, when a transformer mixes information between tokens, when LoRA fine-tunes a 70-billion-parameter model on one GPU, and when PyTorch throws RuntimeError: mat1 and mat2 shapes cannot be multiplied, you are looking at linear algebra. This lesson builds the parts an ML engineer uses daily, with the shapes, the numbers, and the code.

The picture to keep in your head

Think of a vector as an arrow, or equivalently as a list of measurements of one thing: a house described by (square metres, bedrooms, age), or a word described by 768 learned numbers. Think of a matrix as a machine that takes in arrows and puts out new arrows: it can stretch, rotate, squash, project, or mix the coordinates. A neural network layer is exactly such a machine (plus a bias and a nonlinearity). Training a network means adjusting the numbers inside the machines until the arrows that come out are useful.

Two questions recur in every section below: how big is it? (norms) and which way does it point? (directions, angles, eigenvectors). Keep both in mind.

Vectors

A vector $\mathbf{x} \in \mathbb{R}^n$ is an ordered list of $n$ real numbers $\mathbf{x} = (x_1, x_2, \dots, x_n)$. In ML, $n$ is the number of features (for tabular data) or the embedding dimension (for learned representations). We write vectors in bold lowercase and, by convention, treat them as columns.

Two operations define the vector space:

• Addition: $(\mathbf{x} + \mathbf{y})_i = x_i + y_i$. Residual connections in transformers are literally $\mathbf{x} + f(\mathbf{x})$.
• Scalar multiplication: $(c\,\mathbf{x})_i = c\,x_i$. A gradient-descent step is $\mathbf{w} - \eta \nabla L$: scale a vector, then subtract.

Norms: how big is a vector?

A norm measures length. The three you will meet constantly:

For $\mathbf{x} = (3, -4)$: $\|\mathbf{x}\|_1 = 7$, $\|\mathbf{x}\|_2 = \sqrt{9 + 16} = 5$, $\|\mathbf{x}\|_\infty = 4$.

Where they appear: the L2 norm of the weights is the L2 regularization (weight decay) penalty; the L1 norm encourages sparse weights (Lasso); gradient clipping rescales the gradient whenever its L2 norm exceeds a threshold, a standard safeguard in LLM training; and the L∞ norm bounds the perturbation size in adversarial-robustness work. For matrices, the Frobenius norm $\|A\|_F = \sqrt{\sum_{i,j} A_{ij}^2}$ is the L2 norm of all entries flattened.

The dot product and cosine similarity

The dot product (inner product) of two vectors of the same length is

where $\theta$ is the angle between them. The second equality is the geometric meaning: the dot product is large when the vectors are long and point the same way, zero when they are perpendicular (orthogonal), and negative when they point in opposite directions.

Dividing out the lengths leaves pure direction, the cosine similarity:

Worked example. Let $\mathbf{a} = (1, 2, 2)$ and $\mathbf{b} = (2, 1, 2)$.

1. Dot product: $1\cdot 2 + 2 \cdot 1 + 2 \cdot 2 = 2 + 2 + 4 = 8$.
2. Lengths: $\|\mathbf{a}\|_2 = \sqrt{1 + 4 + 4} = 3$, $\|\mathbf{b}\|_2 = \sqrt{4 + 1 + 4} = 3$.
3. Cosine similarity: $8 / (3 \cdot 3) = 8/9 \approx 0.889$, an angle of about $27°$.

Where it appears. Embedding models map text, images or audio to vectors so that semantically similar inputs land close together in direction. Retrieval-augmented generation (RAG) ranks documents by the cosine similarity between the query embedding and each document embedding. If the embeddings are L2-normalized (scaled to length 1), cosine similarity equals the plain dot product, which is why vector databases often store normalized vectors and use fast inner-product search. Attention in transformers scores how much token $i$ should attend to token $j$ using a dot product $\mathbf{q}_i^\top \mathbf{k}_j$ between a query and a key vector, divided by $\sqrt{d_k}$ so the scores' variance does not grow with the key dimension $d_k$.

Matrices as transformations

A matrix $A \in \mathbb{R}^{m \times n}$ has $m$ rows and $n$ columns. Multiplying it by a vector $\mathbf{x} \in \mathbb{R}^n$ gives $\mathbf{y} = A\mathbf{x} \in \mathbb{R}^m$ with

Two equivalent readings, both useful:

• Row view: each output $y_i$ is the dot product of row $i$ of $A$ with $\mathbf{x}$. A linear layer with $m$ output neurons computes $m$ dot products, one per neuron, each measuring how well the input matches that neuron's weight pattern.
• Column view: $\mathbf{y}$ is a weighted sum of the columns of $A$, with weights $x_1, \dots, x_n$. The set of all such outputs is the column space of $A$; its dimension is the rank.

A matrix is a linear map: $A(c\mathbf{x} + \mathbf{y}) = cA\mathbf{x} + A\mathbf{y}$. Rotations, scalings, reflections and projections are all linear. Stacking linear maps gives another linear map ($B(A\mathbf{x}) = (BA)\mathbf{x}$), which is precisely why neural networks need nonlinear activation functions between layers: without them, a 100-layer network would collapse to a single matrix.

Matrix multiplication and the shape rule

The product $C = AB$ of $A \in \mathbb{R}^{m \times k}$ and $B \in \mathbb{R}^{k \times n}$ is $C \in \mathbb{R}^{m \times n}$ with

The inner dimensions must match ($k$ and $k$), and the result takes the outer dimensions ($m \times n$). Matrix multiplication is associative, $(AB)C = A(BC)$, but generally not commutative: $AB \ne BA$, and often only one of them is even defined. Computing $AB$ costs about $2mkn$ floating-point operations (one multiply and one add per term), which is where the "2 × parameters FLOPs per token" rule of thumb for transformers comes from.

Worked example.

Shapes: $(2 \times 2)(2 \times 3) = (2 \times 3)$.

Batch × features × weights

In practice, data arrives as a design matrix $X \in \mathbb{R}^{B \times d_{\text{in}}}$: one row per example in the batch, one column per feature. A linear (fully connected) layer with weight matrix $W \in \mathbb{R}^{d_{\text{in}} \times d_{\text{out}}}$ and bias $\mathbf{b} \in \mathbb{R}^{d_{\text{out}}}$ computes

Read it as "each of the $B$ rows gets transformed independently by the same $W$." For a transformer MLP block with $B = 32$ sequences of $T = 512$ tokens and model width $d = 768$, the activations are $(32, 512, 768)$, the up-projection weight is $(768, 3072)$, and the output is $(32, 512, 3072)$: the leading dimensions just ride along. Note that PyTorch's nn.Linear stores its weight as $(d_{\text{out}}, d_{\text{in}})$ and computes x @ W.T + b; the maths is identical, only the storage convention differs.

The single most useful debugging habit in deep learning is writing the shape of every tensor next to every line of code.

Transpose

The transpose $A^\top$ swaps rows and columns: $(A^\top)_{ij} = A_{ji}$, so an $m \times n$ matrix becomes $n \times m$. Key identity: $(AB)^\top = B^\top A^\top$ (order reverses). Transposes are everywhere in backpropagation: if $Y = XW$, the gradient of the loss with respect to $X$ is $\frac{\partial L}{\partial Y} W^\top$ and with respect to $W$ is $X^\top \frac{\partial L}{\partial Y}$. You will derive this in the calculus lesson; here, notice that the transposes are exactly what make the shapes line up. A matrix with $A = A^\top$ is symmetric; covariance matrices and Hessians are symmetric.

Identity and inverse

The identity $I_n$ has ones on the diagonal and zeros elsewhere; $I\mathbf{x} = \mathbf{x}$. A square matrix $A$ is invertible if there is a matrix $A^{-1}$ with $A^{-1}A = AA^{-1} = I$. This exists exactly when the determinant $\det A \ne 0$, equivalently when $A$ has full rank, equivalently when no non-zero vector is squashed to zero.

For a $2 \times 2$ matrix,

Worked example. $A = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}$ has $\det A = 2\cdot1 - 1\cdot1 = 1$, so $A^{-1} = \begin{pmatrix} 1 & -1 \\ -1 & 2 \end{pmatrix}$. Check: $\begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} 1 & -1 \\ -1 & 2 \end{pmatrix} = \begin{pmatrix} 2-1 & -2+2 \\ 1-1 & -1+2 \end{pmatrix} = I$.

Where it appears: the closed-form least-squares solution for linear regression is $\mathbf{w} = (X^\top X)^{-1} X^\top \mathbf{y}$ (the normal equations). In code you almost never form an inverse explicitly; you call a solver (np.linalg.solve, np.linalg.lstsq), which is faster and far more numerically stable. When $X^\top X$ is nearly singular (highly correlated features), adding $\lambda I$ before solving is ridge regression, a direct link between linear algebra and regularization.

Eigenvalues and eigenvectors

Most vectors change direction when a matrix acts on them. An eigenvector of a square matrix $A$ is a non-zero vector that only gets scaled:

where the scalar $\lambda$ is the eigenvalue. Eigenvectors are the matrix's "natural axes"; eigenvalues say how much it stretches along each.

Worked example. $S = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$. Eigenvalues solve $\det(S - \lambda I) = 0$:

For $\lambda_1 = 3$: $(S - 3I)\mathbf{v} = \mathbf{0}$ gives $-v_1 + v_2 = 0$, so $\mathbf{v}_1 = \tfrac{1}{\sqrt2}(1, 1)$. Check: $S(1,1)^\top = (3, 3)^\top = 3\,(1,1)^\top$. For $\lambda_2 = 1$: $\mathbf{v}_2 = \tfrac{1}{\sqrt2}(1, -1)$. So $S$ stretches the diagonal direction by 3 and leaves the anti-diagonal alone.

Symmetric matrices are especially well behaved (the spectral theorem): their eigenvalues are real and their eigenvectors can be chosen orthonormal, so $S = Q \Lambda Q^\top$ with $Q$ orthogonal and $\Lambda$ diagonal.

Where eigenvalues appear in AI:

• PCA: the eigenvectors of the data's covariance matrix are the directions of greatest variance; the eigenvalues are those variances.
• Optimization: the eigenvalues of the Hessian (the matrix of second derivatives) describe the curvature of the loss. A large ratio between the largest and smallest eigenvalue (a high condition number) means a long narrow valley where plain gradient descent zig-zags; this motivates momentum and Adam (lesson 5).
• Recurrent networks: repeatedly multiplying by the same matrix makes components along eigenvectors with $|\lambda| > 1$ explode and those with $|\lambda| < 1$ vanish, the classic exploding/vanishing gradient problem.

Singular value decomposition (SVD)

Eigendecomposition only works for (some) square matrices. The SVD works for every matrix $A \in \mathbb{R}^{m \times n}$:

where $U \in \mathbb{R}^{m \times m}$ and $V \in \mathbb{R}^{n \times n}$ are orthogonal (their columns are orthonormal) and $\Sigma$ is $m \times n$ with non-negative singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$ on its diagonal. Geometrically: every linear map is a rotation (possibly with a reflection) by $V^\top$, then an axis-aligned stretch by $\Sigma$, then another rotation (or reflection) by $U$. The number of non-zero singular values is the rank.

Worked check: for $M = \begin{pmatrix} 3 & 0 \\ 4 & 5 \end{pmatrix}$, $M^\top M = \begin{pmatrix} 25 & 20 \\ 20 & 25 \end{pmatrix}$ has eigenvalues $45$ and $5$, so the singular values are $\sqrt{45} \approx 6.708$ and $\sqrt{5} \approx 2.236$. (In general the singular values of $A$ are the square roots of the eigenvalues of $A^\top A$.)

Low-rank approximation

Keep only the top $r$ singular values and vectors:

The Eckart–Young theorem says $A_r$ is the best rank-$r$ approximation of $A$ in Frobenius (and spectral) norm. Storing $A_r$ costs $r(m + n)$ numbers instead of $mn$.

PCA through the SVD

Principal component analysis centers the data matrix $X \in \mathbb{R}^{N \times d}$ (subtract each column's mean to get $X_c$), takes the SVD $X_c = U\Sigma V^\top$, and projects onto the first $k$ right singular vectors: $Z = X_c V_{:, 1:k}$. The variance explained by component $i$ is $\sigma_i^2 / (N - 1)$. This is mathematically the same as the eigendecomposition of the covariance matrix $\frac{1}{N-1} X_c^\top X_c$, but numerically better. PCA is used for visualizing embeddings, compressing features, and whitening inputs.

LoRA: low rank in fine-tuning

LoRA (Hu et al., 2021) freezes a pretrained weight matrix $W_0 \in \mathbb{R}^{d \times k}$ and learns an update in factored form:

The product $BA$ has rank at most $r$. For a $4096 \times 4096$ projection, a full update has $16{,}777{,}216$ parameters; with $r = 8$, LoRA trains $8 \times (4096 + 4096) = 65{,}536$, a factor of 256 fewer for that matrix. The bet, which the paper supports empirically, is that the change needed to adapt a model to a task lives in a low-dimensional subspace. The authors report that for GPT-3 175B, LoRA cut trainable parameters by 10,000 times and GPU memory by 3 times compared with full fine-tuning with Adam. Since $B$ is initialized to zero, $\Delta W = 0$ at the start, so training begins exactly at the pretrained model.

Tensors and broadcasting

A tensor, in deep-learning usage, is simply an $n$-dimensional array (this is looser than the physics/maths definition). Scalars are rank-0, vectors rank-1, matrices rank-2. Typical shapes: an image batch in PyTorch is $(B, C, H, W)$; token activations are $(B, T, d)$; attention scores are $(B, h, T, T)$ for $h$ heads.

Broadcasting lets you combine arrays of different shapes without copying data. NumPy and PyTorch use the same rule: align shapes from the rightmost dimension; two dimensions are compatible if they are equal or one of them is 1 (a missing dimension counts as 1); size-1 dimensions are virtually stretched.

The last row is the classic bug: to subtract a per-row value you need shape $(32, 1)$, e.g. x - x.max(axis=1, keepdims=True).

Code: the whole lesson in NumPy

import numpy as np

# Cosine similarity between "embeddings"
a, b = np.array([1., 2., 2.]), np.array([2., 1., 2.])
cos = a @ b / (np.linalg.norm(a) * np.linalg.norm(b))
print(round(cos, 3))                        # 0.889

# Batch x features @ weights
rng = np.random.default_rng(0)
X = rng.normal(size=(32, 768))             # 32 examples, 768 features
W = rng.normal(size=(768, 3072)) * 0.02
bias = np.zeros(3072)
Y = X @ W + bias                           # bias broadcasts over rows
print(Y.shape)                             # (32, 3072)

# Solve instead of invert
A = np.array([[2., 1.], [1., 1.]])
print(np.linalg.solve(A, np.array([3., 2.])))   # [1. 1.]

# Eigen-decomposition of a symmetric matrix
vals, vecs = np.linalg.eigh(np.array([[2., 1.], [1., 2.]]))
print(vals)                                # [1. 3.]

# PCA via SVD on correlated 2-D data
Z = rng.normal(size=(500, 2)) @ np.array([[3., 0.], [2.5, 0.5]])
Zc = Z - Z.mean(axis=0)
U, S, Vt = np.linalg.svd(Zc, full_matrices=False)
explained = S**2 / (S**2).sum()
print(explained.round(3))                  # first component dominates
proj = Zc @ Vt[:1].T                       # (500, 1) projection

# Low-rank (LoRA-style) parameter count
d = k = 4096; r = 8
print(d * k, r * (d + k))                  # 16777216 65536

np.linalg.eigh is the routine for symmetric matrices (it returns eigenvalues in ascending order); np.linalg.svd with full_matrices=False returns the compact SVD. The same operations exist in PyTorch as torch.linalg.eigh, torch.linalg.svd and the @ operator, and broadcasting rules are identical.

Common pitfalls

• Shape mismatches hidden by broadcasting. Subtracting a $(32,)$ vector from a $(32, 1)$ column produces a $(32, 32)$ matrix with no error. Assert shapes in critical code.
• Confusing row and column conventions. Textbooks often write $W\mathbf{x}$ with column vectors; code uses row-major batches $XW$ or PyTorch's x @ W.T. Translate deliberately.
• Computing explicit inverses. inv(A) @ b is slower and less accurate than solve(A, b); for least squares use lstsq.
• Cosine on unnormalized vectors with dot-product indexes. If your vector index uses inner product but your vectors are not unit length, longer vectors win regardless of direction.
• Forgetting to center before PCA. Without centering, the first "component" mostly points at the mean.
• Elementwise vs matrix product. * is elementwise in NumPy and PyTorch; @ is matrix multiplication.

In practice

• Write shapes as comments; use einsum (np.einsum('btd,df->btf', X, W)) when an operation involves more than two axes, because the subscripts document themselves.
• Normalize embeddings once at indexing time and use inner-product search for cosine similarity.
• Use float32 or lower-precision types for large matrix products, but accumulate reductions in higher precision; mixed precision is standard on GPUs.
• Reach for SVD when you need to understand or compress a matrix: checking effective rank of a weight update, PCA plots of embeddings, or choosing a LoRA rank.
• When a model diverges, look at norms: the L2 norm of gradients and of activations per layer is the first diagnostic professionals check, and gradient clipping by global norm is a default safety net.