24.5 Further Reading
The sections of this chapter teach the geometry of one inner product, maps read through their columns, the eigenvectors of one square matrix, and the singular value decomposition with its best low-rank approximation. This section collects, in the order of the chapter, topics they had no room for. Each topic states the problem it addresses and what it changes relative to the chapter. Books and courses for the material that the chapter does teach are listed in its introduction (Chapter 24). Nothing else in the book depends on this section.
24.5.1 Geometry and Linear Maps
24.5.1.1 Other Inner Products: Kernel Methods
The nearest-mean rule of Section 24.1.4.2 compares an image with two class means through one dot product, so its decision boundary is a hyperplane. It fails when a class is spread out or curved, as on t-shirts against shirts. Kernel methods keep the rule and change the inner product: \(\mathbf{x}\cdot\mathbf{x}_i\) is replaced by a function \(k(\mathbf{x},\mathbf{x}_i)\) that is an inner product in some other, possibly infinite-dimensional, space. The class mean becomes a mean in that space, and the rule becomes a Parzen-window classifier whose boundary can bend. Choosing which training points to keep, and with what weights, leads to the support vector machine. Schölkopf and Smola (2002) start their book from exactly this classifier and develop the theory of kernels from it.
24.5.1.2 Random Projections
Section 24.1 showed that random directions in \(\mathbb{R}^d\) are nearly orthogonal. Can this shrink the dimension of a dataset without changing its geometry? The Johnson–Lindenstrauss lemma (Johnson and Lindenstrauss 1984) says yes: for \(n\) points, a scaled random Gaussian map to \(k\) of order \(\varepsilon^{-2}\log n\) dimensions preserves all \(\binom{n}{2}\) pairwise distances within a factor \(1\pm\varepsilon\), with high probability. Unlike principal component analysis (Section 24.4.2.2), the map does not look at the data, and the target dimension does not depend on \(d\), so the lemma pays off when \(d\) is far larger than \(\varepsilon^{-2}\log n\). Dasgupta and Gupta (2003) give a short proof with explicit constants, and Vershynin (2018) places the lemma among the concentration results for random vectors.
24.5.1.3 Elimination, Determinants, and Rank in Practice
Section 24.2 defines the rank by independent columns and the determinant by a sum over \(n!\) permutations, and neither definition is ever run. Libraries use Gaussian elimination with partial pivoting, the LU factorization behind solve. The determinant is the product of the pivots, times the sign of the row permutation, at the cost of one factorization. The rank is read from the singular values with a tolerance, as np.linalg.matrix_rank does, or more cheaply from a QR factorization with column pivoting. Trefethen and Bau (1997) explain why partial pivoting is stable in practice although it is not in the worst case, and Higham (2002) gives the full rounding-error analysis of these algorithms.
24.5.1.4 Sparse and Iterative Solvers
A \(10^6\times10^6\) matrix with five nonzeros per row is cheap to store, yet its inverse is typically dense, with \(10^{12}\) entries. Solvers for such matrices never form the inverse, and often not the LU factors either. Sparse direct methods reorder the unknowns so that elimination creates few new nonzeros. Iterative methods touch the matrix only through products \(\mathbf{A}\mathbf{x}\), as power iteration does (Section 24.3.3.2). For a symmetric positive definite matrix, the conjugate gradient method of Hestenes and Stiefel (1952) shrinks the error after \(k\) steps by at least the factor \(2\bigl((\sqrt\kappa-1)/(\sqrt\kappa+1)\bigr)^k\) in the norm defined by the matrix, where \(\kappa\) is the condition number of Section 24.4.3.2. Golub and Van Loan (1996) cover both families.
24.5.1.5 Contraction Order and Tensor Decompositions
The einsum of Section 24.2.4 fixes the result of a product of several factors but not the order of the pairwise contractions. For \(n\times n\) matrices, \((\mathbf{A}\mathbf{B})\mathbf{x}\) costs about \(n^3\) multiplications and \(\mathbf{A}(\mathbf{B}\mathbf{x})\) costs \(2n^2\). Finding the best order for many factors is a search over parenthesizations, which np.einsum(..., optimize=True) and the opt_einsum package (Smith and Gray 2018) carry out. Low rank also extends beyond matrices. A three-way tensor can be written as a sum of outer products \(\mathbf{a}_i\otimes\mathbf{b}_i\otimes\mathbf{c}_i\) (the CP decomposition), as a small core tensor multiplied by a matrix along each mode (Tucker), or as a chain of three-way cores (the tensor train of Oseledets (2011)). The one difference from matrices is that Eckart–Young has no analogue: computing the CP rank is NP-hard, and a best approximation of a given CP rank need not exist. The survey of Kolda and Bader (2009) covers these decompositions and their algorithms.
24.5.2 Eigenvalues and Spectra
24.5.2.1 How Libraries Compute Eigenvalues
No formula in radicals gives the roots of a general polynomial of degree five or more, so every eigensolver iterates, and none forms the characteristic polynomial of Section 24.3.1.2, whose roots are highly sensitive to its coefficients. eig first reduces the matrix by orthogonal similarity transformations to Hessenberg form, zero below the first subdiagonal, at a cost of \(O(n^3)\). The shifted QR algorithm then factors \(\mathbf{H}-\mu\mathbf{I}=\mathbf{Q}\mathbf{R}\), forms \(\mathbf{R}\mathbf{Q}+\mu\mathbf{I}\), and repeats. Each step is a similarity transformation, so the eigenvalues stay the same, while the subdiagonal entries go to zero and the diagonal converges to the spectrum. For a symmetric matrix the Hessenberg form is tridiagonal, which eigh exploits. Trefethen and Bau (1997) derive the algorithm from power iteration applied to whole subspaces, and Golub and Van Loan (1996) give the production details.
24.5.2.2 Pseudospectra and the Kreiss Matrix Theorem
Section 24.3.4.4 showed that the spectrum of a non-normal matrix does not predict its transient growth. The Kreiss matrix theorem makes the pseudospectra quantitative. With the Kreiss constant \(K(\mathbf{A})=\sup_{|z|>1}(|z|-1)\|(z\mathbf{I}-\mathbf{A})^{-1}\|_2\), it brackets the largest power, \(K(\mathbf{A})\le\sup_k\|\mathbf{A}^k\|_2\le e\,n\,K(\mathbf{A})\) for an \(n\times n\) matrix. The lower bound is an exercise of Section 24.3; the upper bound says that a matrix whose pseudospectra stay close to the unit disc cannot grow much. Trefethen and Embree (2005) develop the subject and its applications, among them the transition to turbulence in fluid flows whose linearization has only stable eigenvalues.
24.5.2.3 Mixing Times and the Spectral Gap
The two-page chain of Section 24.3 approaches its stationary distribution at rate \(1/4\), its second eigenvalue. How long does a general chain take? For a reversible chain, let \(\lambda_\star\) be the largest modulus among the eigenvalues other than 1. The number of steps needed to come within total variation distance \(\varepsilon\) of the stationary distribution is bounded above and below by multiples of \(1/(1-\lambda_\star)\), the multiples depending on \(\varepsilon\) and on the smallest stationary probability. Cheeger’s inequality bounds the gap \(1-\lambda_2\) from both sides by the chain’s conductance, the smallest probability flow out of a set relative to the set’s mass, so the gap can be read from the shape of the graph. Levin and Peres (2017) develop this theory and prove the Perron–Frobenius facts that Section 24.3.3.3 cites from Horn and Johnson (2012) .
24.5.2.4 Random Matrices and Trainable Depth
Section 24.3.4.5 showed that the spectral radii of the factors do not control a product of Jacobians. At initialization, what does? Saxe et al. (2014) solve the learning dynamics of deep linear networks exactly and find that orthogonal initial weights make the learning time independent of depth, where scaled Gaussian weights do not. Pennington et al. (2017) use free probability to compute the distribution of all singular values of the input–output Jacobian of a deep nonlinear network. Dynamical isometry, all these singular values near 1, is reachable with orthogonal weights and tanh units but not with ReLU, and networks that have it train faster at large depth.
24.5.2.5 Spectra of Trained Weights
The Marchenko–Pastur law of Section 24.3.2.4 describes the spectrum of a matrix of independent noise, and the weights of an untrained network follow it. Martin and Mahoney (2021) compare the spectra of trained weight matrices with this null model and find that training produces outliers and a heavy, power-law tail. They read the exponent of that tail as a measure of implicit self-regularization. The departure from Marchenko–Pastur is observed across many models; using the exponent to rank models by quality without test data rests on an empirical correlation, not a theorem.
24.5.3 Low Rank Beyond the Truncated SVD
24.5.3.1 Matrix Completion and the Nuclear Norm
A ratings matrix is observed in a few percent of its entries. Eckart–Young (Section 24.4.2.1) needs the whole matrix, so it does not apply. The one change is to minimize the nuclear norm \(\|\mathbf{Z}\|_*=\sum_i\sigma_i\), the convex stand-in for rank in the way that \(\ell_1\) stands in for sparsity, subject to agreeing with the observed entries. Candès and Recht (2009) prove that this convex program recovers a low-rank matrix exactly, with high probability, from far fewer than all \(n^2\) entries chosen at random, provided its singular vectors are spread out over the coordinates. Soft-Impute (Mazumder et al. 2010), the algorithm of the matrix-completion exercise of Section 24.4, solves the penalized form by alternating a fill-in step with soft-thresholding of the singular values.
24.5.3.2 CUR Decompositions and Sketching
The factors of the truncated SVD mix all rows and columns, so they are hard to interpret, and a full SVD of an \(m\times n\) matrix costs \(O(mn\min(m,n))\). A CUR decomposition keeps actual columns \(\mathbf{C}\) and rows \(\mathbf{R}\) of the data and fits a small core \(\mathbf{U}\) between them. Section 24.4.2.1 showed the naive choice, the columns of largest norm, losing to the truncation by a wide margin. Mahoney and Drineas (2009) sample columns and rows by their leverage scores, the squared lengths of their components along the top singular vectors, and come within a small constant factor of the Eckart–Young error with high probability. The range finder of Section 24.4.4.3 is one instance of sketching: multiply by a random matrix and work with the small result. Woodruff (2014) applies the same idea to least squares and regression, which are solved to within a factor \(1+\varepsilon\) on a sketch with far fewer rows than the data.
24.5.3.3 Adapter Variants
LoRA (Section 24.4.4.1) starts its update \(\mathbf{B}\mathbf{A}\) at zero, so early fine-tuning moves along random directions. PiSSA (Meng et al. 2024) changes the initialization: \(\mathbf{B}\mathbf{A}\) starts as the rank-\(k\) truncated SVD of the pretrained weight, and the residual is frozen, so training adjusts the principal directions first. DoRA (Liu et al. 2024) changes the parameterization instead. It writes each column of the weight as a length times a unit direction, trains the lengths directly, and applies the low-rank update to the directions only.
24.5.3.4 Low Rank in Attention
Attention over \(n\) tokens forms an \(n\times n\) matrix, and several methods taught in the attention chapters replace it, or the cache behind it, by a low-rank factorization. Linear attention (Section 10.5) writes the scores as \(\phi(\mathbf{q})\cdot\phi(\mathbf{k})\), so the matrix has rank at most the feature dimension, and multi-head latent attention (Section 11.3) (DeepSeek-AI 2024) caches one low-dimensional latent per token. Linformer (Wang et al. 2020) starts from the observation that trained attention matrices have rapidly decaying singular values and projects the keys and values along the sequence axis to a fixed length.
24.5.3.5 Steepest Descent under Other Norms
Section 24.4.4.2 identified the polar factor \(\mathbf{U}\mathbf{V}^\top\) as the steepest-descent direction when steps are measured in the spectral norm. Bernstein and Newhouse (2024) turn this into a way of classifying optimizers. With their moving averages switched off, Adam becomes sign descent, which is steepest descent under the \(\ell_\infty\) norm, and Shampoo becomes steepest descent under the spectral norm. The norm, chosen per layer, is then the design decision, and the update rule follows from it. The optimizer built on the spectral norm is taught in Section 9.7.