math.NAAug 21, 2026

Advanced Linear Algebra with Applications - Part I (Numerical linear algebra for PDEs, machine learning, and data assimilation)

Authors: Victorita Dolean, Jemima Tabeart

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

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 {x:∥x∥=1}\{x:\|x\|=1\} for the ℓ1\ell^{1} , ℓ2\ell^{2} , and ℓ∞\ell^{\infty} norms in R2\mathbb{R}^{2} : the same vector has a different "size" depending on which norm is used.
Figure 2.1 : Zoom-in of u(x)u(x) around a point with neighbors x±hx\pm h
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 u(i,j)u(i,j)
Figure 2.4 : Visualisation of fill-in with LU vs inverse. Top left: sparsity pattern original matrix, AA . Top right: sparsity pattern of A−1A^{-1} . Lower left: sparsity pattern of LL . Lower right: sparsity pattern of UU .
Figure 4.1 : Steepest descent versus Conjugate Gradient on the level sets of J(x)=12x⊤Ax−b⊤xJ(x)=\tfrac{1}{2}x^{\top}Ax-b^{\top}x with A=diag(1,9)A=\mathrm{diag}(1,9) 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 x∗x^{\ast} ; the number of steps grows like κ(A)\kappa(A) , so the picture gets dramatically worse as the contours become more elongated. CG replaces orthogonality of the residuals by AA -conjugacy of the search directions : no direction is ever revisited, and in exact arithmetic nn steps suffice.
Figure 4.2 : The optimality property of CG. The ellipses are level sets of the energy norm ∥x−x∗∥A\|x-x^{\ast}\|_{A} . At step kk the iterate is the point of the affine space x(0)+Kk(A,r(0))x^{(0)}+\mathcal{K}_{k}(A,r^{(0)}) (dashed line) closest to x∗x^{\ast} in the AA -norm, i.e. the AA -orthogonal projection of x∗x^{\ast} 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 κ(A)\sqrt{\kappa(A)} rather than κ(A)\kappa(A) . For the 1D Poisson matrix, where κ(A)=O(n2)\kappa(A)=\mathcal{O}(n^{2}) , this is the difference between O(n2)\mathcal{O}(n^{2}) and O(n)\mathcal{O}(n) iterations. The bound is only an upper estimate: when the eigenvalues of AA are clustered, CG typically converges considerably faster than the curve shown.
MethodConvergenceIterations (PDE)Remarks
Jacobi(1−cn−2)k\big(1-cn^{-2}\big)^{k}O(n2)\mathcal{O}(n^{2})Very slow; mainly used as smoother
Gauss–Seidel(1−cn−2)k\big(1-cn^{-2}\big)^{k}O(n2)\mathcal{O}(n^{2})Roughly twice as fast as Jacobi
SOR (optimal ω\omega )(1−cn−1)k\big(1-cn^{-1}\big)^{k}O(n)\mathcal{O}(n)Improved, but still asymptotically slow
Steepest descent(κ−1κ+1)k\big(\frac{\kappa-1}{\kappa+1}\big)^{k}O(κ)\mathcal{O}(\kappa)Gradient method; zig-zag convergence
Conjugate Gradient2(κ−1κ+1)k2\Big(\tfrac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\Big)^{k}O(κ)\mathcal{O}(\sqrt{\kappa})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 AA is never factorised; it is only applied to vectors. After mm steps the problem has been compressed into an m×mm\times m symmetric tridiagonal matrix TmT_{m} whose eigenvalues (the Ritz values ) approximate the extreme eigenvalues of AA , and whose eigenvectors lift back to approximate eigenvectors of AA through QmQ_{m} .
Figure 4.5 : Left: the truncated SVD cuts small- σi\sigma_{i} 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 k⋆k^{\star} .
Figure 5.1 : The Arnoldi decomposition AVk=Vk+1H‾kAV_{k}=V_{k+1}\underline{H}_{k} , drawn for k=4k=4 . Applying AA to the orthonormal basis VkV_{k} of Kk\mathcal{K}_{k} never leaves Kk+1\mathcal{K}_{k+1} , so the result factors through exactly one extra basis vector. Column jj of H‾k\underline{H}_{k} is nonzero only in its first j+1j+1 entries, which is what upper Hessenberg means; in particular the first column has just two nonzeros. The matrix is (k+1)×k(k+1)\times k , the last row containing the single entry hk+1,kh_{k+1,k} that carries vk+1v_{k+1} . When AA is symmetric, H‾k\underline{H}_{k} 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 ∥q(A)∥2=max⁡j∣q(λj)∣\|q(A)\|_{2}=\max_{j}|q(\lambda_{j})| , 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 F(A)={x∗Ax:∥x∥2=1}F(A)=\{x^{\ast}Ax:\|x\|_{2}=1\} , a convex set containing σ(A)\sigma(A) . If F(A)F(A) stays in a disk that excludes the origin, the Crouzeix–Palencia inequality [ 8 ] gives a genuine convergence rate; if the origin lies inside F(A)F(A) , no such bound is available, and GMRES can stagnate for n−1n-1 steps even with all eigenvalues equal to 11 (see the example at the end of this section).
α\alpha∣λ2(G)∣|\lambda_{2}(G)|power iterations to 10−810^{-8}
0.500.500.3540.3541717
0.850.850.6010.6013636
0.950.950.6720.6724545
0.990.990.7000.7005252
Table 5.1 : The damping factor on the three-node example of this section, where ∣λ2(G)∣=α/2|\lambda_{2}(G)|=\alpha/\sqrt{2} exactly. The last column is the number of iterations observed to gain eight digits, in agreement with the asymptotic estimate 8log⁡10/log⁡(1/∣λ2(G)∣)8\log 10/\log(1/|\lambda_{2}(G)|) . The value α=0.85\alpha=0.85 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 ∥Ax−b∥2\|Ax-b\|_{2} is the point of the column space R(A)\mathcal{R}(A) closest to bb , so the residual is orthogonal to R(A)\mathcal{R}(A) , which leads exactly to the normal equations A⊤Ax=A⊤bA^{\top}Ax=A^{\top}b . LSQR realises the same projection without forming A⊤AA^{\top}A .
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 AA and one with A⊤A^{\top} . The coefficients assemble into a lower bidiagonal BkB_{k} with AVk=Uk+1BkAV_{k}=U_{k+1}B_{k} . This is Lanczos applied implicitly to A⊤AA^{\top}A 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 A⊤bA^{\top}b and grows by powers of A⊤AA^{\top}A , components associated with large singular values are resolved first. Stopping early therefore acts as a spectral filter, suppressing the small- σi\sigma_{i} 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 AA (left) are spread, whereas those of M−1AM^{-1}A (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: D1D_{1} (top) and D2D_{2} (bottom). On the overlap, D1=D2=12D_{1}=D_{2}=\tfrac{1}{2} and outside the overlap they sum to one: R1TD1R1+R2TD2R2=IR_{1}^{T}D_{1}R_{1}+R_{2}^{T}D_{2}R_{2}=I .
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 m=31m=31 ; the thick blue curve is the same error after five weighted-Jacobi sweeps ( ω=23\omega=\tfrac{2}{3} ); the dashed orange curve is the piecewise-linear interpolant through its values at the 1515 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 O(n)\mathcal{O}(n) work and reduces the error by a factor independent of hh . Used as a preconditioner for CG, this turns the O(h−1)\mathcal{O}(h^{-1}) iteration count of Chapter 4 into a constant.

Explore similar work

CardsList
Oct 3, 2026

The Numerical Linear Algebra of Large Language Models

Numerical Linear Algebra (NLA) has consistently played a vital role in advancing science by providing tools to solve fundamental problems encountered in scientific and engineering applications. Over the decades, it has continually evolved to meet the demands driven by successive waves of scientific discovery. For instance, during the 1950s and 1960s, substantial efforts were devoted to developing methods for solving eigenvalue problems that emerged from the rapidly growing field of aerodynamics. This led to the discovery of the LR and QR algorithms. Later the attention turned to the solution of sparse linear systems that were common in applications like computational aerodynamics. Today we are experiencing yet another wave of major scientific advancement and NLA is once more at the heart of its development. This Machine Learning (ML) wave is proving to be utterly disruptive in science and engineering. Many tools in ML particularly Large Language Models (LLMs) are grounded in matrix and tensor methods. As we are approaching Artificial General Intelligence (AGI), it is clear that matrix methods will be called to play an even more significant role. For the numerical linear practitioner the speed of the current change makes it particularly challenging to adapt. This is a survey article that centers on machine learning techniques, with a particular focus on large language models. It has two main objectives. The first is to clarify the core concepts behind Large Language Models in a manner accessible to specialists in numerical methods. The second is to examine the key Numerical Linear Algebra concepts employed by LLM techniques, while also highlighting several significant recent contributions of NLA to the field.

Abdelkader Baggag, Yousef Saad

Qatar Computing Research Institute, Hamad Bin Khalifa University, HBKU Research Complex, Doha, Qatar. · Dept. Computer Science and Engineering, University of Minnesota, 200 Union Street S.E., Minneapolis, MN, USA.

Topic pending
AskPDF
Feb 10, 2026

Learning to Discover Iterative Spectral Algorithms

We introduce AutoSpec, a neural network framework for discovering iterative spectral algorithms for large-scale numerical linear algebra and numerical optimization. Our self-supervised models adapt to input operators using coarse spectral information (e.g., eigenvalue estimates and residual norms), and predict recurrence coefficients for computing or applying a matrix polynomial tailored to a downstream task. The effectiveness of AutoSpec relies on three ingredients: an architecture whose inference pass implements short, executable numerical linear algebra recurrences; efficient training on small synthetic problems with transfer to large-scale real-world operators; and task-defined objectives that enforce the desired approximation or preconditioning behavior across the range of spectral profiles represented in the training set. We apply AutoSpec to discovering algorithms for representative tasks on spd matrices: accelerating matrix function approximation; accelerating sparse linear solvers; and spectral filtering/preconditioning for eigenvalue computations. On real-world matrices, the learned procedures deliver up to order-of-magnitude improvements in accuracy and/or reductions in iteration count, relative to spectrum-agnostic baselines. We find clear connections to classical theory: the induced polynomials may exhibit equioscillation behavior characteristic of Chebyshev polynomial approximation. The code is available at: https://github.com/zihanghliu/AutoSpec .

Zihang Liu, Oleg Balabanov, Yaoqing Yang +1

1International Computer Science Institute · University of California, Berkeley · 4Lawrence Berkeley National Laboratory +1

Automated Algorithm DiscoveryLearning to OptimizeSpectral Methods
AskPDF
Apr 25, 2026

Well-Conditioned Oblivious Perturbations in Linear Space

Perturbing a deterministic nn-dimensional matrix with small Gaussian noise is a cornerstone of smoothed analysis of algorithms [Spielman and Teng, JACM 2004], as it reduces the condition number of the input to O(n)O(n), and with it the complexity of many matrix algorithms. However, when deployed algorithmically, these perturbations are expensive due to the cost of generating and storing n2n^2 Gaussian random variables. We propose a perturbation that requires generating and storing O(n)O(n) random numbers in O(log⁡n)O(\log n) bits of precision, and reduces the condition number of any deterministic matrix to O(n)O(n), matching Gaussian perturbations. Our result in particular implies a better complexity for the perturbed conjugate gradient algorithm, showing that we can solve an n×nn\times n linear system in linear space to within an arbitrarily small constant backward error using O(n)O(n) matrix-vector products. In our construction, we introduce the concept of a pattern matrix, which is a dense deterministic matrix that maps all sparse vectors into dense vectors, and we combine it with a sparse perturbation whose entries are dependent and located in a non-uniform fashion. In order to analyze this construction, we develop new techniques for lower bounding the smallest singular value of a random matrix with dependent entries.

Shabarish Chenakkod, Michał Dereziński, Xiaoyu Dong +1

University of Michigan, Ann Arbor, MI, USA · National University of Singapore, Singapore

Random Matrix Theory
AskPDF