math.NAAug 21, 2026
SaveAdvanced Linear Algebra with Applications - Part I (Numerical linear algebra for PDEs, machine learning, and data assimilation)
Abstract
These lecture notes form the first part of a master's-level course on advanced numerical linear algebra. Their aim is not only to present the classical algorithms, but to show why the subject has become considerably more central than it was a generation ago. Numerical linear algebra grew up alongside the numerical solution of partial differential equations, and for a long time that is where its large sparse systems came from. Ranking the nodes of a network, assimilating observations into a weather forecast, and fitting a model to a large noisy data set now lead to problems of the same kind: too large to factorise, structured, and accessible only through matrix-vector products. Strikingly few ideas are needed for all of them. Each chapter therefore develops a standard topic and then puts it to work outside its original setting. We treat norms, factorisations, conditioning and floating-point arithmetic; sparse matrices arising from finite differences, from graphs and from machine learning; stationary iterations and the smoothing property; the conjugate gradient and Lanczos methods, with spectral clustering and regularisation by early stopping; Arnoldi and GMRES, with PageRank and large least squares; and finally preconditioning, Schwarz domain decomposition and multigrid. We assume a first course in linear algebra. Every section closes with a summary of what should be retained and every chapter with exercises, several drawn from past examinations. Accompanying Python code reproduces the numerical illustrations.
Figures & tables
Figure 1.1 : Unit balls for the , , and norms in : the same vector has a different "size" depending on which norm is used.
Figure 2.1 : Zoom-in of around a point with neighbors
Figure 2.2 : Left: Function values on discretised mesh. Right: Stencil for central differences in 1D
Figure 2.3 : Left: Indexing of grid. Right: stencil to update
Figure 2.4 : Visualisation of fill-in with LU vs inverse. Top left: sparsity pattern original matrix, . Top right: sparsity pattern of . Lower left: sparsity pattern of . Lower right: sparsity pattern of .
Figure 4.1 : Steepest descent versus Conjugate Gradient on the level sets of with and the same starting point (iterates computed exactly). Because consecutive steepest-descent residuals are orthogonal, the iterates zig-zag across the valley and creep towards ; the number of steps grows like , so the picture gets dramatically worse as the contours become more elongated. CG replaces orthogonality of the residuals by -conjugacy of the search directions : no direction is ever revisited, and in exact arithmetic steps suffice.
Figure 4.2 : The optimality property of CG. The ellipses are level sets of the energy norm . At step the iterate is the point of the affine space (dashed line) closest to in the -norm, i.e. the -orthogonal projection of onto it.
Figure 4.3 : Error decay of steepest descent and CG for the same SPD matrix. Both are linear (geometric) rates, but CG’s contraction factor is governed by rather than . For the 1D Poisson matrix, where , this is the difference between and iterations. The bound is only an upper estimate: when the eigenvalues of are clustered, CG typically converges considerably faster than the curve shown.
| Method | Convergence | Iterations (PDE) | Remarks |
|---|---|---|---|
| Jacobi | Very slow; mainly used as smoother | ||
| Gauss–Seidel | Roughly twice as fast as Jacobi | ||
| SOR (optimal ) | Improved, but still asymptotically slow | ||
| Steepest descent | Gradient method; zig-zag convergence | ||
| Conjugate Gradient | Optimal Krylov solver |
Table 4.1 : Comparison of classical stationary methods, steepest descent, and Conjugate Gradient for SPD PDE systems.
Figure 4.4 : The Lanczos idea in one line. The large sparse symmetric matrix is never factorised; it is only applied to vectors. After steps the problem has been compressed into an symmetric tridiagonal matrix whose eigenvalues (the Ritz values ) approximate the extreme eigenvalues of , and whose eigenvectors lift back to approximate eigenvectors of through .
Figure 4.5 : Left: the truncated SVD cuts small- directions off abruptly, ridge damps them smoothly, and the CG filter lies between the two. Right: semi-convergence — the unresolved part (dotted) falls while the admitted noise (dashed) grows, so the total (solid) is smallest at .
Figure 5.1 : The Arnoldi decomposition , drawn for . Applying to the orthonormal basis of never leaves , so the result factors through exactly one extra basis vector. Column of is nonzero only in its first entries, which is what upper Hessenberg means; in particular the first column has just two nonzeros. The matrix is , the last row containing the single entry that carries . When is symmetric, is additionally symmetric and therefore tridiagonal, and Arnoldi reduces to the Lanczos algorithm of Chapter 4.
Figure 5.2 : Why eigenvalues are not enough. For a normal matrix the eigenvalues alone control , so a tight cluster away from the origin guarantees rapid GMRES convergence. For a non-normal matrix the eigenvector basis is ill-conditioned and the relevant object is the field of values , a convex set containing . If stays in a disk that excludes the origin, the Crouzeix–Palencia inequality [ 8 ] gives a genuine convergence rate; if the origin lies inside , no such bound is available, and GMRES can stagnate for steps even with all eigenvalues equal to (see the example at the end of this section).
| power iterations to | ||
|---|---|---|
Table 5.1 : The damping factor on the three-node example of this section, where exactly. The last column is the number of iterations observed to gain eight digits, in agreement with the asymptotic estimate . The value used by Google is a compromise between fidelity to the link structure and cost, not a mathematical constant.
Figure 5.3 : Least squares as orthogonal projection. The minimiser of is the point of the column space closest to , so the residual is orthogonal to , which leads exactly to the normal equations . LSQR realises the same projection without forming .
Figure 5.4 : Golub–Kahan bidiagonalisation. Two orthonormal sequences are grown in alternation, one in the row space and one in the column space, each step costing one product with and one with . The coefficients assemble into a lower bidiagonal with . This is Lanczos applied implicitly to but the squaring never happens numerically, which is precisely why LSQR is stable where CG on the normal equations is not.
Figure 5.5 : LSQR as an iterative truncated SVD. Because the Krylov space is built from and grows by powers of , components associated with large singular values are resolved first. Stopping early therefore acts as a spectral filter, suppressing the small- directions that magnify noise in an ill-posed problem. The iteration count plays the role of a regularisation parameter, the same mechanism as truncating the SVD, but obtained without ever computing one.
Figure 5.6 : After truncated SVD, documents cluster in a low-dimensional semantic space.
Figure 6.1 : Effect of preconditioning on the spectrum: eigenvalues of (left) are spread, whereas those of (right) are tightly clustered around 1.
Figure 6.2 : 1D domain split into two overlapping subdomains.
Figure 6.3 : Alternating Schwarz converges geometrically in 1D. The contraction factor improves when the overlap increases.
Figure 6.4 : Partition of unity in the overlap: (top) and (bottom). On the overlap, and outside the overlap they sum to one: .
Figure 6.5 : Low-frequency and high-frequency error modes on a 1D grid.
Figure 6.6 : Why a coarse space can carry the correction. The thin grey curve is a random initial error on a grid with ; the thick blue curve is the same error after five weighted-Jacobi sweeps ( ); the dashed orange curve is the piecewise-linear interpolant through its values at the coarse nodes only.
Figure 6.7 : The two-grid method (left) and the V-cycle obtained by applying it recursively (right). Each descending edge is a restriction of the residual, each ascending edge a prolongation of the correction, and relaxation is applied before and after every transfer. Recursion is what makes the method optimal: the coarse solve is itself replaced by a two-grid step until the problem is small enough to solve directly, so one cycle costs work and reduces the error by a factor independent of . Used as a preconditioner for CG, this turns the iteration count of Chapter 4 into a constant.