Linear Algebra for Deep Learning: Revision Notes

These are my own study notes on the linear algebra you need before reading deep learning papers. They follow the topics in Chapter 2 of Deep Learning by Goodfellow, Bengio and Courville (online chapter). I wrote everything in my own words and added derivations, small numeric examples, NumPy checks and self-test questions. Read the original chapter too. It is free and short.

Table of contents

1. Objects: scalars, vectors, matrices, tensors
2. Products of matrices and vectors
3. Identity and inverse
4. Linear dependence, span and solving Ax = b
5. Norms
6. Special matrices and vectors
7. Eigendecomposition
8. Singular value decomposition
9. The Moore-Penrose pseudoinverse
10. Trace
11. Determinant
12. Worked example: deriving PCA
13. Cheat sheet
14. Further references

1. Objects: scalars, vectors, matrices, tensors

Things worth remembering:

• Vectors are columns by default. $\mathbf{x} \in \mathbb{R}^n$ is an $n \times 1$ matrix. A row vector is written $\mathbf{x}^\top$.
• Index sets. For $S = {1, 3}$, $\mathbf{x}_S = (x_1, x_3)$. A minus sign means “all except”: $\mathbf{x}_{-1}$ drops $x_1$, and $\mathbf{x}_{-S}$ drops every index in $S$.
• Rows and columns. $\mathbf{A}_{i,:}$ is row $i$. $\mathbf{A}_{:,j}$ is column $j$. For an expression such as $f(\mathbf{A})$, write the entry as $f(\mathbf{A})_{i,j}$.
• Tensors are just arrays with more than two axes. A batch of RGB images is a 4-axis tensor of shape (batch, channels, height, width).

Transpose

Picture it as a mirror image across the main diagonal (top left to bottom right). An $m \times n$ matrix becomes $n \times m$. A scalar is its own transpose: $a^\top = a$. This small fact is used in many proofs below.

Element-wise operations and broadcasting

• Matrices of the same shape add entry by entry: $C_{i,j} = A_{i,j} + B_{i,j}$.
• Scalars act on every entry: $(a\mathbf{B} + c)_{i,j} = a B_{i,j} + c$.
• Broadcasting (a deep learning convention, not standard math): adding a vector to a matrix adds it to every row. \(\mathbf{C} = \mathbf{A} + \mathbf{b} \quad\Longleftrightarrow\quad C_{i,j} = A_{i,j} + b_j\) This is exactly what happens with a bias term: $\mathbf{H} = \mathbf{X}\mathbf{W} + \mathbf{b}$, where each row of $\mathbf{X}$ is one example.

import numpy as np
A = np.zeros((3, 2)); b = np.array([1., 2.])
print(A + b)   # every row becomes [1, 2]

2. Products of matrices and vectors

Matrix product

For $\mathbf{A} \in \mathbb{R}^{m \times n}$ and $\mathbf{B} \in \mathbb{R}^{n \times p}$:

The inner dimensions must match. Entry $(i, j)$ is the dot product of row $i$ of $\mathbf{A}$ with column $j$ of $\mathbf{B}$.

Cost. Computing $\mathbf{AB}$ naively takes $m \cdot n \cdot p$ multiply-adds. This is the number to use when counting FLOPs of a linear layer (often reported as $2mnp$ FLOPs, counting multiply and add separately).

Hadamard product

Element-wise, same shapes. Do not confuse it with the matrix product. In NumPy, A * B is Hadamard and A @ B is the matrix product.

Dot product

It is a matrix product of a $1 \times n$ and an $n \times 1$ matrix, so the result is $1 \times 1$, a scalar.

Algebraic properties

Proof that the dot product commutes, using only the transpose rules:

The first step uses “a scalar equals its transpose”. The second uses the product rule for transposes.

Associativity matters for speed. For $\mathbf{x} \in \mathbb{R}^n$ and $\mathbf{A}, \mathbf{B} \in \mathbb{R}^{n \times n}$, computing $\mathbf{A}(\mathbf{B}\mathbf{x})$ costs $O(n^2)$, while $(\mathbf{AB})\mathbf{x}$ costs $O(n^3)$. Same answer, very different cost. This is why LoRA computes $\mathbf{B}(\mathbf{A}\mathbf{x})$ instead of forming $\mathbf{BA}$ during training.

Linear systems

A system of $m$ linear equations in $n$ unknowns:

Row $i$ reads $\mathbf{A}_{i,:}\mathbf{x} = b_i$, that is $A_{i,1}x_1 + \dots + A_{i,n}x_n = b_i$.

3. Identity and inverse

Identity. $\mathbf{I}_n \in \mathbb{R}^{n \times n}$ has ones on the diagonal and zeros elsewhere. For every $\mathbf{x} \in \mathbb{R}^n$, $\mathbf{I}_n\mathbf{x} = \mathbf{x}$.

Inverse. $\mathbf{A}^{-1}$ is the matrix with

If it exists, then

Practical warning. The inverse is mostly a tool for reasoning. In code, never compute inv(A) @ b. Use np.linalg.solve(A, b) (LU factorization) or a least-squares routine. Forming the inverse explicitly is slower and loses more precision, especially when $\mathbf{A}$ is ill-conditioned.

A = np.array([[3., 1.], [1., 2.]]); b = np.array([9., 8.])
x = np.linalg.solve(A, b)        # preferred
print(x, np.allclose(A @ x, b))  # [2. 3.] True

4. Linear dependence, span and solving Ax = b

How many solutions can there be?

Exactly zero, one, or infinitely many. Never exactly two.

Proof. Suppose $\mathbf{x}$ and $\mathbf{y}$ are two different solutions. For any $\alpha \in \mathbb{R}$, let $\mathbf{z} = \alpha\mathbf{x} + (1 - \alpha)\mathbf{y}$. Then

Each $\alpha$ gives a different point on the line through $\mathbf{x}$ and $\mathbf{y}$, so there are infinitely many solutions.

The column picture

The most useful way to read $\mathbf{A}\mathbf{x}$:

$\mathbf{A}\mathbf{x}$ is a linear combination of the columns of $\mathbf{A}$, with weights given by $\mathbf{x}$. Each $x_j$ says “how far to move along column $j$”.

• The span of a set of vectors is every linear combination of them.
• The column space (or range) of $\mathbf{A}$ is the span of its columns.
• $\mathbf{A}\mathbf{x} = \mathbf{b}$ has a solution if and only if $\mathbf{b}$ is in the column space of $\mathbf{A}$.

When is there a solution for every b?

We need the column space to be all of $\mathbb{R}^m$.

• Necessary: $n \ge m$ (at least as many columns as dimensions of $\mathbf{b}$). With fewer columns, the column space is at most an $n$-dimensional subspace of $\mathbb{R}^m$. Example: a $3 \times 2$ matrix maps into a plane inside $\mathbb{R}^3$, so most $\mathbf{b}$ are unreachable.
• Not sufficient: columns can be redundant. The matrix \(\begin{pmatrix} 1 & 1 \\ 2 & 2 \end{pmatrix}\) has two columns but both point the same way. Its column space is a line, not $\mathbb{R}^2$.

Linear independence

A set of vectors is linearly independent if no vector in it is a linear combination of the others. Equivalently, $\sum_j c_j \mathbf{v}_j = \mathbf{0}$ forces all $c_j = 0$.

So the precise condition is: the columns of $\mathbf{A}$ must contain a set of exactly $m$ linearly independent vectors. (A set of $m$-dimensional vectors cannot have more than $m$ independent members.)

The number of linearly independent columns is the rank of $\mathbf{A}$. Solvable for every $\mathbf{b}$ means $\operatorname{rank}(\mathbf{A}) = m$.

When is the solution unique?

Uniqueness needs at most $m$ columns. Otherwise there are extra directions that can be traded against each other (a nonzero null space), giving infinitely many solutions.

Putting both conditions together: solvable for every $\mathbf{b}$ with a unique solution requires $\mathbf{A}$ to be square ($m = n$) with all columns linearly independent. Such a matrix is invertible.

A square matrix with linearly dependent columns is singular. A singular or non-square $\mathbf{A}$ can still give solutions for some $\mathbf{b}$, but not via matrix inversion.

Left and right inverses. So far we used $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ (left inverse). One can also ask for $\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}$ (right inverse). For square matrices the two are the same matrix.

Summary table

5. Norms

A norm measures the size of a vector. Formally, $f$ is a norm if for all $\mathbf{x}, \mathbf{y}$ and scalars $\alpha$:

1. $f(\mathbf{x}) = 0 \Rightarrow \mathbf{x} = \mathbf{0}$
2. $f(\mathbf{x} + \mathbf{y}) \le f(\mathbf{x}) + f(\mathbf{y})$ (triangle inequality)

3. $f(\alpha\mathbf{x}) =

   \alpha

   f(\mathbf{x})$ (absolute homogeneity)

Non-negativity, $f(\mathbf{x}) \ge 0$, follows from these.

The $L^p$ family

For $p < 1$ the triangle inequality fails, so it is not a norm.

Why squared $L^2$ versus $L^1$?

• Squared $L^2$ is easy to work with: its gradient is $2\mathbf{x}$, and each partial derivative depends only on its own coordinate. But near zero it grows very slowly. A coordinate of $0.01$ contributes $10^{-4}$. So it barely pushes small values all the way to zero.
• $L^1$ grows at the same rate everywhere. Moving a coordinate from $0.01$ to $0$ reduces it by exactly $0.01$, the same as moving from $5.01$ to $5$. So it does push small values to exactly zero. Use it when the difference between “exactly zero” and “small but nonzero” matters.

Gradient view. $\nabla \tfrac{1}{2}|\mathbf{x}|_2^2 = \mathbf{x}$ (shrinks proportionally), while the subgradient of $|\mathbf{x}|_1$ is $\operatorname{sign}(\mathbf{x})$ (constant push toward zero).

The “L0 norm” is not a norm

Counting the nonzero entries is often called the $L^0$ norm. It is not a norm: scaling $\mathbf{x}$ by $\alpha \ne 0$ does not change the count, so homogeneity fails. Because the count is non-differentiable and combinatorial, $L^1$ is usually used as a convex substitute. (Magnitude pruning and SparseGPT impose an $L^0$ budget directly, by choosing which weights to zero.)

Dot product and angle

Normalize both vectors and the dot product is the cosine similarity. This is the score in CLIP, SimCLR and most contrastive losses.

x = np.array([3., -4., 0.])
print(np.linalg.norm(x, 1), np.linalg.norm(x), np.linalg.norm(x, np.inf))
# 7.0 5.0 4.0

6. Special matrices and vectors

Diagonal matrices

$\mathbf{D}$ is diagonal if $D_{i,j} = 0$ whenever $i \ne j$. $\operatorname{diag}(\mathbf{v})$ is the square diagonal matrix with $\mathbf{v}$ on the diagonal.

• Cheap multiplication: $\operatorname{diag}(\mathbf{v})\mathbf{x} = \mathbf{v} \odot \mathbf{x}$. That is $O(n)$, not $O(n^2)$.
• Cheap inverse: exists iff every $v_i \ne 0$, and then $\operatorname{diag}(\mathbf{v})^{-1} = \operatorname{diag}\big((1/v_1, \dots, 1/v_n)^\top\big)$.
• Restricting a model to diagonal matrices is a common trick to reduce cost (diagonal covariance in VAEs, Adam’s per-parameter scaling, LayerNorm’s scale vector).

Rectangular diagonal matrices also exist. For $\mathbf{D} \in \mathbb{R}^{m \times n}$:

• if $m > n$ (tall), $\mathbf{Dx}$ scales each entry and then appends zeros;
• if $m < n$ (wide), $\mathbf{Dx}$ scales the first $m$ entries and drops the rest.

They have no inverse but still multiply cheaply. They appear in the SVD.

Symmetric matrices

$\mathbf{A} = \mathbf{A}^\top$. They arise when entries come from a function that does not care about argument order, for example a distance matrix $A_{i,j} = d(\mathbf{x}^{(i)}, \mathbf{x}^{(j)})$, a covariance matrix, a Gram matrix $\mathbf{X}^\top\mathbf{X}$ or a Hessian.

Unit, orthogonal and orthonormal vectors

• Unit vector: $|\mathbf{x}|_2 = 1$.
• Orthogonal: $\mathbf{x}^\top\mathbf{y} = 0$. If both are nonzero, they are at 90 degrees. In $\mathbb{R}^n$ at most $n$ nonzero vectors can be mutually orthogonal.
• Orthonormal: orthogonal and unit norm.

Orthogonal matrices

A square matrix whose rows are mutually orthonormal and whose columns are mutually orthonormal:

The inverse costs nothing: just transpose.

Facts worth knowing (beyond the chapter):

• Orthogonal matrices preserve lengths and angles: $|\mathbf{Q}\mathbf{x}|_2^2 = \mathbf{x}^\top\mathbf{Q}^\top\mathbf{Q}\mathbf{x} = |\mathbf{x}|_2^2$. They are rotations and reflections.
• $\det(\mathbf{Q}) = \pm 1$.
• Naming trap: “orthogonal matrix” means orthonormal columns. A matrix whose columns are orthogonal but not unit length has no special name.
• Non-square case: a tall $\mathbf{U} \in \mathbb{R}^{n \times l}$ with orthonormal columns satisfies $\mathbf{U}^\top\mathbf{U} = \mathbf{I}_l$ but not $\mathbf{U}\mathbf{U}^\top = \mathbf{I}_n$. Instead $\mathbf{U}\mathbf{U}^\top$ is the projection onto the span of the columns. This is exactly the PCA decoder matrix in Section 12.

7. Eigendecomposition

Why decompose? Writing 12 as $2 \times 2 \times 3$ tells you things the decimal representation hides (divisible by 3, not by 5). Decomposing a matrix similarly reveals properties that do not depend on how we happen to write it down.

Eigenvectors and eigenvalues

A nonzero vector $\mathbf{v}$ is an eigenvector of square $\mathbf{A}$ with eigenvalue $\lambda$ if

$\mathbf{A}$ only stretches $\mathbf{v}$; it does not change its direction (a negative $\lambda$ flips it).

• Left eigenvectors satisfy $\mathbf{v}^\top\mathbf{A} = \lambda\mathbf{v}^\top$. “Eigenvector” alone means right eigenvector.
• If $\mathbf{v}$ is an eigenvector, so is $s\mathbf{v}$ for any $s \ne 0$, with the same $\lambda$. Convention: use unit eigenvectors.
• Finding them (not in the chapter, but asked in interviews): $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$ with $\mathbf{v} \ne \mathbf{0}$ means $(\mathbf{A} - \lambda\mathbf{I})$ is singular, so $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$. This polynomial in $\lambda$ is the characteristic polynomial.

Worked example.

For $\lambda = 3$: $\mathbf{v} = \tfrac{1}{\sqrt 2}(1, 1)^\top$. For $\lambda = 1$: $\mathbf{v} = \tfrac{1}{\sqrt 2}(1, -1)^\top$. They are orthogonal, as expected for a symmetric matrix.

The decomposition

If $\mathbf{A} \in \mathbb{R}^{n \times n}$ has $n$ linearly independent eigenvectors $\mathbf{v}^{(1)}, \dots, \mathbf{v}^{(n)}$, stack them as columns of $\mathbf{V}$ and the eigenvalues in $\boldsymbol{\lambda}$. Then

Why: $\mathbf{A}\mathbf{V} = \mathbf{V}\operatorname{diag}(\boldsymbol{\lambda})$ is just the $n$ eigenvector equations side by side. Multiply both sides on the right by $\mathbf{V}^{-1}$.

Caveats:

• Not every matrix has one. Example: $\begin{pmatrix} 0 & 1 \ 0 & 0 \end{pmatrix}$ has only eigenvalue $0$ and only one independent eigenvector.
• Real matrices can have complex eigenvalues. Example: a 90 degree rotation $\begin{pmatrix} 0 & -1 \ 1 & 0 \end{pmatrix}$ has $\lambda = \pm i$. No real direction stays fixed under a rotation.

Real symmetric matrices: the spectral theorem

Every real symmetric matrix can be decomposed using only real numbers and an orthogonal matrix:

$\mathbf{Q}$ is orthogonal (columns are orthonormal eigenvectors) and $\boldsymbol{\Lambda}$ is diagonal (eigenvalues). Column $i$ of $\mathbf{Q}$ goes with $\Lambda_{i,i}$.

Geometric picture. Write $\mathbf{A}\mathbf{x} = \mathbf{Q}\boldsymbol{\Lambda}(\mathbf{Q}^\top\mathbf{x})$:

1. $\mathbf{Q}^\top\mathbf{x}$: express $\mathbf{x}$ in the eigenvector basis.
2. $\boldsymbol{\Lambda}$: scale coordinate $i$ by $\lambda_i$.
3. $\mathbf{Q}$: go back to the original basis.

Equivalently, $\mathbf{A} = \sum_i \lambda_i \mathbf{q}^{(i)}\mathbf{q}^{(i)\top}$: a weighted sum of rank-one projections.

Uniqueness. If two or more eigenvectors share an eigenvalue, any orthonormal basis of their span works too, so $\mathbf{Q}$ is not unique. Convention: sort eigenvalues in descending order. Under that convention, the decomposition is unique if all eigenvalues are distinct (up to the sign of each eigenvector).

What eigenvalues tell you

• Singular iff some eigenvalue is zero. (A zero eigenvalue means $\mathbf{A}\mathbf{v} = \mathbf{0}$ for some $\mathbf{v} \ne \mathbf{0}$, so information is lost.)
• Inverse of a symmetric invertible matrix: $\mathbf{A}^{-1} = \mathbf{Q}\boldsymbol{\Lambda}^{-1}\mathbf{Q}^\top$. Same eigenvectors, reciprocal eigenvalues.
• Powers: $\mathbf{A}^k = \mathbf{Q}\boldsymbol{\Lambda}^k\mathbf{Q}^\top$.

Optimizing a quadratic form on the unit sphere

For real symmetric $\mathbf{A}$, consider

Result: the maximum of $f$ is $\lambda_{\max}$ and the minimum is $\lambda_{\min}$, attained at the corresponding eigenvectors. If $\mathbf{x}$ is any eigenvector, $f(\mathbf{x})$ equals its eigenvalue.

Proof (derive this yourself). Let $\mathbf{y} = \mathbf{Q}^\top\mathbf{x}$. Since $\mathbf{Q}$ is orthogonal, $|\mathbf{y}|_2 = 1$. Then

The weights $y_i^2$ are non-negative and sum to 1, so $f$ is a weighted average of the eigenvalues. A weighted average lies between the smallest and largest value, with equality when all weight sits on one eigenvalue, that is when $\mathbf{y}$ is a standard basis vector, that is when $\mathbf{x}$ is an eigenvector. $\blacksquare$

The ratio $\mathbf{x}^\top\mathbf{A}\mathbf{x} / \mathbf{x}^\top\mathbf{x}$ is the Rayleigh quotient. This result is the core of PCA (Section 12) and of curvature arguments in optimization: for a Hessian $\mathbf{H}$, $\mathbf{d}^\top\mathbf{H}\mathbf{d}$ is the curvature along unit direction $\mathbf{d}$, and the extremes are the extreme eigenvalues.

Definiteness

The eigenvalue and quadratic-form columns are equivalent, by the weighted-average argument above.

Where this matters:

• Any Gram matrix $\mathbf{X}^\top\mathbf{X}$ is PSD, because $\mathbf{v}^\top\mathbf{X}^\top\mathbf{X}\mathbf{v} = |\mathbf{X}\mathbf{v}|_2^2 \ge 0$. So covariance matrices are PSD.
• A Hessian that is PD at a critical point means a strict local minimum. Indefinite means a saddle point.
• $\mathbf{X}^\top\mathbf{X} + \alpha\mathbf{I}$ with $\alpha > 0$ is PD (eigenvalues shift up by $\alpha$), hence invertible. This is why ridge regression always has a unique solution.

8. Singular value decomposition (SVD)

The eigendecomposition only works for some square matrices. The SVD works for every real matrix, square or not.

Geometric picture. $\mathbf{A}\mathbf{x} = \mathbf{U}(\mathbf{D}(\mathbf{V}^\top\mathbf{x}))$: rotate or reflect in the input space, scale each axis by $\sigma_i$ (and change dimension), rotate or reflect in the output space. Every linear map is “rotate, stretch, rotate”.

Connection to eigendecomposition

Both are eigendecompositions of symmetric PSD matrices. Read off:

• right-singular vectors ($\mathbf{V}$) are eigenvectors of $\mathbf{A}^\top\mathbf{A}$;
• left-singular vectors ($\mathbf{U}$) are eigenvectors of $\mathbf{A}\mathbf{A}^\top$;
• the nonzero singular values are the square roots of the nonzero eigenvalues of $\mathbf{A}^\top\mathbf{A}$ (equivalently of $\mathbf{A}\mathbf{A}^\top$): $\sigma_i = \sqrt{\lambda_i}$.

Useful extras (beyond the chapter)

• Rank = number of nonzero singular values.
• Rank-one expansion: $\mathbf{A} = \sum_{i=1}^{r} \sigma_i \mathbf{u}^{(i)}\mathbf{v}^{(i)\top}$.
• Best low-rank approximation (Eckart-Young): keeping the top $k$ terms gives the rank-$k$ matrix closest to $\mathbf{A}$ in Frobenius (and spectral) norm. The error is $\sqrt{\sum_{i>k}\sigma_i^2}$ in Frobenius norm. This is the reason low-rank adapters (LoRA) and low-rank compression make sense when singular values decay fast.
• $|\mathbf{A}|_F^2 = \sum_i \sigma_i^2$ and the spectral norm $|\mathbf{A}|_2 = \sigma_{\max}$.
• Condition number $\sigma_{\max}/\sigma_{\min}$: large means solving $\mathbf{Ax} = \mathbf{b}$ amplifies errors.

A = np.random.randn(5, 3)
U, s, Vt = np.linalg.svd(A)                     # full: U 5x5, Vt 3x3
eig = np.linalg.eigvalsh(A.T @ A)[::-1]         # descending
print(np.allclose(s**2, eig))                   # True

9. The Moore-Penrose pseudoinverse

Problem. For non-square $\mathbf{A}$ there is no inverse. We would still like a matrix $\mathbf{B}$ so that $\mathbf{x} = \mathbf{B}\mathbf{y}$ solves $\mathbf{A}\mathbf{x} = \mathbf{y}$ as well as possible.

• Tall $\mathbf{A}$ (more rows than columns, more equations than unknowns): usually no exact solution.
• Wide $\mathbf{A}$ (more columns than rows): usually infinitely many solutions.

Definition

Read it as “ridge regression with vanishing regularization”. The $\alpha\mathbf{I}$ makes the matrix invertible (Section 7), and the limit removes the bias.

How it is actually computed

where $\mathbf{D}^+$ is built from $\mathbf{D}$ by taking the reciprocal of each nonzero diagonal entry and then transposing (so $\mathbf{D}^+$ is $n \times m$). Zero singular values stay zero.

Sanity check. If $\mathbf{A}$ is square and invertible, then $\mathbf{D}^+ = \mathbf{D}^{-1}$ and $\mathbf{V}\mathbf{D}^{-1}\mathbf{U}^\top = (\mathbf{U}\mathbf{D}\mathbf{V}^\top)^{-1} = \mathbf{A}^{-1}$.

What it gives you

Least squares link. When $\mathbf{A}$ has full column rank, $\mathbf{A}^+ = (\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top$. Derivation: minimize $|\mathbf{Ax} - \mathbf{y}|_2^2$, set the gradient $2\mathbf{A}^\top(\mathbf{Ax} - \mathbf{y})$ to zero, get the normal equations $\mathbf{A}^\top\mathbf{A}\mathbf{x} = \mathbf{A}^\top\mathbf{y}$.

Minimum norm link. When $\mathbf{A}$ has full row rank, $\mathbf{A}^+ = \mathbf{A}^\top(\mathbf{A}\mathbf{A}^\top)^{-1}$. The solution lies in the row space of $\mathbf{A}$, so it has no component in the null space. Any other solution adds a null-space component, which can only increase the norm (Pythagoras). This is also the classic explanation for why gradient descent from zero on an overparameterized linear regression converges to the minimum-norm interpolator.

A = np.random.randn(2, 4); y = np.random.randn(2)   # wide
x = np.linalg.pinv(A) @ y
print(np.allclose(A @ x, y))                        # True, exact solution

10. Trace

The trace is useful because it lets you write sums and norms as matrix expressions, which are easier to manipulate.

Properties

Cyclic invariance in general:

You may move the last factor to the front (cycle), but you may not swap arbitrary factors: $\operatorname{Tr}(\mathbf{ABC}) \ne \operatorname{Tr}(\mathbf{ACB})$ in general.

It holds even when shapes change. For $\mathbf{A} \in \mathbb{R}^{m \times n}$, $\mathbf{B} \in \mathbb{R}^{n \times m}$: $\mathbf{AB}$ is $m \times m$ and $\mathbf{BA}$ is $n \times n$, yet

The double sum is the proof: both traces sum the same products.

The trace trick (very common): a scalar quadratic form can be rewritten as a trace and cycled.

This is how you take expectations of quadratic forms, $\mathbb{E}[\mathbf{x}^\top\mathbf{A}\mathbf{x}] = \operatorname{Tr}(\mathbf{A}\,\mathbb{E}[\mathbf{x}\mathbf{x}^\top])$, and how you derive the Gaussian MLE for the covariance.

11. Determinant

$\det(\mathbf{A})$ maps a square matrix to a real number.

• $\det = 0$: at least one direction is squashed flat, all volume is lost, and the map cannot be undone (singular).

• $

  \det

  = 1$: volume is preserved (rotations, reflections, shears).

• The sign tells you whether orientation is flipped.

Extras worth knowing:

• $\det(\mathbf{AB}) = \det(\mathbf{A})\det(\mathbf{B})$, $\det(\mathbf{A}^\top) = \det(\mathbf{A})$, $\det(\mathbf{A}^{-1}) = 1/\det(\mathbf{A})$.
• $2 \times 2$: $\det\begin{pmatrix} a & b \ c & d\end{pmatrix} = ad - bc$.

• The determinant appears in the Gaussian density normalizer $

  2\pi\boldsymbol{\Sigma}

  ^{-1/2}$ and in the change-of-variables formula behind normalizing flows, $\log p_X(\mathbf{x}) = \log p_Z(f(\mathbf{x})) + \log

  \det \mathbf{J}_f(\mathbf{x})

  $.

• Use np.linalg.slogdet, not det, for anything large: determinants overflow or underflow quickly.

12. Worked example: deriving PCA

This section uses nearly everything above. Derive it on paper at least twice.

Setup

We have $m$ points ${\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(m)}}$ in $\mathbb{R}^n$. We want lossy compression: store a code $\mathbf{c}^{(i)} \in \mathbb{R}^l$ with $l < n$ for each point.

• Encoder: $f(\mathbf{x}) = \mathbf{c}$.
• Decoder: $g(\mathbf{c}) \approx \mathbf{x}$.
• Reconstruction: $r(\mathbf{x}) = g(f(\mathbf{x}))$.

Design choices

1. Linear decoder: $g(\mathbf{c}) = \mathbf{D}\mathbf{c}$ with $\mathbf{D} \in \mathbb{R}^{n \times l}$.
2. Orthogonal columns of $\mathbf{D}$. (Note: $\mathbf{D}$ is not square, so it is not an “orthogonal matrix”, but $\mathbf{D}^\top\mathbf{D} = \mathbf{I}_l$.)
3. Unit-norm columns. Without this the problem has no unique solution: scale a column of $\mathbf{D}$ by $s$ and the matching code entries by $1/s$, and the reconstruction is unchanged.

So the constraint is $\mathbf{D}^\top\mathbf{D} = \mathbf{I}_l$.

We assume the data is centered (mean subtracted). Without centering, the method finds directions of largest second moment about the origin, not largest variance.

Step 1: the optimal code for a fixed D

Choose the code that minimizes reconstruction error:

Square the norm. $|\cdot|_2 \ge 0$ and $t \mapsto t^2$ is increasing on $[0, \infty)$, so the minimizer is the same:

Expand (distributivity, and $\mathbf{x}^\top g(\mathbf{c}) = g(\mathbf{c})^\top\mathbf{x}$ because a scalar is its own transpose):

Drop $\mathbf{x}^\top\mathbf{x}$, which does not depend on $\mathbf{c}$. Substitute $g(\mathbf{c}) = \mathbf{Dc}$ and use $\mathbf{D}^\top\mathbf{D} = \mathbf{I}_l$:

Set the gradient to zero:

The objective is a convex quadratic in $\mathbf{c}$ (Hessian $2\mathbf{I}$, which is PD), so this stationary point is the global minimum.

So the encoder is a single matrix-vector product, $f(\mathbf{x}) = \mathbf{D}^\top\mathbf{x}$, and

$\mathbf{D}\mathbf{D}^\top$ is the orthogonal projection onto the column space of $\mathbf{D}$. PCA reconstruction is “project onto a subspace”.

Step 2: choose D

$\mathbf{D}$ is shared across all points, so minimize the total error over the dataset, as a Frobenius norm over all reconstruction errors:

The case $l = 1$

Now $\mathbf{D}$ is a single unit vector $\mathbf{d}$. Drop the square root (same argument as before):

$\mathbf{d}^\top\mathbf{x}^{(i)}$ is a scalar, so it can move to the left and be transposed: $\mathbf{d}\mathbf{d}^\top\mathbf{x}^{(i)} = (\mathbf{x}^{(i)\top}\mathbf{d})\mathbf{d}$.

Matrix form. Stack the points as rows of the design matrix $\mathbf{X} \in \mathbb{R}^{m \times n}$, so $\mathbf{X}_{i,:} = \mathbf{x}^{(i)\top}$. Row $i$ of $\mathbf{X}\mathbf{d}\mathbf{d}^\top$ is $(\mathbf{x}^{(i)\top}\mathbf{d})\mathbf{d}^\top$, the transposed reconstruction of point $i$. So

Rewrite the Frobenius norm as a trace and expand. Using $|\mathbf{M}|_F^2 = \operatorname{Tr}(\mathbf{M}^\top\mathbf{M})$:

Simplify.

• Drop $\operatorname{Tr}(\mathbf{X}^\top\mathbf{X})$, which does not depend on $\mathbf{d}$.
• The two middle terms are equal (cycle the trace, or transpose it), giving $-2\operatorname{Tr}(\mathbf{X}^\top\mathbf{X}\mathbf{d}\mathbf{d}^\top)$.
• Cycle the last term: $\operatorname{Tr}(\mathbf{X}^\top\mathbf{X}\mathbf{d}\mathbf{d}^\top\mathbf{d}\mathbf{d}^\top)$. Now apply the constraint $\mathbf{d}^\top\mathbf{d} = 1$ to get $\operatorname{Tr}(\mathbf{X}^\top\mathbf{X}\mathbf{d}\mathbf{d}^\top)$.

The constraint can only be applied after the cycling brings $\mathbf{d}^\top\mathbf{d}$ together. Before that, nothing simplifies. This is the step people get stuck on.

The objective becomes

Minimizing the negative is maximizing the positive, and one more cycle turns the trace into a scalar:

Solve. This is exactly the quadratic form on the unit sphere from Section 7, with the symmetric PSD matrix $\mathbf{X}^\top\mathbf{X}$. The maximum is the largest eigenvalue, attained at its eigenvector. So $\mathbf{d}$ is the eigenvector of $\mathbf{X}^\top\mathbf{X}$ with the largest eigenvalue: the first principal component.

Alternative via Lagrange multipliers (good to show in an interview): $\mathcal{L} = \mathbf{d}^\top\mathbf{X}^\top\mathbf{X}\mathbf{d} - \lambda(\mathbf{d}^\top\mathbf{d} - 1)$. Setting $\nabla_{\mathbf{d}}\mathcal{L} = 0$ gives $\mathbf{X}^\top\mathbf{X}\mathbf{d} = \lambda\mathbf{d}$, so every stationary point is an eigenvector, and the objective equals $\mathbf{d}^\top(\lambda\mathbf{d}) = \lambda$. Pick the largest.

The general case $l > 1$

The optimal $\mathbf{D}$ has as columns the $l$ eigenvectors of $\mathbf{X}^\top\mathbf{X}$ with the $l$ largest eigenvalues. The book leaves the proof as an exercise. A short version:

The same expansion gives $\mathbf{D}^* = \arg\max_{\mathbf{D}} \operatorname{Tr}(\mathbf{D}^\top\mathbf{X}^\top\mathbf{X}\mathbf{D})$ subject to $\mathbf{D}^\top\mathbf{D} = \mathbf{I}_l$. Write $\mathbf{X}^\top\mathbf{X} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top$ and $\mathbf{W} = \mathbf{Q}^\top\mathbf{D}$ (also orthonormal columns). Then

Each $p_k \in [0, 1]$ (it is a squared length of a row of a matrix with orthonormal columns, which is at most 1), and $\sum_k p_k = |\mathbf{W}|_F^2 = l$. Maximizing $\sum_k\lambda_k p_k$ under these constraints puts $p_k = 1$ on the $l$ largest eigenvalues and $0$ elsewhere. That is achieved when the columns of $\mathbf{D}$ span the top-$l$ eigenvectors. (This is Ky Fan’s theorem. An induction proof adds one direction at a time, each orthogonal to the previous ones.)

Interpretation and links

• If $\mathbf{X}$ is centered, $\tfrac{1}{m}\mathbf{X}^\top\mathbf{X}$ is the sample covariance. PCA directions are the directions of largest variance, and the eigenvalue is the variance along that direction. Minimizing reconstruction error and maximizing projected variance are the same problem.
• Via SVD: if $\mathbf{X} = \mathbf{U}\mathbf{S}\mathbf{V}^\top$, the principal components are the top columns of $\mathbf{V}$, and the eigenvalues of $\mathbf{X}^\top\mathbf{X}$ are $\sigma_i^2$. In practice compute the SVD of $\mathbf{X}$ rather than forming $\mathbf{X}^\top\mathbf{X}$, which squares the condition number.
• The leftover error equals the sum of the discarded eigenvalues (Eckart-Young again).
• A linear autoencoder trained with squared error learns the same subspace as PCA (but not necessarily orthonormal columns).

rng = np.random.default_rng(0)
X = rng.normal(size=(500, 5)) @ rng.normal(size=(5, 5))
X -= X.mean(0)                                   # center
evals, evecs = np.linalg.eigh(X.T @ X)
D = evecs[:, ::-1][:, :2]                        # top-2 eigenvectors
R = X @ D @ D.T                                  # reconstruction (rows)
err = np.linalg.norm(X - R, 'fro')**2
print(np.isclose(err, evals[::-1][2:].sum()))    # True: error = discarded eigenvalues

13. Cheat sheet

14. Further references

• Goodfellow, Bengio, Courville. Deep Learning, Chapter 2. deeplearningbook.org
• Petersen and Pedersen. The Matrix Cookbook. A dense list of identities, useful for lookup.
• Shilov. Linear Algebra (Dover). A full treatment if you want proofs.
• Strang. Introduction to Linear Algebra, and his MIT 18.06 lectures, for the four fundamental subspaces.
• Parr and Howard. “The Matrix Calculus You Need for Deep Learning”. The natural next step after these notes.