Linear Algebra

The math of arrows, grids, and lists of numbers — and the single most important subject for anyone who wants to understand modern machine learning. By the end of this page you should be able to read a matrix expression, know what it does geometrically, and see why every line of a neural network is quietly doing linear algebra.

1. Why linear algebra

Here is the headline: linear algebra is the language every modern ML model, every 3D graphics pipeline, and every signal processing system is written in. When someone says "GPT-4 has 1.7 trillion parameters," they mean that a collection of matrices has 1.7 trillion numbers in it. When someone says "the embedding of dog," they mean a vector. When someone says "attention," they mean dot products. When someone says "a forward pass of a neural network," they mean a sequence of matrix multiplications with some non-linear squashing in between.

If you only remember three sentences from this page, remember these:

• Vectors are lists of numbers with direction and length. Every word, every image, every user in a recommender system eventually becomes a vector.
• Matrices are functions that send vectors to other vectors, in a very restricted way (they preserve lines and the origin). "Restricted" is what makes them tractable and fast.
• The dot product measures alignment. Two vectors that point the same way have a big dot product. Two that are perpendicular have zero. Attention is literally this, scaled.

If you already know what a vector is, skim sections 3–6 and dive into eigenvalues and SVD. If you've never seen any of this — or haven't seen it since school — read straight through. Nothing here assumes calculus; we'll point to calculus.html only where we mention derivatives, and you can ignore those footnotes on a first pass.

2. Vocabulary cheat-sheet

The notation is ugly and inconsistent across textbooks. Here is the set of conventions this page uses, and which we recommend you steal.

If a symbol appears on this page and you can't find it in this table, we'll define it inline where it appears. If we don't, report the bug.

3. Vectors

A vector is two things that happen to be the same thing:

• An ordered list of numbers, like $(3, -1, 4, 0.5)$. Programmers know this as a 1D array.
• An arrow in space with a length and a direction, starting at the origin.

When you do linear algebra, you want both pictures in your head at once. The list is how you compute; the arrow is how you reason.

A 2D vector $\mathbf{v} = (v_1, v_2)$ is an arrow from $(0,0)$ to the point $(v_1, v_2)$ in the plane. A 3D vector $(v_1, v_2, v_3)$ is an arrow in space. A 768-dimensional vector — the default size of a GPT embedding — is an arrow in a space you cannot visualize, but that still obeys the same rules. That's the whole trick: once you know how 2D works, the algebra carries you up into any dimension for free.

Vector addition

You add vectors component-by-component:

Scalar multiplication

Multiplying a vector by a plain number (a scalar) stretches or squashes it:

These two operations — addition and scalar multiplication — are the only operations that define what "linear" means. Everything else in this subject is built out of them.

Norms: how long is a vector?

A norm is any reasonable way to measure the length of a vector. There are several, and they show up in different ML contexts.

The default, and the one you already know from Pythagoras, is the Euclidean norm or $\ell_2$ norm:

The $\ell_1$ norm adds absolute values instead:

And the $\ell_\infty$ norm takes the max:

All three answer the question "how big is this vector?" but give different answers. Pick the one that fits your problem. For most things in ML, $\ell_2$ is the default.

4. Dot products

The dot product is the single most important operation in this subject. Learn it and you've learned 40% of linear algebra.

There are two definitions. They're equal. The trick is that each definition tells you something different about why it matters.

The algebraic definition

The geometric definition

The fact that these two formulas are equal is not obvious. It follows from the law of cosines, and it's essentially the entire reason linear algebra works: the thing you can compute (the sum) equals the thing that has geometric meaning (the angle × length).

Rearranging:

Projections

If $\mathbf{u}$ is a unit vector (length 1), the dot product $\mathbf{a}\cdot\mathbf{u}$ is the length of the shadow $\mathbf{a}$ casts onto the line through $\mathbf{u}$. This is the projection of $\mathbf{a}$ onto $\mathbf{u}$:

Why does attention appear here? In a transformer, the "score" between a query and a key is literally $\mathbf{q}\cdot\mathbf{k}$. High dot product → "these two tokens are relevant to each other" → high attention weight. That's the whole idea. The softmax is window dressing.

5. Matrices

A matrix is a rectangular grid of numbers. If it has $m$ rows and $n$ columns we say it is $m \times n$ and write $A \in \mathbb{R}^{m\times n}$. The entry in row $i$, column $j$ is written $A_{ij}$.

Matrix as function: a vector-transforming machine

Here is the point of view that turns all of this into something meaningful, rather than a spreadsheet exercise:

And not just any function — a linear function. "Linear" means two properties hold:

How do you know a matrix when it's given to you as a function? Apply it to the standard basis vectors and read off the answers:

6. Matrix multiplication

Matrix multiplication is the operation people get most wrong most often. We'll show two pictures, both equally valid. Use whichever makes more sense at the time.

The row · column picture

If $A$ is $m \times k$ and $B$ is $k \times n$, the product $C = AB$ is $m \times n$, with entries:

The column picture (this is the one to remember)

Here is the much more useful way to see it. If $B$ has columns $\mathbf{b}_1, \dots, \mathbf{b}_n$, then:

A worked example

Let's compute a small one by hand. Take:

Row × column:

• $C_{11}$ = row 1 of $A$ · column 1 of $B$ = $(1,2) \cdot (5,7) = 5+14 = 19$
• $C_{12}$ = row 1 of $A$ · column 2 of $B$ = $(1,2) \cdot (6,8) = 6+16 = 22$
• $C_{21}$ = row 2 of $A$ · column 1 of $B$ = $(3,4) \cdot (5,7) = 15+28 = 43$
• $C_{22}$ = row 2 of $A$ · column 2 of $B$ = $(3,4) \cdot (6,8) = 18+32 = 50$

Now compute $BA$ and watch it come out different:

What does multiply (associativity and distributivity)

Matrix multiplication is associative and distributive, but not commutative:

7. Interactive: 2×2 transformations

Time to build intuition with your hands. Below is a 2D playground. Adjust the four entries of a 2×2 matrix $M$ with the sliders. The cyan outline is the original unit square with its basis vectors drawn; the violet outline is the image of the unit square under $M$ — what it becomes after you apply the transformation. You'll see the columns of $M$ become the new positions of $\hat x$ and $\hat y$. That is not a coincidence; that's section 5.

Watch for:

• What a rotation looks like: try $\begin{pmatrix}\cos 30° & -\sin 30°\\ \sin 30° & \cos 30°\end{pmatrix} \approx \begin{pmatrix}0.87 & -0.5\\ 0.5 & 0.87\end{pmatrix}$.
• What a scale looks like: set $M = \begin{pmatrix}2 & 0\\0 & 0.5\end{pmatrix}$.
• What a shear looks like: set $M = \begin{pmatrix}1 & 1\\0 & 1\end{pmatrix}$.
• A collapse: $M = \begin{pmatrix}1 & 1\\2 & 2\end{pmatrix}$ — notice the unit square becomes a line segment. The determinant is 0. Rank dropped.

8. Special matrices

Most matrices you meet in the wild are one of a handful of types. Naming them matters, because the algorithms that operate on them are specialized to the type.

Identity matrix $I$

The matrix that does nothing: $I\mathbf{x} = \mathbf{x}$ for any $\mathbf{x}$. It has 1s on the diagonal and 0s everywhere else. The multiplicative identity: $AI = IA = A$.

Zero matrix $\mathbf{0}$

All entries 0. Sends every vector to the zero vector. The additive identity.

Diagonal matrix

Only the diagonal entries can be non-zero. A diagonal matrix $D$ scales each coordinate independently:

Triangular matrix

All entries above (or all below) the diagonal are 0. Upper triangular looks like:

Symmetric matrix

A matrix is symmetric if $A = A^T$, i.e. $A_{ij} = A_{ji}$.

Orthogonal matrix

A square matrix $Q$ is orthogonal if $Q^T Q = I$, i.e. $Q^T = Q^{-1}$. Equivalently, its columns are mutually perpendicular unit vectors.

Rotation matrix

In 2D, the matrix that rotates by angle $\theta$ counterclockwise is:

Permutation matrix

An identity matrix with its rows rearranged. When you multiply $P\mathbf{x}$, you get $\mathbf{x}$ with its components permuted. Permutation matrices are orthogonal.

Projection matrix

A square matrix $P$ with $P^2 = P$. Applied twice, it does the same thing as once — because after the first application you're already in the target subspace. The "shine a light onto this subspace" operator.

9. Linear systems $A\mathbf{x}=\mathbf{b}$

The original motivation for linear algebra, and still the most common practical problem: you're given a matrix $A$ and a vector $\mathbf{b}$, and you want to find the vector $\mathbf{x}$ such that $A\mathbf{x} = \mathbf{b}$.

Three outcomes

There are exactly three things that can happen when you try to solve $A\mathbf{x} = \mathbf{b}$:

1. Unique solution. There's exactly one $\mathbf{x}$. This is the best case. It happens when $A$ is square and "full rank" (see §11).
2. No solution. The equation is inconsistent. Geometrically: $\mathbf{b}$ is not in the span of $A$'s columns — there's no way to combine the columns to hit it. Think "three equations, two unknowns, incompatible."
3. Infinitely many solutions. There's at least one $\mathbf{x}$, but you can slide it along some direction and stay a solution. This happens when $A$ has fewer effective constraints than unknowns.

Geometrically, each row of $A\mathbf{x}=\mathbf{b}$ is a constraint ("stay on this hyperplane"). The solution set is the intersection of all the constraints. Three planes in 3D generically meet at a point; if two are parallel and different you get no solution; if all three coincide you get a whole plane of solutions.

Gaussian elimination

The canonical algorithm for solving $A\mathbf{x}=\mathbf{b}$, going back 2000+ years in China and named after Gauss in the 1800s. The idea is boring and unstoppable: add multiples of one equation to another until the system is upper triangular, then back-substitute.

Three operations are allowed (they don't change the solution set):

• Swap two rows.
• Multiply a row by a non-zero scalar.
• Add a multiple of one row to another.

The resulting form is called row echelon form. If you also clean up the upper part to get 1s on the pivots and 0s above them, it's reduced row echelon form (RREF). The row of a pivot tells you which variables are "leading"; the rest are free parameters.

10. Vector spaces

We've been throwing around the word "vector" informally. Here's the formal version. A vector space is a set $V$ with two operations — addition and scalar multiplication — that obey the eight rules you'd expect (associativity, commutativity of addition, distributivity, existence of a zero element, etc.). For our purposes $V = \mathbb{R}^n$ is the only example that matters.

Span

The span of a set of vectors is the set of all their linear combinations:

Linear independence

A set of vectors $\{\mathbf{v}_1, \dots, \mathbf{v}_k\}$ is linearly independent if the only way to get the zero vector as a combination is to use all zero coefficients:

Basis and dimension

A basis of a vector space is a linearly independent set that spans the whole space. In $\mathbb{R}^n$, the "standard basis" is the set of unit vectors along each axis: $\mathbf{e}_1 = (1,0,\dots,0)$, $\mathbf{e}_2 = (0,1,\dots,0)$, and so on. Any vector in $\mathbb{R}^n$ is a unique linear combination of these.

The dimension of a vector space is the number of vectors in any basis. For $\mathbb{R}^n$ that's $n$. Different bases can have different vectors — they just always have the same count.

Why care? Because the whole point of choosing a basis is that expressing a vector or a matrix in the right basis can make a hard problem trivial. That's the entire content of diagonalization: find a basis in which $A$ looks like a diagonal matrix.

11. Rank and nullspace

Given a matrix $A$, there are four natural subspaces associated with it. They are called the four fundamental subspaces (Gil Strang's phrase) and they tell you everything about what $A$ does.

• Column space $\text{Col}(A)$ — the span of the columns of $A$. This is the set of possible outputs of $A\mathbf{x}$. If $\mathbf{b}$ is in the column space, $A\mathbf{x} = \mathbf{b}$ has a solution; otherwise it doesn't.
• Null space $\text{Null}(A)$ — the set of $\mathbf{x}$ with $A\mathbf{x} = \mathbf{0}$. These are the input vectors the matrix annihilates.
• Row space $\text{Row}(A) = \text{Col}(A^T)$.
• Left null space $\text{Null}(A^T)$.

The rank of $A$ is the dimension of its column space — equivalently, the number of linearly independent columns. Equivalently again, the number of linearly independent rows. That these are all equal is already a small miracle.

Rank–nullity theorem

A square $n\times n$ matrix is full rank (rank $n$) if and only if its null space is just $\{\mathbf{0}\}$, if and only if it is invertible, if and only if its determinant is non-zero, if and only if $A\mathbf{x}=\mathbf{b}$ has a unique solution for every $\mathbf{b}$. These are all the same statement.

12. Determinants

The determinant $\det(A)$ of a square matrix is a single number that says "by how much does $A$ scale oriented volume?" That's it. That's the whole story.

For a 2×2 matrix:

For a 3×3 matrix, the determinant is the signed volume of the parallelepiped spanned by the columns. For an $n\times n$ matrix, it's the signed $n$-dimensional volume of the $n$-parallelepiped spanned by the columns. Every textbook defines it differently — by cofactor expansion, by permutations, by LU decomposition — but they all compute the same number, and the number always means this.

Key properties

• $\det(I) = 1$ — the identity scales volume by 1.
• $\det(AB) = \det(A)\det(B)$ — composing transformations multiplies their volume factors.
• $\det(A^T) = \det(A)$.
• $\det(A) = 0 \iff A$ is singular $\iff$ columns are linearly dependent $\iff$ rank < $n$.
• $\det(A^{-1}) = 1/\det(A)$.

If you see $\det(A) = 0$ anywhere in ML, something collapsed. Usually it's a signal to reach for a pseudoinverse or add regularization.

13. Eigenvalues and eigenvectors

Pick any square matrix $A$. Most vectors, when you apply $A$, come out pointing in a different direction. But for some special vectors, $A$ just stretches them — it doesn't rotate them. Those special directions are eigenvectors, and the stretch factors are eigenvalues.

How to find them by hand

Rewrite the eigenvalue equation:

So the eigenvalues are the roots of the characteristic polynomial:

Worked example: a 2×2 case

Take $A = \begin{pmatrix} 4 & 1 \\ 2 & 3 \end{pmatrix}$. Then:

Solve the quadratic: $\lambda^2 - 7\lambda + 10 = 0$ factors as $(\lambda - 5)(\lambda - 2)$. So $\lambda_1 = 5$, $\lambda_2 = 2$. To find the eigenvectors, plug each $\lambda$ back into $(A-\lambda I)\mathbf{v}=\mathbf{0}$:

• For $\lambda = 5$: $(A - 5I)\mathbf{v} = \begin{pmatrix}-1 & 1\\ 2 & -2\end{pmatrix}\mathbf{v} = 0$. The rows are multiples of each other, so one equation is enough: $-v_1 + v_2 = 0$, hence $v_2 = v_1$. Any vector along $(1,1)$ works. Take $\mathbf{v}_1 = (1,1)$.
• For $\lambda = 2$: $(A - 2I)\mathbf{v} = \begin{pmatrix}2&1\\2&1\end{pmatrix}\mathbf{v} = 0$ gives $2v_1 + v_2 = 0$, i.e. $v_2 = -2v_1$. Take $\mathbf{v}_2 = (1,-2)$.

Check: $A\mathbf{v}_1 = \begin{pmatrix}4+1\\2+3\end{pmatrix} = \begin{pmatrix}5\\5\end{pmatrix} = 5\mathbf{v}_1$. Good. $A\mathbf{v}_2 = \begin{pmatrix}4-2\\2-6\end{pmatrix} = \begin{pmatrix}2\\-4\end{pmatrix} = 2\mathbf{v}_2$. Good.

Once you've done this a few times, you'll never want to do it for anything bigger than 3×3 again. Everything larger gets computed numerically.

14. Diagonalization

Suppose $A$ is a square matrix with $n$ linearly independent eigenvectors $\mathbf{v}_1, \dots, \mathbf{v}_n$ and eigenvalues $\lambda_1, \dots, \lambda_n$. Pack the eigenvectors into the columns of a matrix $P$ and the eigenvalues onto the diagonal of a matrix $D$. Then:

Not every matrix is diagonalizable — if it has fewer than $n$ independent eigenvectors it isn't, and you need Jordan normal form. But every symmetric matrix is diagonalizable by an orthogonal matrix (the spectral theorem), and even matrices that aren't diagonalizable over the reals can be diagonalized in the generalized sense of SVD — see the next section.

Diagonalization is why eigenvalues matter in practice. If $A = PDP^{-1}$, then:

15. Orthogonal bases and Gram–Schmidt

A basis is orthogonal if all pairs of basis vectors are perpendicular: $\mathbf{v}_i \cdot \mathbf{v}_j = 0$ for $i \ne j$. It's orthonormal if additionally each basis vector has length 1. Orthonormal bases are magical: once you have one, everything is easy.

In an orthonormal basis $\{\mathbf{q}_1, \dots, \mathbf{q}_n\}$, any vector $\mathbf{x}$ is:

Gram–Schmidt is the recipe for turning an arbitrary linearly independent set $\{\mathbf{a}_1, \dots, \mathbf{a}_k\}$ into an orthonormal set $\{\mathbf{q}_1, \dots, \mathbf{q}_k\}$ spanning the same space. You process the vectors one at a time, subtracting off the part that points along the directions you've already made orthonormal:

16. The Singular Value Decomposition

Now the big one. The SVD is the most important factorization in linear algebra and, by extension, in machine learning. It is to linear algebra what the Taylor series is to calculus: it takes any object of the subject and writes it in a standard, useful form.

Why SVD is the most important factorization in ML

Four huge applications, all of which reduce to "look at the SVD":

• PCA (Principal Component Analysis). Given a data matrix $X$, the top singular vectors of (a centered version of) $X$ point in the directions of maximum variance. Dimensionality reduction, denoising, feature engineering — all this, from SVD. See next section.
• Low-rank approximation. The best rank-$k$ approximation of $A$ in both Frobenius and spectral norm is the one obtained by keeping only the top $k$ singular values. This is Eckart–Young. Low-rank approximations are how you compress weight matrices (think LoRA), how you do matrix completion (Netflix), and how you regularize fat models.
• Pseudoinverse. When $A$ is not invertible — rectangular or rank-deficient — the Moore–Penrose pseudoinverse $A^+ = V\Sigma^+ U^T$ gives you the "best available" solution to $A\mathbf{x}=\mathbf{b}$. $\Sigma^+$ is $\Sigma$ with non-zero singular values inverted. This is what numpy.linalg.lstsq is doing under the hood.
• Condition number. $\kappa(A) = \sigma_1/\sigma_n$. Tells you how numerically dangerous a linear system is. Close to 1: safe. Huge: forget it. This is the single most useful diagnostic in numerical linear algebra.

The Eckart–Young truncation:

17. Tensors

You've seen scalars, vectors, and matrices. A tensor is just the natural generalisation: an array of numbers indexed by more than two indices, equipped with rules for how it transforms when you change coordinates. Vectors are rank-1 tensors. Matrices are rank-2 tensors. A 3-D grid of numbers is a rank-3 tensor, and so on.

Einstein summation notation

Writing $\sum_j A_{ij} v_j$ for every matrix-vector product gets tedious. Einstein's convention: if an index appears twice in a term, it is summed over. You never write the $\Sigma$.

NumPy and PyTorch both implement this with einsum: torch.einsum('bik,bkj->bij', A, B) is a batched matrix multiply. Reading einsum strings is a core skill in ML engineering.

Tensor operations

Outer product

$u \otimes v = u_i v_j$ — takes two rank-1 tensors and produces a rank-2 tensor (a matrix). Generalises: $u \otimes v \otimes w = u_i v_j w_k$ gives rank-3. The tensor product of spaces is built from this.

Contraction

Sum over a pair of matching indices. Matrix-vector product is a contraction. The trace is a contraction of both indices of a matrix. Each contraction reduces the rank by 2.

Transpose / permute

Reorder the indices: $A_{ij} \to A_{ji}$ is the usual transpose. For higher-rank tensors, a permutation like $(0, 2, 1, 3)$ reorders the four axes. In PyTorch: x.permute(0, 2, 1, 3).

Reshape / view

Change the shape without touching the data. A $(B, H \times W, C)$ tensor can be reshaped to $(B, H, W, C)$ — same underlying array, different index interpretation. Reshape does not change values.

Frobenius norm

$\|T\|_F = \sqrt{\sum_{i_1 \cdots i_k} T_{i_1 \cdots i_k}^2}$ — the generalisation of the matrix Frobenius norm to any rank. Compute it in NumPy with np.linalg.norm(T).

Tensors in ML practice

Modern deep learning frameworks (PyTorch, JAX, TensorFlow) are built around tensors. Here is what each axis typically means:

The multi-head attention formula $\text{softmax}\!\left(\frac{QK^T}{\sqrt{d_k}}\right)V$ is really a rank-4 einsum: 'bhid,bhjd->bhij' for the dot-product step, followed by 'bhij,bhjd->bhid' for the weighted sum — all batched over both the batch dimension $b$ and the heads dimension $h$ simultaneously.

The connection to matrices

Any rank-$k$ tensor with shape $(n_1, n_2, \ldots, n_k)$ can be unfolded (also called matricised) into a matrix by grouping some indices as rows and the rest as columns. This is the key insight behind tensor decompositions:

• Tucker decomposition — generalises SVD to higher rank tensors: $T_{ijk} = \sum_{pqr} G_{pqr} U_{ip} V_{jq} W_{kr}$, where $G$ is a small core tensor and $U, V, W$ are factor matrices.
• CP decomposition — expresses a tensor as a sum of outer products (rank-1 tensors): $T \approx \sum_{r=1}^R \lambda_r\, a_r \otimes b_r \otimes c_r$. Analogue of SVD where $k = 1$ for every term.
• Tensor train / matrix product state — chain of small matrices, used in physics and to compress huge embedding tables.

// ── Tensor fundamentals ──────────────────────────────────────
// Outer product (rank-1 ⊗ rank-1 → rank-2)
function outer(u[n], v[m]) -> T[n][m]:
    for i in 0..n:
        for j in 0..m:
            T[i][j] = u[i] * v[j]
    return T

// Einstein summation: matrix multiply
// 'ik,kj->ij'  means sum over k
function matmul_einsum(A[m][k], B[k][n]) -> C[m][n]:
    for i in 0..m:
        for j in 0..n:
            C[i][j] = sum(A[i][p] * B[p][j]  for p in 0..k)

// Tensor contraction: contract axis-1 of T3 with axis-0 of M
// T3[i][j][k], M[j][l]  →  R[i][l][k]  (sum over j)
function contract(T3[I][J][K], M[J][L]) -> R[I][L][K]:
    for i,l,k:
        R[i][l][k] = sum(T3[i][j][k] * M[j][l]  for j in 0..J)

// Reshape: (A*B, C) → (A, B, C)  -- no data copy
function reshape(T[N][C], A, B):
    assert A * B == N
    return view T as T2[A][B][C]

// Frobenius norm over arbitrary shape
function frobenius(T):
    return sqrt(sum(x*x for x in flatten(T)))

import numpy as np

# ── Outer product ─────────────────────────────────────────────
u = np.array([2.0, 1.0, -1.0])
v = np.array([1.0, 3.0, 2.0])

outer = np.outer(u, v)          # shape (3, 3)
# Equivalently: np.einsum('i,j->ij', u, v)
print("u ⊗ v:\n", outer)

# ── Einstein summation ────────────────────────────────────────
A = np.random.randn(4, 5)
B = np.random.randn(5, 6)

# Matrix multiply
C = np.einsum('ik,kj->ij', A, B)            # (4, 6)
assert np.allclose(C, A @ B)

# Dot product
a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
dot = np.einsum('i,i->', a, b)              # scalar
print("dot:", dot)

# Trace
M = np.eye(4) * 3
tr = np.einsum('ii->', M)                   # 12.0

# Batch matmul (N, I, K) × (N, K, J) → (N, I, J)
N = 8
X = np.random.randn(N, 4, 5)
Y = np.random.randn(N, 5, 6)
Z = np.einsum('bik,bkj->bij', X, Y)        # (8, 4, 6)
assert np.allclose(Z, X @ Y)

# ── Rank-3 tensor operations ──────────────────────────────────
T = np.random.randn(2, 3, 4)   # shape (2,3,4)

# Permute axes
T_perm = np.transpose(T, (1, 0, 2))   # (3, 2, 4)

# Reshape / view
T_flat = T.reshape(2, 12)              # (2, 12)
T_back = T_flat.reshape(2, 3, 4)      # back to (2,3,4)

# Frobenius norm
fnorm = np.linalg.norm(T)

# ── CNN-style 4-D tensor (N, C, H, W) ────────────────────────
batch = np.random.randn(16, 3, 32, 32)  # 16 RGB 32×32 images
print("Batch shape:", batch.shape)

# Sum over spatial dims → (N, C)
spatial_mean = batch.mean(axis=(2, 3))
print("Spatial mean shape:", spatial_mean.shape)

# ── Tucker-like low-rank outer product sum ────────────────────
# Rank-1 approximation: T ≈ a ⊗ b ⊗ c
a_ = np.array([1.0, -1.0])
b_ = np.array([0.5, 1.0, -0.5])
c_ = np.array([2.0, 0.0, -1.0, 1.0])
rank1 = np.einsum('i,j,k->ijk', a_, b_, c_)   # (2,3,4)

// ── Outer product ─────────────────────────────────────────────
function outer(u, v) {
  return u.map(ui => v.map(vj => ui * vj));
}

// ── Einsum helpers (basic cases) ─────────────────────────────
// Matrix multiply: 'ik,kj->ij'
function matmul(A, B) {
  const [m, k] = [A.length, A[0].length];
  const n = B[0].length;
  return Array.from({length: m}, (_, i) =>
    Array.from({length: n}, (_, j) =>
      Array.from({length: k}, (_, p) => A[i][p] * B[p][j])
        .reduce((s, x) => s + x, 0)
    )
  );
}

// Dot product: 'i,i->'
function dot(a, b) {
  return a.reduce((s, ai, i) => s + ai * b[i], 0);
}

// Batch matmul: 'bik,bkj->bij'
function batchMatmul(A, B) {
  return A.map((a, b) => matmul(a, B[b]));
}

// ── Rank-3 tensor indexing ────────────────────────────────────
// Shape (D0, D1, D2) stored as nested arrays
function zeros(d0, d1, d2) {
  return Array.from({length: d0}, () =>
    Array.from({length: d1}, () => new Float32Array(d2))
  );
}

// Outer product of three vectors → rank-3 tensor
function outer3(a, b, c) {
  return a.map(ai =>
    b.map(bi =>
      c.map(ci => ai * bi * ci)
    )
  );
}

// Frobenius norm
function frobeniusNorm(T) {
  let sum = 0;
  function recurse(x) {
    if (Array.isArray(x)) x.forEach(recurse);
    else sum += x * x;
  }
  recurse(T);
  return Math.sqrt(sum);
}

// ── Demo ──────────────────────────────────────────────────────
const u = [2, 1, -1], v = [1, 3, 2];
console.log("u ⊗ v:", outer(u, v));

const A = [[1,2],[3,4]], B = [[5,6],[7,8]];
console.log("A @ B:", matmul(A, B));

const t3 = outer3([1,-1], [0.5,1,-0.5], [2,0,-1,1]);
console.log("rank-3 shape:", t3.length, t3[0].length, t3[0][0].length);
console.log("Frobenius norm:", frobeniusNorm(t3).toFixed(4));

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/* ── Outer product ────────────────────────────────────────────
   u[m], v[n]  →  out[m*n]  (row-major) */
void outer(const double *u, int m,
           const double *v, int n,
           double *out) {
    for (int i = 0; i < m; i++)
        for (int j = 0; j < n; j++)
            out[i*n + j] = u[i] * v[j];
}

/* ── Matrix multiply: C[m][n] = A[m][k] * B[k][n] ─────────── */
void matmul(const double *A, const double *B, double *C,
            int m, int k, int n) {
    memset(C, 0, m * n * sizeof(double));
    for (int i = 0; i < m; i++)
        for (int p = 0; p < k; p++)
            for (int j = 0; j < n; j++)
                C[i*n + j] += A[i*k + p] * B[p*n + j];
}

/* ── Rank-3 tensor contraction ────────────────────────────────
   T[I][J][K], M[J][L]  →  R[I][L][K]  (sum over j) */
void contract3(const double *T, int I, int J, int K,
               const double *M, int L,
               double *R) {
    memset(R, 0, I * L * K * sizeof(double));
    for (int i = 0; i < I; i++)
      for (int j = 0; j < J; j++)
        for (int k = 0; k < K; k++)
          for (int l = 0; l < L; l++)
            R[(i*L + l)*K + k] +=
              T[(i*J + j)*K + k] * M[j*L + l];
}

/* ── Frobenius norm ───────────────────────────────────────────*/
double frob(const double *T, int n) {
    double s = 0;
    for (int i = 0; i < n; i++) s += T[i]*T[i];
    return sqrt(s);
}

int main(void) {
    double u[] = {2, 1, -1}, v[] = {1, 3, 2};
    double M[9];
    outer(u, 3, v, 3, M);
    printf("u⊗v[0][1] = %.1f  (expect 6)\n", M[0*3+1]);

    double A[] = {1,2,3,4}, B[] = {5,6,7,8}, C[4];
    matmul(A, B, C, 2, 2, 2);
    printf("(A@B)[0][0] = %.0f  (expect 19)\n", C[0]);
    printf("Frobenius |M| = %.4f\n", frob(M, 9));
    return 0;
}

#include <iostream>
#include <vector>
#include <cmath>
#include <numeric>
#include <cassert>

// ── Simple rank-N tensor (flat storage + shape) ───────────────
struct Tensor {
    std::vector<int> shape;
    std::vector<double> data;

    Tensor(std::vector<int> sh) : shape(sh) {
        int n = 1;
        for (int d : sh) n *= d;
        data.assign(n, 0.0);
    }

    // Linear index from multi-index
    int idx(std::vector<int> ix) const {
        int lin = 0, stride = 1;
        for (int d = shape.size()-1; d >= 0; --d) {
            lin += ix[d] * stride;
            stride *= shape[d];
        }
        return lin;
    }

    double& at(std::vector<int> ix) { return data[idx(ix)]; }
    double  at(std::vector<int> ix) const { return data[idx(ix)]; }

    double frobenius() const {
        double s = 0;
        for (double x : data) s += x*x;
        return std::sqrt(s);
    }
};

// ── Outer product: rank-1 ⊗ rank-1 → rank-2 ──────────────────
Tensor outer(const std::vector<double>& u,
             const std::vector<double>& v) {
    int m = u.size(), n = v.size();
    Tensor T({m, n});
    for (int i = 0; i < m; i++)
        for (int j = 0; j < n; j++)
            T.at({i,j}) = u[i] * v[j];
    return T;
}

// ── Matrix multiply via einsum 'ik,kj->ij' ───────────────────
Tensor matmul(const Tensor& A, const Tensor& B) {
    int m = A.shape[0], k = A.shape[1], n = B.shape[1];
    assert(B.shape[0] == k);
    Tensor C({m, n});
    for (int i = 0; i < m; i++)
      for (int p = 0; p < k; p++)
        for (int j = 0; j < n; j++)
            C.at({i,j}) += A.at({i,p}) * B.at({p,j});
    return C;
}

int main() {
    std::vector<double> u = {2, 1, -1}, v = {1, 3, 2};
    auto M = outer(u, v);
    std::cout << "M[0][1] = " << M.at({0,1}) << "  (expect 6)\n";
    std::cout << "Frobenius = " << M.frobenius() << "\n";

    Tensor A({2,2}), B({2,2});
    A.at({0,0})=1; A.at({0,1})=2; A.at({1,0})=3; A.at({1,1})=4;
    B.at({0,0})=5; B.at({0,1})=6; B.at({1,0})=7; B.at({1,1})=8;
    auto C = matmul(A, B);
    std::cout << "(A@B)[0][0] = " << C.at({0,0}) << "  (expect 19)\n";
    return 0;
}

import java.util.Arrays;

public class Tensors {

    // ── Outer product u[m] ⊗ v[n] → M[m][n] ─────────────────
    static double[][] outer(double[] u, double[] v) {
        double[][] M = new double[u.length][v.length];
        for (int i = 0; i < u.length; i++)
            for (int j = 0; j < v.length; j++)
                M[i][j] = u[i] * v[j];
        return M;
    }

    // ── Matrix multiply A[m][k] @ B[k][n] → C[m][n] ─────────
    static double[][] matmul(double[][] A, double[][] B) {
        int m = A.length, k = A[0].length, n = B[0].length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++)
          for (int p = 0; p < k; p++)
            for (int j = 0; j < n; j++)
                C[i][j] += A[i][p] * B[p][j];
        return C;
    }

    // ── Frobenius norm ────────────────────────────────────────
    static double frobenius(double[][] M) {
        double s = 0;
        for (double[] row : M)
            for (double x : row) s += x * x;
        return Math.sqrt(s);
    }

    // ── Rank-3 outer product u⊗v⊗w ───────────────────────────
    static double[][][] outer3(double[] u, double[] v, double[] w) {
        double[][][] T = new double[u.length][v.length][w.length];
        for (int i = 0; i < u.length; i++)
          for (int j = 0; j < v.length; j++)
            for (int k = 0; k < w.length; k++)
                T[i][j][k] = u[i] * v[j] * w[k];
        return T;
    }

    public static void main(String[] args) {
        double[] u = {2, 1, -1}, v = {1, 3, 2};
        double[][] M = outer(u, v);
        System.out.printf("M[0][1] = %.1f  (expect 6.0)%n", M[0][1]);
        System.out.printf("Frobenius = %.4f%n", frobenius(M));

        double[][] A = {{1,2},{3,4}}, B = {{5,6},{7,8}};
        double[][] C = matmul(A, B);
        System.out.printf("(A@B)[0][0] = %.0f  (expect 19)%n", C[0][0]);

        double[] w = {0.5, -1, 2};
        double[][][] T3 = outer3(u, v, w);
        System.out.printf("T3[0][1][2] = %.1f  (expect %.1f)%n",
            T3[0][1][2], u[0]*v[1]*w[2]);
    }
}

package main

import (
	"fmt"
	"math"
)

// ── Outer product ─────────────────────────────────────────────
func outer(u, v []float64) [][]float64 {
	m, n := len(u), len(v)
	M := make([][]float64, m)
	for i := range M {
		M[i] = make([]float64, n)
		for j := range M[i] {
			M[i][j] = u[i] * v[j]
		}
	}
	return M
}

// ── Matrix multiply ───────────────────────────────────────────
func matmul(A, B [][]float64) [][]float64 {
	m, k, n := len(A), len(A[0]), len(B[0])
	C := make([][]float64, m)
	for i := range C {
		C[i] = make([]float64, n)
	}
	for i := 0; i < m; i++ {
		for p := 0; p < k; p++ {
			for j := 0; j < n; j++ {
				C[i][j] += A[i][p] * B[p][j]
			}
		}
	}
	return C
}

// ── Frobenius norm ────────────────────────────────────────────
func frobenius(M [][]float64) float64 {
	var s float64
	for _, row := range M {
		for _, x := range row {
			s += x * x
		}
	}
	return math.Sqrt(s)
}

// ── Rank-3 outer product ──────────────────────────────────────
func outer3(a, b, c []float64) [][][]float64 {
	T := make([][][]float64, len(a))
	for i := range T {
		T[i] = make([][]float64, len(b))
		for j := range T[i] {
			T[i][j] = make([]float64, len(c))
			for k := range T[i][j] {
				T[i][j][k] = a[i] * b[j] * c[k]
			}
		}
	}
	return T
}

func main() {
	u := []float64{2, 1, -1}
	v := []float64{1, 3, 2}
	M := outer(u, v)
	fmt.Printf("M[0][1] = %.1f  (expect 6.0)\n", M[0][1])
	fmt.Printf("Frobenius = %.4f\n", frobenius(M))

	A := [][]float64{{1, 2}, {3, 4}}
	B := [][]float64{{5, 6}, {7, 8}}
	C := matmul(A, B)
	fmt.Printf("(A@B)[0][0] = %.0f  (expect 19)\n", C[0][0])

	w := []float64{0.5, -1, 2}
	T3 := outer3(u, v, w)
	fmt.Printf("T3[0][1][2] = %.1f  (expect %.1f)\n",
		T3[0][1][2], u[0]*v[1]*w[2])
}

18. PCA in 30 lines

Let's walk through PCA (Principal Component Analysis) as a single, concrete application of SVD — no more abstract fanfare.

You have $N$ data points, each a $d$-dimensional vector. You stack them into a data matrix $X$ of shape $N \times d$. You want to represent each point with just $k \ll d$ numbers, while losing as little information as possible.

The recipe:

1. Center the data. Subtract the mean of each column. Call the result $\tilde X$. This ensures the data cloud is centered at the origin.
2. Compute the SVD: $\tilde X = U\Sigma V^T$.
3. The columns of $V$ are the principal components. $V_1$ is the direction of maximum variance in the data; $V_2$ is the direction of second-most variance, perpendicular to $V_1$; and so on.
4. Project: the $k$-dimensional representation of each data point is $\tilde X V_{:,1:k}$, where you keep only the first $k$ columns of $V$.
5. Reconstruct (optional): $\tilde X_{\text{approx}} = \tilde X V_{:,1:k} V_{:,1:k}^T$.

The variance captured by the $i$-th principal component is $\sigma_i^2 / (N-1)$. Sum the top $k$ squared singular values and divide by the total sum of squared singular values; that's the fraction of variance you kept. If it's 99% with $k=20$ on 1000-dimensional data, you just compressed by 50× with almost no loss.

19. Applications to ML

Here's the quick survey. Every item on this list reuses the stuff above.

Embeddings are vectors

Every embedding model maps things (words, sentences, images, users, molecules) to vectors in $\mathbb{R}^d$, typically with $d$ between 128 and 4096. You then compare them with cosine similarity — which we showed is dot product normalized by norms. The entire business of semantic search is linear algebra over these embeddings.

A neural network layer is $W\mathbf{x}+\mathbf{b}$

The fully connected ("dense") layer, the workhorse of every deep network, does:

Attention is $\text{softmax}(QK^T/\sqrt{d_k})V$

Self-attention computes, for each query vector, a weighted sum of value vectors, where the weights are softmaxed dot products between the query and keys. Three matrix multiplies plus a scaling. See self-attention for the full mechanics — the math is entirely the stuff on this page.

Softmax, layer norm, residuals

• Softmax is coordinate-wise exponentiation then division by a sum. It turns a vector of "scores" into a probability distribution. Not linear, but it operates on vectors from linear-algebraic pipelines.
• Layer normalization rescales each vector to unit variance and zero mean — a norm-based operation.
• Residual connections add the input back to the output: $\mathbf{y} = \mathbf{x} + f(\mathbf{x})$. This is just vector addition, but it changes the optimization landscape dramatically.

Weight matrices, and why they're often nearly low-rank

Empirically, the weight matrices of trained transformers are approximately low-rank — their spectrum decays quickly. That's the observation LoRA exploits: instead of fine-tuning the full $W$, add a low-rank update $\Delta W = BA$ where $A \in \mathbb{R}^{r \times d}$ and $B \in \mathbb{R}^{d \times r}$ with $r$ small (8 to 64). You update $r \cdot 2d$ parameters instead of $d^2$. Same SVD intuition, applied to updates instead of weights.

20. NumPy code

Three tabs: core matrix operations, eigendecomposition, and SVD + PCA. No library beyond NumPy.

import numpy as np

# Vectors
a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])

# Dot product — four equivalent forms
print(np.dot(a, b))          # 32.0
print(a @ b)                 # 32.0  (Python 3.5+ operator)
print((a * b).sum())         # 32.0
print(np.einsum("i,i->", a, b))# 32.0

# Norms
print(np.linalg.norm(a))                  # L2 norm ≈ 3.742
print(np.linalg.norm(a, ord=1))           # L1 norm = 6.0
print(np.linalg.norm(a, ord=np.inf))      # L-inf norm = 3.0

# Matrices
A = np.array([[1.0, 2.0], [3.0, 4.0]])
B = np.array([[5.0, 6.0], [7.0, 8.0]])

# Matrix multiplication — use @, not *
C = A @ B                             # [[19, 22], [43, 50]]
print(C)

# A * B is ELEMENT-WISE, not matrix multiplication. Classic bug.
elementwise = A * B                   # [[ 5, 12], [21, 32]]

# Transpose, inverse, determinant
print(A.T)                          # [[1, 3], [2, 4]]
print(np.linalg.inv(A))              # [[-2, 1], [1.5, -0.5]]
print(np.linalg.det(A))              # -2.0

# Solving Ax = b — NEVER compute inv(A) @ b. Use solve().
b_vec = np.array([5.0, 11.0])
x = np.linalg.solve(A, b_vec)        # [1., 2.] (unique solution exists)
print("check residual:", A @ x - b_vec)

import numpy as np

# A general 2x2 matrix
A = np.array([[4.0, 1.0],
              [2.0, 3.0]])

# Eigenvalues and eigenvectors (columns of V are the eigenvectors)
vals, V = np.linalg.eig(A)
print("eigenvalues:", vals)            # [5., 2.]
print("eigenvectors (cols):", V)

# Sanity check: A v = lambda v for each eigenpair
for i in range(len(vals)):
    lhs = A @ V[:, i]
    rhs = vals[i] * V[:, i]
    print(np.allclose(lhs, rhs))       # True

# Reconstruct A = P D P^-1 (diagonalization)
D = np.diag(vals)
A_rebuilt = V @ D @ np.linalg.inv(V)
print(np.allclose(A_rebuilt, A))        # True

# For SYMMETRIC matrices, use eigh — it is faster and gives
# real eigenvalues and orthonormal eigenvectors.
S = np.array([[2.0, 1.0],
              [1.0, 3.0]])
w, Q = np.linalg.eigh(S)
print(w)                              # real, sorted ascending
print(Q.T @ Q)                         # ≈ identity: Q is orthogonal

# Power iteration: find the dominant eigenvalue without eig().
# Used for PageRank, and as a building block for Lanczos / Arnoldi.
def power_iter(M, steps=100):
    v = np.random.randn(M.shape[0])
    for _ in range(steps):
        v = M @ v
        v /= np.linalg.norm(v)
    lam = v @ (M @ v)                 # Rayleigh quotient
    return lam, v

lam, v = power_iter(A)
print("dominant lambda ≈", lam)       # ≈ 5.0

import numpy as np

# =========  SVD of an arbitrary (non-square) matrix  =========
A = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0],
              [7.0, 8.0, 9.0],
              [10., 11., 12.]])   # 4x3, rank 2 (columns are collinear in a 2-plane)

U, s, Vt = np.linalg.svd(A, full_matrices=False)
print("singular values:", s)         # e.g. [25.5, 1.29, ~0]
print("rank ≈", np.sum(s > 1e-10))   # 2

# Reconstruct A from its SVD
A_rebuilt = U @ np.diag(s) @ Vt
print(np.allclose(A_rebuilt, A))      # True

# Rank-1 approximation (Eckart–Young)
A1 = s[0] * np.outer(U[:, 0], Vt[0])
print("rank-1 error:", np.linalg.norm(A - A1))

# =========  PCA in a dozen lines  =========
def pca(X, k):
    # X: (N, d) data matrix
    mu = X.mean(axis=0, keepdims=True)
    Xc = X - mu                       # center the data
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    components = Vt[:k]               # (k, d) — principal axes
    scores     = Xc @ components.T    # (N, k) — compressed coords
    var_explained = (s[:k] ** 2) / (s ** 2).sum()
    return scores, components, var_explained, mu

# Try it on random 100-d data with hidden 5-d structure
rng = np.random.default_rng(0)
latent = rng.standard_normal((1000, 5))
W      = rng.standard_normal((5, 100))
X      = latent @ W + 0.01 * rng.standard_normal((1000, 100))

scores, comps, var_exp, mu = pca(X, k=5)
print("variance captured by top-5:", var_exp.sum())   # >0.999

# Reconstruct
X_hat = scores @ comps + mu
print("rmse:", np.sqrt(((X - X_hat) ** 2).mean()))

21. Summary

• Vectors are lists of numbers with geometric meaning. You add them component-wise and stretch them with scalars. The three most useful norms are $\ell_1$, $\ell_2$, and $\ell_\infty$.
• The dot product is the simultaneously algebraic ($\sum a_i b_i$) and geometric ($\|a\|\|b\|\cos\theta$) measure of alignment. It's what makes attention, embeddings search, and projection work.
• A matrix is a linear function from one vector space to another. Its columns are the images of the basis vectors. Matrix multiplication is function composition, which is why it isn't commutative.
• $A\mathbf{x} = \mathbf{b}$ has either one, zero, or infinitely many solutions, and Gaussian elimination tells you which.
• Rank is "how many independent directions the matrix sees"; nullity is "how many it kills"; rank-nullity says they add up to $n$.
• The determinant is the signed volume scaling factor. Zero determinant = rank-deficient = non-invertible.
• Eigenvectors are directions a matrix preserves (up to scaling); eigenvalues are the scale factors. A matrix is diagonalizable when it has a full set of independent eigenvectors, and in that basis it acts like a diagonal matrix.
• SVD is the ultimate factorization: every matrix is a rotation, a scaling, and another rotation. It gives you PCA, low-rank approximation, the pseudoinverse, and the condition number.
• Tensors are generalisations of vectors and matrices to arbitrary rank. Einstein notation (repeated index = sum) is the shorthand used everywhere in ML code. einsum in NumPy/PyTorch implements it directly.
• Every modern ML technique you've heard of — embeddings, attention, weight matrices, PCA, LoRA — is linear algebra with a few layers of glue on top.

See also