Linear Algebra

A neural network is a pipeline of linear maps interleaved with nonlinearities. Strip the nonlinearities and what remains, the weight matrices, the activations, the gradients, is pure linear algebra: a composition of matrix products acting on arrays of numbers. This lesson is the working vocabulary, distilled to the objects a network manipulates and the two decompositions that explain what those objects do to space.¹

Objects: scalars, vectors, matrices, tensors

The data structures form a ladder by number of indices.

A vector $x\in {\mathbb{R}}^n$ is a point in $n$-dimensional space; its entries $x_1,\dots ,x_n$ are coordinates. A matrix $A\in {\mathbb{R}}^{m\times n}$ has entry $A_{i,j}$ in row $i$, column $j$. The transpose mirrors across the main diagonal, $(A^{\top }{)}_{i,j}=A_{j,i}$, turning an $m\times n$ matrix into an $n\times m$ one; for a product it reverses order, $(AB{)}^{\top }=B^{\top }{A}^{\top }$.

Broadcasting lets operations on mismatched shapes proceed by implicitly replicating along missing or size-$1$ axes. Adding a bias vector $b\in {\mathbb{R}}^m$ to every column of $A\in {\mathbb{R}}^{m\times n}$ is the broadcast $C_{i,j}=A_{i,j}+b_i$; Goodfellow writes this $C=A+b$, with $b$ added to each row. The rule is mechanical: align shapes from the right, and any axis of size $1$ (or absent) is stretched to match.² The same tensor objects and shape rules are what every framework's array type exposes directly.³

Two products

There are two different multiplications of arrays, and conflating them is the most common shape bug in practice.

A single dense layer is the matrix product $z=Wx+b$ with $W\in {\mathbb{R}}^{m\times n}$; the elementwise activation and the gradient masks of backpropagation are Hadamard products. The dot product of two vectors is the special case $x^{\top }y={\sum }_i{x}_i{y}_i$, a scalar measuring alignment.⁴

For example, for $A\in {\mathbb{R}}^{2\times 3}$ and $B\in {\mathbb{R}}^{3\times 2}$ the shared inner dimension $3$ is summed away, leaving a $2\times 2$ result:

Entry $C_{i,j}$ is the dot product of row $i$ of the left factor with column $j$ of the right factor. When a layer processes a batch of $N$ examples at once, the same rule applies with the batch stacked as rows: $X\in {\mathbb{R}}^{N\times n}$ times $W^{\top }\in {\mathbb{R}}^{n\times m}$ yields ${\mathbb{R}}^{N\times m}$, one output row per example, and the bias broadcasts down the $N$ rows. Reading a layer as a shape transformation is the fastest way to catch a bug before it runs:

Linear systems: $Ax=b$

A matrix times a vector is a system of linear equations. Solving $Ax=b$ asks which input $x$ the map $A$ sends to a target $b$. The answer is governed by three intertwined notions: span, linear independence, and rank.

For $Ax=b$ to have a solution for every $b\in {\mathbb{R}}^m$, the columns must span ${\mathbb{R}}^m$, requiring $n\ge m$ and $\mathrm{Rank}(A)=m$. For the solution to be unique, the columns must be linearly independent, requiring $n\le m$. Both at once force $A$ square ($m=n$) and full rank, the condition under which the inverse exists.

A matrix as a transformation of space

Geometrically, $A$ is a linear map: it sends the unit square (spanned by the basis vectors $e_1,e_2$) to the parallelogram spanned by $A$'s columns. Full rank means the parallelogram has nonzero area: the map does not collapse a dimension.

The signed area of that parallelogram is the determinant, met below. It records exactly how much $A$ scales volumes, and it is zero precisely when $A$ is singular.

Norms: measuring size

A norm assigns a nonnegative length to a vector, with $\Vert x\Vert =0$ only at $x=0$, scaling $\Vert \alpha x\Vert =\left|\alpha \right|\Vert x\Vert$, and the triangle inequality $\Vert x+y\Vert \le \Vert x\Vert +\Vert y\Vert$. The $L^p$ family covers most of what a network needs.

The squared $L^2$ norm $x^{\top }x$ is preferred in objectives because its gradient $2x$ is cheap and it is differentiable everywhere, unlike $\Vert x{\Vert }_2$, whose derivative misbehaves at the origin. The $L^1$ norm is what induces sparsity: because it grows at the same rate ($1$ per unit) no matter how close a coordinate is to zero, it pushes coordinates exactly to zero, which $L^2$ does not.

Unit balls

The set $\left\{x:\Vert x\Vert =1\right\}$, the surface of the unit ball, pictures each norm. The shapes explain the sparsity behavior: regularizing toward the origin pulls a solution to the boundary of one of these balls, and the $L^1$ diamond's corners sit on the axes, where coordinates vanish.

Special matrices

Three structured forms recur because each makes a hard operation cheap.

An orthogonal matrix is a rotation or reflection: it moves points without distorting them, since $\Vert Ax{\Vert }_2^2=x^{\top }{A}^{\top }Ax=x^{\top }x=\Vert x{\Vert }_2^2$. Diagonal matrices are the cheap case of every operation (multiply, invert, take powers), and that cheapness drives the decompositions below, which aim to turn a general matrix into a diagonal one in the right coordinates.

Eigendecomposition

An eigenvector of a square $A$ is a direction the map only stretches, never rotates; the stretch factor is its eigenvalue.

Stacking the eigenvectors as columns of $V$ and the eigenvalues on the diagonal of $\mathrm{diag}(\lambda )$, the relations $A{v}_i={\lambda }_i{v}_i$ assemble into the eigendecomposition:

This says $A$ does nothing more than: change to the eigenbasis ($V^{-1}$), scale each axis ($\mathrm{diag}(\lambda )$), change back ($V$). The picture is a space stretched along its eigenvector axes.

Eigenvalues come from the condition that $Av=\lambda v$ have a nonzero solution, i.e. that $(A-\lambda I)$ collapse a direction to zero, i.e. that $\det (A-\lambda I)=0$ — the characteristic polynomial. For $A=\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]$,

so ${\lambda }_1=3$, ${\lambda }_2=1$. Substituting ${\lambda }_1=3$ into $(A-3I)v=0$ gives $-v_1+v_2=0$, the direction $(1,1)$; ${\lambda }_2=1$ gives $(1,-1)$. Normalizing, the eigenvectors are $\frac{1}{\sqrt{2}}(1,1)$ and $\frac{1}{\sqrt{2}}(1,-1)$ — orthogonal, as the spectral theorem promises for a symmetric $A$. The map stretches the $(1,1)$ direction by $3$ and leaves $(1,-1)$ fixed. The trace $2+2=4$ equals ${\lambda }_1+{\lambda }_2$, and the determinant $4-1=3$ equals ${\lambda }_1{\lambda }_2$, the two scalar checks below.

Quadratic forms and positive-definiteness

The eigenvalues classify the quadratic form $f(x)=x^{\top }Ax$ (for symmetric $A$). In the eigenbasis, writing $y=Q^{\top }x$,

so the sign of the form is decided entirely by the signs of the eigenvalues.

Positive-definiteness is the multivariate second derivative is positive: a critical point where the Hessian is positive definite is a local minimum, where it is negative definite a maximum, and where it is indefinite a saddle point. This is the classification that drives the entire study of the optimization landscape.

Singular value decomposition

Eigendecomposition needs a square, diagonalizable matrix. The singular value decomposition has no such restriction; it factors any matrix.⁵

The columns of $U$ are left-singular vectors (eigenvectors of $A{A}^{\top }$), the columns of $V$ are right-singular vectors (eigenvectors of $A^{\top }A$), and the singular values are the square roots of the shared nonzero eigenvalues, ${\sigma }_i=\sqrt{{\lambda }_i(A^{\top }A)}$. Geometrically the factorization reads right-to-left as rotate, scale, rotate: $V^{\top }$ rotates the input, $D$ scales along axes, $U$ rotates the output. A unit circle becomes an ellipse whose axis lengths are the singular values.

In the caption, $s_1,s_2$ stand for the singular values ${\sigma }_1,{\sigma }_2$, the ellipse's semi-axis lengths. (Greek sigma is spelled out in the figure to keep the node text renderable.)

The Moore–Penrose pseudoinverse

When $A$ is not square, $Ax=b$ is over- or under-determined and has no exact inverse. The SVD supplies the next best thing: the pseudoinverse$A^{+}=V{D}^{+}{U}^{\top }$, where $D^{+}$ inverts each nonzero singular value ($1/{\sigma }_i$) and transposes. It returns the least-error, smallest-norm answer.

The least-squares case is the normal-equations solution $A^{+}=(A^{\top }A{)}^{-1}{A}^{\top }$ seen for linear regression, now derived through the SVD and valid even when $A^{\top }A$ is singular.

Trace and determinant

Two scalar summaries of a square matrix close the toolkit.

Worked example: principal components analysis

PCA is the canonical application of this machinery: a derivation in which the answer to a geometric optimization is a set of eigenvectors.⁶ The setup: $m$ points $x^{(1)},\dots ,x^{(m)}$ in ${\mathbb{R}}^n$ (assume centered, mean zero), to be compressed to $\ell <n$ dimensions with the least reconstruction error.

Encode each point by projecting onto an orthonormal set of directions stacked as columns of $D\in {\mathbb{R}}^{n\times \ell }$ ($D^{\top }D=I_{\ell }$), and decode by projecting back:

Choose $D$ to minimize the total squared reconstruction error, ${\sum }_i\Vert x^{(i)}-D{D}^{\top }{x}^{(i)}{\Vert }_2^2$.

The covariance $X^{\top }X$ is symmetric, so the spectral theorem guarantees a real, orthonormal eigenbasis: the principal axes. The first principal component is the direction of greatest variance; projecting onto it loses the least information.

The compression keeps the top $\ell$ axes and discards the rest. Because variance along the discarded axes is smallest, the data, projected onto an $\ell$-dimensional subspace, loses the least possible squared error, which is exactly the SVD's low-rank approximation guarantee restated for data.

Low rank as a design principle

Low rank has also become a design principle for large models, not only an analysis tool.

Low-rank adaptation. Fine-tuning a large model by updating every weight is expensive. LoRA (Hu et al., 2021) freezes a pretrained weight matrix $W$ and learns only a low-rank correction $W+BA$, where $B\in {\mathbb{R}}^{m\times r}$ and $A\in {\mathbb{R}}^{r\times n}$ with $r\ll \min (m,n)$. The update $BA$ has rank at most $r$, so it costs $r(m+n)$ parameters instead of $mn$ — for a $4096\times 4096$ layer and $r=8$, that is $65,536$ trainable numbers instead of $16.8$ million. The hypothesis, borne out empirically, is that the adaptation a task needs lies on a low-dimensional subspace even when the model does not: the SVD's rank structure applied to a weight update.

Randomized SVD. The exact SVD costs $O(mn\min (m,n))$, prohibitive for the huge matrices in modern pipelines. Randomized SVD (Halko, Martinsson, Tropp, 2011) projects the matrix onto a small random subspace, computes an exact decomposition of the tiny result, and recovers the top-$k$ singular vectors in $O(mnk)$ — the standard route to PCA on data too large for the full factorization.

Attention as a matrix product. Transformer attention is $\mathrm{softmax}(Q{K}^{\top }/\sqrt{d})V$: two matrix products with a row-wise softmax between them (Vaswani et al., 2017). Every idea in this lesson — shapes that must line up, the contraction that sums an inner dimension away, the $1/\sqrt{d}$ scaling that keeps the dot products in a stable range — appears in this one operation at the center of the field.⁷⁸

Takeaways

• Every quantity in a network is a tensor; layers are matrix products ($Wx+b$), gates and gradient masks are Hadamard products ($A\odot B$).
• $Ax=b$ is solvable iff $b$ lies in the column space; the inverse exists iff $A$ is square and full rank. Geometrically $A$ maps the unit square to a parallelogram whose signed area is $\det A$.
• Norms measure size; $L^1$ (diamond, sparsity), $L^2$ (circle, distance), $L^{\infty }$ (square, max), Frobenius (matrix). The unit-ball corners of $L^1$ explain why it zeroes coordinates.
• Eigendecomposition $A=V\mathrm{diag}(\lambda )V^{-1}$ diagonalizes a square map along its invariant axes; the spectral theorem makes this orthogonal for symmetric $A$, and the eigenvalue signs classify the quadratic form $x^{\top }Ax$ (positive-definite $\Rightarrow$ local minimum).
• SVD $A=UD{V}^{\top }$ exists for any matrix and reads as rotate–scale–rotate; it yields the pseudoinverse $A^{+}$, the least-squares / minimum-norm solver for non-invertible systems.
• PCA is the worked example: minimizing reconstruction error is maximizing projected variance, whose optimum is the top eigenvectors of $X^{\top }X$, the same statement as the SVD's low-rank approximation.