Chapter 5: Linear Algebra

Independence formalizes the idea of no wasted vectors. The formula says the only way to combine $\mathbf{v}_1, \ldots, \mathbf{v}_k$ with coefficients $c_i$ into the zero vector is to set every $c_i$ to zero. If some nontrivial combination gives zero, then any vector with a nonzero coefficient can be solved for in terms of the others, so it was redundant. In practice you stack the vectors as the columns of a matrix $A$ and check whether $\text{rank}(A) = k$; for a square matrix, $\det(A) \neq 0$ is equivalent. In $\mathbb{R}^n$ at most $n$ vectors can be independent. Floating point blurs the test, because nearly dependent columns make the rank depend on the tolerance.

Independence prevents redundancy and spanning guarantees reach, and a basis $\mathcal{B}$ has both. In the formula, $|\mathcal{B}|$ is the number of vectors in any basis of $V$. The theorem behind it, the exchange lemma, says this count never depends on which basis you pick, so the dimension is well defined. The standard basis of $\mathbb{R}^n$ uses $\mathbf{e}_i$, which has a 1 in slot $i$ and zeros elsewhere. For a subspace defined by equations, count the free variables: each one contributes one basis vector. A common mistake is to call any spanning set a basis. If it has more vectors than the dimension, it is dependent.

A basis turns abstract vectors into lists of numbers. Given a basis $\mathcal{B}$ with vectors $\mathbf{b}_i$, the coefficients $c_i$ in the formula are the coordinates of $\mathbf{v}$, and independence makes them unique. To compute them, place the basis vectors in the columns of a matrix $B$ and solve $B\mathbf{c} = \mathbf{v}$. For an orthonormal basis this collapses to dot products, $c_i = \mathbf{b}_i \cdot \mathbf{v}$. Changing the basis changes the numbers, not the vector: PCA rewrites each data point in the eigenvector basis of the covariance matrix. The usual pitfall is mixing coordinates taken from two different bases, which silently compares incompatible numbers.

The column space answers the question of what a matrix can output. Writing $A\mathbf{x}$ as $x_1\mathbf{a}_1 + \cdots + x_n\mathbf{a}_n$, where $\mathbf{a}_j$ is column $j$, shows that every output is a combination of the columns, so the column space is their span. Consequently $A\mathbf{x} = \mathbf{b}$ is solvable if and only if $\mathbf{b}$ lies in that span. Its dimension is the rank, the number of pivots found by row reduction, and the pivot columns of the original $A$ form a basis. A 3×2 matrix can reach at most a plane in $\mathbb{R}^3$, so most right-hand sides have no exact solution, which is the situation least squares will handle.

The null space (or kernel) of $A$ is the set of all vectors that $A$ maps to zero. If the null space contains only the zero vector, then $A$ is injective (one-to-one), and solutions to $A\mathbf{x} = \mathbf{b}$ are unique when they exist. A larger null space means more 'degrees of freedom' in the solution. In neural networks, the null space of a layer's weight matrix represents input directions that the layer ignores — understanding this helps diagnose vanishing gradients and dead neurons.

Every column of $A$ ends up doing one of two jobs. Either it adds a new independent direction to the output, raising the rank, or it is a combination of earlier columns, which creates a free variable and one more dimension of null space. The theorem makes this bookkeeping exact: $\text{rank}(A) + \text{nullity}(A) = n$, where $n$ counts columns, not rows. The practical reading is simple. If the rank equals $n$, the null space contains only the zero vector and solutions are unique; any shortfall in rank shows up as a family of solutions with that many dimensions. A frequent slip is subtracting the rank from the number of rows instead.

Orthogonality is the cleanest possible relationship between two directions. The formula says $\mathbf{u} \perp \mathbf{v}$ exactly when the dot product, the sum of $u_i v_i$, is zero; in two or three dimensions that means a right angle, because $\mathbf{u} \cdot \mathbf{v} = \|\mathbf{u}\|\|\mathbf{v}\|\cos\theta$. For orthogonal vectors Pythagoras holds, so squared lengths simply add. A set of mutually orthogonal nonzero vectors is automatically independent: dot any combination equal to zero with $\mathbf{u}_j$, and every term except $c_j\|\mathbf{u}_j\|^2$ vanishes, forcing $c_j = 0$. The zero vector is orthogonal to everything, so it must be excluded before calling such a set a basis.

To find the closest point to $\mathbf{v}$ on the line through $\mathbf{u}$, look for a multiple $t\mathbf{u}$ whose error $\mathbf{v} - t\mathbf{u}$ is orthogonal to $\mathbf{u}$. Setting $(\mathbf{v} - t\mathbf{u}) \cdot \mathbf{u} = 0$ and solving for $t$ gives exactly the formula. The numerator measures how much $\mathbf{v}$ points along $\mathbf{u}$, and the denominator removes the effect of the length of $\mathbf{u}$, so rescaling $\mathbf{u}$ leaves the projection unchanged. The result splits $\mathbf{v}$ into a part along $\mathbf{u}$ plus a perpendicular residual. Forgetting the denominator is a common bug; it is only safe when $\mathbf{u}$ has unit length.

For a subspace $W$ spanned by the columns of $A$, the projection has the form $A\hat{\mathbf{x}}$ for some coefficient vector $\hat{\mathbf{x}}$, and the residual $\mathbf{b} - A\hat{\mathbf{x}}$ must be orthogonal to every column, which means $A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}$. Solving gives the coefficients in the formula, and the projection matrix is $P = A(A^TA)^{-1}A^T$. It satisfies $P^2 = P$, since projecting twice changes nothing, and $P^T = P$. The formula requires independent columns; otherwise $A^TA$ is singular. When the columns are orthonormal, $A^TA = I$ and $P$ simplifies to $AA^T$, which is why orthonormal bases are prized.

Given a subspace $W$ of $V$, its orthogonal complement collects every vector whose dot product with all of $W$ is zero. It is itself a subspace, and the dimensions add up, as $\dim(W) + \dim(W^\perp) = \dim(V)$. Every vector splits uniquely as $\mathbf{w} + \mathbf{w}^\perp$, where $\mathbf{w}$ is the projection onto $W$ and $\mathbf{w}^\perp$ is the residual. The most useful instance comes from $A\mathbf{x} = \mathbf{0}$, which says $\mathbf{x}$ has zero dot product with each row of $A$, so the null space is the orthogonal complement of the row space. To test membership it is enough to check a spanning set of $W$, not every vector in it.

Gram-Schmidt builds the orthonormal basis one vector at a time. For the $k$-th input $\mathbf{v}_k$, the formula subtracts its projection onto each earlier direction $\mathbf{u}_j$, leaving $\mathbf{u}_k$, the part of $\mathbf{v}_k$ that is genuinely new; dividing by its length gives the unit vector $\mathbf{e}_k$. Because only earlier directions are removed, the first $k$ outputs span the same space as the first $k$ inputs. If some $\mathbf{u}_k$ comes out zero, then $\mathbf{v}_k$ depended on its predecessors. The classical version projects the original $\mathbf{v}_k$ and loses orthogonality in floating point when inputs are nearly parallel. The modified version subtracts from the running residual instead.

QR records the bookkeeping of Gram-Schmidt. Each column $\mathbf{a}_k$ of $A$ is a combination of $\mathbf{e}_1$ through $\mathbf{e}_k$ only, with coefficients $r_{jk} = \mathbf{e}_j \cdot \mathbf{a}_k$ and diagonal entry $r_{kk} = \|\mathbf{u}_k\|$. Stacking the $\mathbf{e}_j$ as the columns of $Q$ and the coefficients into $R$ gives $A = QR$, where $Q^TQ = I$ and $R$ is upper triangular because later directions never appear in earlier columns. To solve least squares, multiply by $Q^T$ and back-substitute in $R\mathbf{x} = Q^T\mathbf{b}$. Production libraries compute QR with Householder reflections rather than Gram-Schmidt, because reflections keep $Q$ orthogonal to machine precision.

Least squares minimizes $\|A\mathbf{x} - \mathbf{b}\|^2$, the squared length of the residual. Two derivations lead to the same answer. With calculus, the gradient $2A^T(A\mathbf{x} - \mathbf{b})$ must vanish. With geometry, the best $A\hat{\mathbf{x}}$ is the projection of $\mathbf{b}$ onto the column space, so the residual is orthogonal to every column. Either way we get the normal equations in the formula, an $n \times n$ system however many rows $A$ has. If the columns of $A$ are independent, $A^TA$ is invertible and $\hat{\mathbf{x}}$ is unique. There are two failure modes: dependent columns make $A^TA$ singular, and forming $A^TA$ squares the condition number, which QR avoids.

The pseudoinverse extends inversion to every matrix. Write the SVD $A = U\Sigma V^T$, with orthonormal $U$ and $V$ and the singular values $\sigma_i$ on the diagonal of $\Sigma$. Then $\Sigma^+$ replaces each nonzero $\sigma_i$ by $1/\sigma_i$, leaves the zeros alone, and transposes the shape, and the formula assembles $A^+$. The vector $A^+\mathbf{b}$ minimizes the residual and, among all minimizers, has the smallest norm. For full column rank it equals $(A^TA)^{-1}A^T$. In floating point, singular values near zero must be cut off by a tolerance, NumPy's rcond, because inverting a tiny $\sigma_i$ turns noise into enormous coefficients.

Linear regression is least squares under statistical names. Each row of the design matrix $X$ is one example and each column one feature, plus a column of ones for the intercept; $\mathbf{y}$ holds the targets and $\hat{\boldsymbol{\beta}}$ the fitted weights, so the formula is the normal equations in new letters. Under Gaussian noise this estimate is also the maximum-likelihood one, which ties it back to the probability chapter. Ridge regression adds the penalty $\lambda\|\boldsymbol{\beta}\|^2$, giving $(X^TX + \lambda I)^{-1}X^T\mathbf{y}$. Adding $\lambda I$ lifts every eigenvalue of $X^TX$ by $\lambda$, so the system stays invertible for any positive $\lambda$, even with collinear features, at the price of shrinking the weights.

General square matrices can have complex eigenvalues or too few eigenvectors, but symmetric ones are always well behaved. The theorem says a real symmetric $A$ has real eigenvalues $\lambda_i$ and orthonormal eigenvectors, collected as the columns of $Q$, so $Q^TQ = I$ and $\Lambda$ holds the eigenvalues on its diagonal. Read $Q\Lambda Q^T\mathbf{x}$ from right to left: rotate into eigenvector coordinates, scale each coordinate by its eigenvalue, and rotate back. Eigenvectors for distinct eigenvalues are automatically orthogonal. Covariance matrices are symmetric, so PCA is exactly this factorization. For non-symmetric matrices the guarantee fails: the rotation [[0, -1], [1, 0]] has eigenvalues $\pm i$.

A quadratic form assigns each vector the number $\mathbf{x}^TA\mathbf{x}$, a polynomial whose terms all have degree two. Only the symmetric part of $A$ matters, so we take $A$ symmetric and split each cross term evenly; a term $4xy$ puts 2 in both off-diagonal slots. Substituting $\mathbf{x} = Q\mathbf{z}$ from the spectral theorem gives $f = \sum_i \lambda_i z_i^2$, a pure sum of squares. The eigenvalue signs now decide the shape: all positive gives a bowl, all negative a dome, mixed signs a saddle, and a zero eigenvalue a flat trough. Judging from the diagonal alone fails, since [[1, 3], [3, 1]] has a positive diagonal yet an eigenvalue of -2.

Positive definiteness is the matrix version of a positive number. The formula lists two equivalent tests: $\mathbf{x}^TA\mathbf{x} > 0$ for every nonzero $\mathbf{x}$, or all eigenvalues positive. Further equivalents are Sylvester's criterion, which requires all leading principal minors to be positive, and the existence of a Cholesky factor $A = LL^T$, which is also the cheapest numerical test: attempt it and see whether it fails. Semi-definite matrices allow zero eigenvalues. Covariance and kernel matrices are always semi-definite, since $\mathbf{x}^T\Sigma\mathbf{x}$ is a variance. At a critical point of a loss, a positive definite Hessian certifies a strict local minimum, while a single negative eigenvalue marks a saddle.

The Rayleigh quotient measures how much a matrix 'stretches' in a given direction. Its maximum is the largest eigenvalue (achieved at the corresponding eigenvector), and its minimum is the smallest. This is exactly how PCA works: the first principal component maximizes the Rayleigh quotient of the covariance matrix. It also appears in spectral clustering, Fisher's linear discriminant, and the power iteration method used to find dominant singular vectors in large-scale SVD (as used by recommendation engines).

Before trusting a computed answer, ask how sensitive the problem itself is. The condition number $\kappa(A) = \|A\| \cdot \|A^{-1}\|$, which in the 2-norm equals the ratio of the largest to the smallest singular value, bounds the amplification: a relative change $\epsilon$ in $\mathbf{b}$ can move $\mathbf{x}$ by up to $\kappa(A)\epsilon$ in relative terms. With unit roundoff near $10^{-16}$, a condition number of $10^k$ leaves about $16 - k$ trustworthy digits. Orthogonal matrices have $\kappa = 1$, the ideal; the 8×8 Hilbert matrix has $\kappa \approx 1.5 \times 10^{10}$. Conditioning belongs to the problem, not the algorithm, so no solver can recover digits an ill-conditioned matrix has already destroyed.

IEEE 754 float64 stores 53 significant bits, so the gap between 1 and the next representable number, machine epsilon, is $2^{-52} \approx 2.2 \times 10^{-16}$, while float32 has $\varepsilon \approx 1.2 \times 10^{-7}$. Each operation is exact up to a relative error of about $\varepsilon$, but errors compound. Subtracting nearly equal numbers cancels leading digits, and adding a tiny number to a huge one loses it entirely: in float64, $10^{16} + 1$ evaluates to $10^{16}$. Three rules follow. Compare floats with a tolerance, never with ==. Solve with LU, QR, or Cholesky rather than forming $A^{-1}$. In FP16 training, scale the loss so small gradients do not underflow to zero.

Direct factorizations cost $O(n^3)$ time and fill in sparse matrices, which is hopeless for a million unknowns. Iterative methods instead refine a guess using only products $A\mathbf{v}$. The formula shows the template: compute the residual $\mathbf{b} - A\mathbf{x}_k$, apply an easy approximate inverse $M^{-1}$ called the preconditioner, and correct $\mathbf{x}_k$. Jacobi takes $M$ to be the diagonal of $A$. Conjugate gradient, for symmetric positive definite $A$, chooses optimal search directions and needs roughly $\sqrt{\kappa}$ iterations per digit gained, and GMRES handles general matrices. Without preconditioning, a large $\kappa$ makes them stall for thousands of steps, so a good $M$ often matters more than the method.

Large matrices from graphs, meshes, and text are mostly zeros. Sparse formats keep only the nnz nonzero entries, so storage and a matrix-vector product cost $O(\text{nnz})$ instead of $O(n^2)$. COO stores (row, column, value) triples and is easy to build. CSR packs each row's column indices and values contiguously plus a row-pointer array, which makes row slicing and $A\mathbf{v}$ fast; CSC does the same by columns. A TF-IDF matrix with a 50,000-word vocabulary and about 200 distinct words per document is over 99.5 percent zeros. The trap is fill-in: LU can turn a sparse matrix dense, so iterative solvers suit sparse storage.