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

arXiv:2608.21234v1 [math.NA] 21 Aug 2026

Advanced Linear Algebra with Applications Part I

PDE discretisation

Sparse matrices

Graph Laplacians

Victorita Dolean

Spectral analysis in ML

& Jemima Tabeart

Numerical linear algebra for PDEs, machine learning, and data assimilation

Lecture notes August 24, 2026

ii

Contents Preface

v

1 Solving Linear Systems: A Unified Foundation 1.1 Norms and Matrix Decompositions . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Vector and Matrix Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Eigenvalues and Matrix Decompositions . . . . . . . . . . . . . . . . . . . 1.1.3 Matrix Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Floating-Point Arithmetic, Error, Conditioning, and Factorizations . . . . . . . . 1.2.1 Floating-Point Arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Error, Conditioning and Stability of algorithms . . . . . . . . . . . . . . . 1.2.3 Gaussian Elimination and LU Factorization . . . . . . . . . . . . . . . . . 1.2.4 Cholesky Factorization for SPD Matrices . . . . . . . . . . . . . . . . . . . 1.3 Stationary Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Examples of Stationary Methods and their convergence . . . . . . . . . . 1.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 2 2 3 3 4 4 6 7 8 8 10 13

2 Discretization of 1D and 2D PDEs (and Beyond) 2.1 Finite difference approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Sparse Matrices: PDEs, Graphs, and Machine Learning . . . . . . . . . . . . . . 2.2.1 The 2D Poisson equation: structure, stencils and blocks . . . . . . . . . . 2.2.2 Graph Laplacians . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Sparse Matrices in Machine Learning . . . . . . . . . . . . . . . . . . . . . 2.2.4 Sparsity in data assimilation . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Sparse Matrix Toolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Sparsity Patterns and storage formats . . . . . . . . . . . . . . . . . . . . 2.3.2 Spectral Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15 15 18 18 19 21 22 22 22 25 28

3 Applications of Basic Iterative Methods 3.1 Basic Iterative Methods for PDE discretisations . . . . . . . . . . . . . . . . . . . 3.1.1 Convergence of the Jacobi iteration . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Convergence of the Gauss–Seidel iteration . . . . . . . . . . . . . . . . . . 3.1.3 Convergence of the SOR iteration . . . . . . . . . . . . . . . . . . . . . . . 3.2 Basic Iterative Methods for Graph Laplacians . . . . . . . . . . . . . . . . . . . . 3.3 Basic Iterative Methods in Machine Learning . . . . . . . . . . . . . . . . . . . . 3.3.1 Richardson iteration and gradient descent . . . . . . . . . . . . . . . . . . 3.3.2 Applications in machine learning . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Smoothing and spectral bias . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31 31 32 32 33 34 39 39 40 40 42

iii

iv

CONTENTS

4 Krylov Methods for Symmetric Systems 4.1 The Conjugate Gradient method . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The Lanczos Method and the Symmetric Eigenproblem . . . . . . . . . . . . . . . 4.2.1 Lanczos and CG are the same recurrence . . . . . . . . . . . . . . . . . . . 4.2.2 Application: computing the Fiedler vector . . . . . . . . . . . . . . . . . . 4.3 Krylov Methods and Machine Learning (optional) . . . . . . . . . . . . . . . . . . 4.3.1 Least squares and the SVD . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Spectral filters and the truncated SVD . . . . . . . . . . . . . . . . . . . . 4.3.3 Ridge regression: a smooth filter . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Early stopping: the filter generated by CG . . . . . . . . . . . . . . . . . . 4.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45 46 50 53 54 55 55 55 56 56 58

5 Krylov Methods for Non-Symmetric Systems 5.1 GMRES and Arnoldi for non-normal matrices . . . . . . . . . . . . . . . . . . . . 5.1.1 GMRES: minimising the residual over the Krylov space . . . . . . . . . . 5.1.2 Convergence: why eigenvalues are not enough . . . . . . . . . . . . . . . . 5.2 Networks: PageRank as an Eigenproblem . . . . . . . . . . . . . . . . . . . . . . 5.2.1 The role of the damping factor . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Beyond the power iteration: the Arnoldi process . . . . . . . . . . . . . . 5.3 Krylov Methods for Large Least Squares (optional) . . . . . . . . . . . . . . . . . 5.3.1 Golub–Kahan bidiagonalisation and LSQR . . . . . . . . . . . . . . . . . . 5.3.2 Early stopping and the truncated SVD . . . . . . . . . . . . . . . . . . . . 5.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61 62 63 65 66 68 68 69 71 72 74

6 Preconditioning and Accelerating Solvers 6.1 Motivation and Basic Preconditioners . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Preconditioning SPD problems . . . . . . . . . . . . . . . . . . . . . . . . 6.1.2 A general framework for preconditioners . . . . . . . . . . . . . . . . . . . 6.1.3 One framework for the whole chapter . . . . . . . . . . . . . . . . . . . . . 6.2 Schwarz Domain Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Two-Level Preconditioning and Multigrid Ideas . . . . . . . . . . . . . . . . . . . 6.3.1 Weighted Jacobi and the smoothing factor . . . . . . . . . . . . . . . . . . 6.3.2 Coarse-grid correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.3 From coarse correction to the two-grid method . . . . . . . . . . . . . . . 6.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77 78 78 79 80 80 84 85 87 88 90

Bibliography

93

Preface These notes form the first part of a longer course on advanced numerical linear algebra, taught at master’s level. Their purpose 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 was where its large sparse systems came from. This is no longer the whole story. Ranking the nodes of a network, assimilating observations into a weather forecast, and fitting a model to a large noisy data set all lead to problems of exactly the same kind: too large to factorise, structured, and accessible only through matrix–vector products. What is striking is how few ideas are needed for all of them. Sparsity, the spectrum, Krylov subspaces and preconditioning recur throughout, and a method designed for a discretised Laplacian is usually the right method for a graph Laplacian or a design matrix as well. The chapters that follow are organised around that observation. Each develops a classical topic and then puts it to work outside its original setting: sparse matrices from finite differences, graphs and machine learning; stationary iterations and the smoothing property; conjugate gradient and Lanczos, with spectral clustering and regularisation by early stopping; Arnoldi and GMRES, with PageRank and large least squares; and finally preconditioning, domain decomposition and multigrid. This first part is deliberately foundational. The material collected here (factorisations and conditioning, iterative methods, Krylov subspaces, preconditioners) is what the second part of the course, not reproduced in this volume, takes for granted throughout. The level is that of a first-year master’s student in mathematics, computer science or a computational discipline. We assume a first course in linear algebra and some familiarity with numerical analysis, and we prove what can be proved compactly, leaving the more technical results to the literature. Every section closes with a summary of what should be retained, and every chapter with exercises, several taken from past examinations. Reading about an iterative method is, however, a poor substitute for watching one converge. The illustrations in these notes are reproducible, and we encourage readers to run and modify the accompanying Python code, available at https://github.com/vicdolean/scicomp_examples. Changing a parameter, a matrix or a stopping criterion and watching what happens to the convergence history is by far the fastest way to develop a feeling for when these methods work and when they do not. A note on the use of AI tools. The structure of these notes, the choice and ordering of the topics, and the mathematical content are the authors’ own. Some of the material is standard and presented here in our own way; some of it is drawn from research we are ourselves engaged in. In preparing the manuscript we used a generative AI assistant (Claude Opus 5, Anthropic) to help identify the seminal references on a given topic, to improve the language, and to improve the graphical content. All statements, proofs and references have been checked by the authors, who take full responsibility for the contents. Victorita Dolean and Jemima M. Tabeart Eindhoven, August 24, 2026 v

vi

Preface

Chapter 1

Solving Linear Systems: A Unified Foundation

Overview

Linear algebra forms the backbone of many modern scientific and engineering computations. From solving systems of equations in structural analysis, optimizing models in machine learning, and forecasting weather through data assimilation, linear algebraic methods are pervasive. However, in practical computations, the mathematical representations of real numbers must confront the limitations of computer arithmetic. Floating-point numbers can only approximate real numbers, leading to small but inevitable rounding errors, and even mathematically exact operations can be extremely sensitive to those tiny errors depending on the structure of A.

This chapter lays the numerical foundations used throughout the rest of the course: how to measure the size of vectors, matrices, and errors (norms); how a problem’s own sensitivity to perturbations is quantified (conditioning); how floating-point arithmetic and algorithms can amplify or control error (stability); and the two standard families of solvers for Ax = b — direct factorizations (LU, Cholesky) and stationary iterative methods (Jacobi, Gauss–Seidel, SOR, Richardson). The Scale of Linear Algebra in Real Life Weather forecasts use systems with 109 unknowns – solving linear systems is essential every few hours, and must be done in minutes. Google’s PageRank was originally computed by solving a giant eigenvalue problem. Deep Learning has a lot of linear algebra: backpropagation is mostly matrix multiplies and chain-rule gradients. 1

2

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION Learning objectives By the end of this chapter, you should be able to: • Define and compute common vector and matrix norms, and explain why all norms are equivalent in finite dimension. • Relate eigenvalues, eigendecompositions, the Schur decomposition, and the SVD to the spectral radius and to matrix norms. • Explain the sources of floating-point error (rounding, truncation, cancellation) and use the condition number κ(A) to relate backward and forward error. • Compute the LU and Cholesky factorizations of a matrix, and state why partial pivoting is needed for backward stability. • Derive the Jacobi, Gauss–Seidel, SOR, and Richardson iterations from a matrix splitting A = M − N , and state the convergence criterion ρ(G) < 1.

1.1

Norms and Matrix Decompositions

The first electronic computer, ENIAC (1940s), already faced severe rounding problems. Engineers had to carefully redesign algorithms to avoid divergence even after a few thousand operations.

1.1.1

Vector and Matrix Norms

Definition 1.1.1 (Vector norms). Given a vector x ∈ Rn , the common norms are: P • ℓ1 norm: ∥x∥1 = ni=1 |xi | • ℓ2 norm (Euclidean): ∥x∥2 =

2 1/2 i=1 |xi |

Pn

• ℓ∞ norm: ∥x∥∞ = maxi |xi | These norms satisfy the following properties: • Non-negativity: ∥x∥ ≥ 0 and ∥x∥ = 0 if and only if x = 0 • Homogeneity: ∥αx∥ = |α|∥x∥ for any scalar α • Triangle inequality: ∥x + y∥ ≤ ∥x∥ + ∥y∥ Norms give us a way to measure the size (or length) of a vector. They are particularly useful in determining convergence of sequences, measuring distances between vectors, and quantifying approximation errors. Proposition 1.1.2 (Norm equivalence). All norms on a finite-dimensional vector space are equivalent. This means that for any two norms ∥ · ∥a and ∥ · ∥b , there exist constants c1 , c2 > 0 such that: c1 ∥x∥a ≤ ∥x∥b ≤ c2 ∥x∥a for all x ∈ Rn . Thus, different norms can be used interchangeably when analyzing asymptotic behavior, even √ though their numerical values may differ. For instance, ∥x∥∞ ≤ ∥x∥2 ≤ n∥x∥∞ and √ √ ∥x∥2 ≤ ∥x∥1 ≤ n∥x∥2 , so here c1 = 1 and c2 = n (see Problem 2).

1.1. NORMS AND MATRIX DECOMPOSITIONS

3

Figure 1.1: Unit balls {x : ∥x∥ = 1} for the ℓ1 , ℓ2 , and ℓ∞ norms in R2 : the same vector has a different "size" depending on which norm is used.

1.1.2

Eigenvalues and Matrix Decompositions

Before turning to matrix norms, it is useful to recall the key spectral objects that will reappear throughout the course: eigenvalues, and the two canonical matrix factorizations built from them. Definition 1.1.3 (Eigenvalues, eigendecomposition, spectral radius). A scalar λ ∈ C and nonzero vector x satisfying Ax = λx are an eigenvalue and eigenvector of A. If A ∈ Rn×n is diagonalizable, it admits an eigendecomposition A = V ΛV −1 ,

Λ = diag(λ1 , . . . , λn ),

where the columns of V are eigenvectors. The spectral radius is ρ(A) = maxi |λi |. Definition 1.1.4 (Schur decomposition and SVD). Every A ∈ Cn×n admits a Schur decomposition A = QT Q∗ with Q unitary (Q∗ Q = QQ∗ = I) and T upper triangular, with the eigenvalues of A on the diagonal of T . Every A ∈ Rm×n (or Cm×n ) admits a singular value decomposition (SVD) A = U ΣV ∗ with U, V unitary and Σ = diag(σ1 , . . . ) where σi ≥ 0. Remark 1.1.5. Every symmetric matrix is diagonalizable with real eigenvalues; if it is also symmetric positive definite (SPD), i.e. x⊤ Ax ≥ 0 with equality only for x = 0, its eigenvalues are strictly positive. Unlike the eigendecomposition, the Schur decomposition and the SVD exist for every matrix, which is why they underpin general-purpose numerical algorithms. Takeaways Eigenvalues and the spectral radius ρ(A) control the long-run behaviour of iterative processes (Section 1.3); the Schur decomposition and the SVD are the two factorizations that always exist and that later chapters (and the induced 2-norm below) build on.

1.1.3

Matrix Norms

Matrix norms are critical for analyzing the stability and conditioning of algorithms and problems involving linear systems, eigenvalue computations, and optimization. Matrix norms are often induced by vector norms.

4

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION Definition 1.1.6 (Induced norm (operator norm)). For a matrix A ∈ Rm×n the induced norm (operator norm) is defined by ∥Ax∥ = max ∥Ax∥ ∥x∥=1 x̸=0 ∥x∥

∥A∥ = sup

The three induced norms used throughout these notes are the following. • Spectral norm (2-norm): ∥A∥2 = σmax (A), where σmax is the largest singular value. • 1-norm ∥A∥1 : maximum absolute column sum • Infinity-norm ∥A∥∞ : maximum absolute row sum Not every matrix norm is induced by a vector norm. The most important exception is the Frobenius norm,  1/2 m X n X ∥A∥F =  |aij |2  i=1 j=1

It is equivalent to the ℓ2 norm of the matrix viewed as a vector of its entries. Proposition 1.1.7 (Note on induced norms). For an induced matrix norm ∥A∥, we have: • ∥AB∥ ≤ ∥A∥ · ∥B∥ • ∥A∥ ≥ ∥Ax∥/∥x∥ for all x ̸= 0 This submultiplicative property makes induced norms useful in bounding errors and proving convergence. Takeaways On finite-dimensional spaces all norms are equivalent: asymptotically they give the same notion of convergence, but numerically their values can differ significantly.

1.2

Floating-Point Arithmetic, Error, Conditioning, and Factorizations

1.2.1

Floating-Point Arithmetic

The Ariane 5 rocket self-destructed (1996) due to a floating-point overflow error when converting a 64-bit number to a 16-bit integer. This incident underscores the risks inherent in floating-point arithmetic and careless type conversions. Definition 1.2.1 (Floating-point model). A real number x is represented as: x ≈ fl(x) = x(1 + δ),

|δ| ≤ εmach

where εmach is the machine epsilon. Floating-point arithmetic is governed by the IEEE 754 standard, which defines how numbers are stored and manipulated. A typical floating-point number consists of: • A sign bit (1 bit)

1.2. FLOATING-POINT ARITHMETIC, ERROR, CONDITIONING, AND FACTORIZATIONS5 • An exponent field (used to scale the number) • A significand or mantissa (which holds the significant digits) The first IEEE standard was introduced in 1985. For a long time single (8 exponent bits and 23-bits for the significand) and double precision (11 exponent bits and 52 significand bits) were the most commonly-used formats. However, in recent years there has been a growth in interest in reduced precisions, such as half precision (5 exponent bits and 10 significand bits), non-standard precisions such as bfloat16 (8 exponent bits and 7 significand bits) and proprietary formats like Nvidia’s TensorFloat-32 (8 exponent bits and 10 significand bits). Much of this development has been driven by machine learning applications, and the desire to accelerate performance. In data assimilation, the computational cost is dominated by evaluations of the PDE model, so large gains can be made by applying the model in single precision [19]. One option is to run an entire algorithm in a reduced precision, which will lead to increased computational efficiency, but may result in degraded performance. This has motivated ‘mixedprecision’ approaches, particularly in the numerical linear algebra community, where different steps of an algorithm are run in different precisions. Such approaches aim to combine the computational benefits of reduced precision, while ensuring rapid convergence and computational stability (see e.g. [22, 1] for comprehensive reviews). Machine Epsilon The machine epsilon εmach is the smallest number such that: 1 + εmach ̸= 1 It gives an upper bound on relative error due to rounding in floating-point arithmetic. For double precision, εmach ≈ 2−53 ≈ 1.11 × 10−16 . Sources of Numerical Errors • Rounding errors: Due to truncating or rounding numbers to fit into finite precision. • Truncation errors: Arise when series are approximated by finite sums. • Cancellation errors: Occur when subtracting two nearly equal numbers, leading to loss of significant digits. Example: Catastrophic Cancellation Suppose we compute x = a − b with a = 1.0000001 and b = 1.0000000. In exact arithmetic: x = 0.0000001 but in floating-point arithmetic, a and b may be stored with limited precision, and the result may be inaccurate or even zero. This is an example of catastrophic cancellation, where most significant digits cancel out. Takeaways • εmach (machine epsilon) quantifies roundoff. For double precision: ≈ 10−16 . • Cancellation (subtracting nearly equal numbers) is the most dangerous source of error. • Floating-point arithmetic is not associative: (x + y) + z ̸= x + (y + z). • Any finite precision implementation can be affected by overflow/rounding errors, but lower precisions are more at risk.

6

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION

1.2.2

Error, Conditioning and Stability of algorithms

Numerical errors in computations arise from various sources and can significantly affect the reliability of results. To understand and mitigate these effects, we analyze the nature of errors and how sensitive a problem is to perturbations. During the Gulf War (1991), the Patriot missile system failed to intercept a missile due to a small rounding error in the system’s clock, leading to a timing drift of 0.34 seconds. This small error had devastating consequences, highlighting the importance of careful error analysis.

The Norwegian oil platform Sleipner A collapsed (1991) due to a small error in the finite element analysis stiffness matrix, demonstrating how critical numerical stability and accuracy are in engineering simulations. Definition 1.2.2 (Forward and backward error). If x is the exact solution and x̂ the computed one: forward error = ∥x̂ − x∥,

backward error = min{∥δA∥ : (A + δA)x̂ = b}.

Definition 1.2.3 (Condition number). For a nonsingular A, κ(A) = ∥A∥ · ∥A−1 ∥. Definition 1.2.4 (Condition number in the 2-norm). For a symmetric positive definite A, κ2 (A) = ∥A∥2 · ∥A−1 ∥2 =

λmax (A) . λmin (A)

Proposition 1.2.5 (Error amplification). A perturbation δb yields ∥δx∥ ∥δb∥ ≤ κ(A) . ∥x∥ ∥b∥ Example 1.2.6 (Error amplification). Consider   1 1 A= , 1 1.0001



 2 b= . 2.0001

This system is very close to being singular, so solving Ax = b produces a result that is highly sensitive to small perturbations in A or in b. Stability at a glance • Stable algorithm: Controls the propagation of rounding errors. • Backward stable algorithm: Solves a nearby (perturbed) problem exactly. • Combined view: Stability of the algorithm + conditioning of the problem together determine the final accuracy.

1.2. FLOATING-POINT ARITHMETIC, ERROR, CONDITIONING, AND FACTORIZATIONS7 Takeaways • Small κ(A) (close to 1): well-conditioned, stable predictions. • Large κ(A): ill-conditioned, even perfect algorithms may give unreliable results. • An algorithm can be stable (small backward error) yet produce inaccurate results if the problem is ill-conditioned (large condition number). • Conversely, for ill-conditioned problems no algorithm can suppress error amplification.

1.2.3

Gaussian Elimination and LU Factorization

Gaussian elimination, though often attributed to Gauss (1820), was already known in ancient China (Chapter 8 of The Nine Chapters on the Mathematical Art, ∼200 AD). Gauss’ major contribution was to interpret elimination as a matrix factorization. It took him more than twenty years to arrive at that interpretation. We wish to solve Ax = b with A ∈ Rn×n , and proceed in three steps. 1. Apply elimination to zero out entries below the diagonal. 2. This produces a factorization A = LU with L unit lower triangular, U upper triangular. 3. Solve Ly = b by forward substitution, then U x = y by backward substitution. Definition 1.2.7 (LU factorization). A matrix A ∈ Rn×n admits an LU factorization if it can be written A = LU, where L is unit lower triangular and U is upper triangular. Example 1.2.8 (Worked 3 × 3 LU factorization, no pivoting). For   2 1 1 A =  4 −6 0 , −2 7 2 elimination proceeds as: ℓ21 = 4/2 = 2, row2 ←row2 − 2 , row1 = [0, −8, −2]; ℓ31 = −2/2 = −1, row3 ←row3 +, row1 = [0, 8, 3]; pivot U22 = −8, ℓ32 = 8/(−8) = −1, row3 ←row3 − (−1) row2 = [0, 0, 1]. This gives     1 0 0 2 1 1 1 0 , L= 2 U = 0 −8 −2 , −1 −1 1 0 0 1 and one checks A = LU . To solve Ax = b, solve Ly = b by forward substitution, then U x = y by backward substitution. Proposition 1.2.9 (Complexity). The cost of Gaussian elimination (LU factorization) is O(n3 ) flops, while forward and backward substitution cost O(n2 ) each. Remark 1.2.10. For large linear systems with n ∼ 106 , O(n3 ) is infeasible. This motivates the use of iterative solvers and preconditioners. Definition 1.2.11 (Partial pivoting). During elimination, rows are swapped to ensure that pivot elements are large in magnitude. This produces a factorization P A = LU,

P a permutation matrix.

8

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION

Theorem 1.2.12 (Backward stability of GEPP). Gaussian elimination with partial pivoting (GEPP) is backward stable: the computed x̂ solves (A + δA)x̂ = b,

∥δA∥ = O(εmach ). ∥A∥

Pitfall Without pivoting, Gaussian elimination may fail catastrophically even for well-conditioned A. [5] develops bounds on the growth of entries when using full pivoting (not commonly used in practice). The appendix contains a nice overview of computational issues with partial pivoting, including in relation to new hardware such as GPUs. Gauss used elimination in 1809 to compute planetary orbits from astronomical data. Today, LU is the default “black box” in numerical linear algebra libraries (LAPACK, and hence NumPy/SciPy). Takeaways • LU factorization provides a direct, stable way to solve Ax = b. • Pivoting is essential for stability. • Complexity O(n3 ) makes LU unsuitable for very large sparse linear systems. • Iterative solvers (next section) are designed to overcome this scaling barrier.

1.2.4

Cholesky Factorization for SPD Matrices

When A is symmetric positive definite, LU factorization can be specialized to exploit the symmetry, halving both storage and computational cost. Definition 1.2.13 (Cholesky factorization). If A ≻ 0 is SPD, then A = LL⊤ for a unique lower triangular L with positive diagonal entries, computed via v P u k−1 X u aik − k−1 j=1 ℓij ℓkj 2 t , i > k. ℓkk = akk − ℓkj , ℓik = ℓkk j=1

Key idea Cholesky needs no pivoting and is backward stable for SPD matrices, uses half the storage of LU (only L is stored), and costs ≈ n3 /3 flops — about half the cost of a general LU factorization. Takeaways If A is SPD, use Cholesky first: it is the fastest, most memory-efficient, backward-stable direct solver available. If some ℓkk ≤ 0 is encountered during the factorization, A is not SPD (or is poorly scaled); switch to LDL⊤ with symmetric pivoting, or regularize via A ← A + λI.

1.3

Stationary Iterative Methods

Direct methods (Gaussian elimination / LU with pivoting) are robust and backward stable, but they become impractical at scale.

9

1.3. STATIONARY ITERATIVE METHODS

• Cost. Factorization costs O(n3 ) flops; each triangular solve costs O(n2 ). For sparse banded matrices with half-bandwidth w, the costs are O(nw2 ) (factor) and O(nw) (solve). • Memory. Fill-in during elimination destroys sparsity; storage can grow from O(nnz(A))1 to orders of magnitude larger. Key idea Iterative methods trade exactness for scalability: they approximate the solution via many cheap updates. One step costs O(nnz(A)), preserves sparsity, and is parallel-friendly.

In the early 19th century, Gauss and Jacobi independently realized that, instead of eliminating variables, one can successively refine a guess. Jacobi’s 1845 paper on planetary orbits is often cited as the first systematic study of iterative solvers. Definition 1.3.1 (Splitting method). Let A = M − N with M invertible. Define the stationary iteration x(k+1) = M −1 N x(k) + M −1 b,

G := M −1 N,

c := M −1 b.

If it converges, the limit x⋆ satisfies (I − G)x⋆ = c and hence Ax⋆ = b. Theorem 1.3.2 (Convergence criterion). The iteration converges for every initial guess x(0) if and only if ρ(G) < 1. Moreover, x(k) → x⋆ = (I − G)−1 c and the solution is unique. Sketch. Let e(k) := x(k) − x⋆ . Subtract the fixed-point relation (I − G)x⋆ = c from the iteration to obtain e(k+1) = G e(k) , hence e(k) = Gk e(0) . (⇒) If ρ(G) ≥ 1, there is an eigenpair (λ, v) with |λ| ≥ 1; choose e(0) = v to get e(k) = λk v, which does not tend to 0. Thus convergence for all x(0) implies ρ(G) < 1. (⇐) If ρ(G) < 1, then for any ε > 0 there exists a (subordinate) norm ∥ · ∥ε with ∥G∥ε ≤ 2 . Therefore ∥e(k) ∥ ≤ ∥G∥k ∥e(0) ∥ ≤ q k ∥e(0) ∥ → 0. Equivalently, the ρ(G) + ε =: q < ε ε ε P1 Pε ∞ j , hence x⋆ = (I − G)−1 c. Neumann series j=0 Gj converges3 and (I − G)−1 = ∞ G j=0 Corollary 1.3.3 (A practical sufficient condition). If there exists any subordinate matrix norm with ∥G∥ < 1, then x(k) → x⋆ and ∥e(k) ∥ ≤ ∥G∥k ∥e(0) ∥,

∥x⋆ − x(k) ∥ ≤

∥G∥k ∥x(1) − x(0) ∥. 1 − ∥G∥

Takeaways To design an effective stationary method: • Choose M easy to invert (triangular / block-diagonal / diagonal) so each step is cheap. • Ensure ρ(M −1 N ) < 1 (often via diagonal dominance or SPD structure). • Prefer splittings that keep ∥G∥ comfortably below 1 to avoid roundoff-limited plateaus. 1

nnz(A) denotes the number of non-zero elements of the matrix A see the theoretical problems at the end of the chapter 3 we admit this result without proof 2

10

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION

1.3.1

Examples of Stationary Methods and their convergence

Write A = D − L − U with L strictly lower, D diagonal, U strictly upper triangular. • Jacobi: M = D, GJ = D−1 (L + U ). Each component is updated using only values from the previous iterate. • Gauss–Seidel: M = D − L, GGS = (D − L)−1 U . New values are used immediately within the same sweep, leading to faster convergence. • Successive Over-Relaxation (SOR): M = ω1 D + L,

Gω = ( ω1 D − L)−1



1−ω ω D−U



.

Here ω > 0 is a relaxation parameter. For ω = 1 we recover Gauss–Seidel. • Richardson’s method: M = α1 I, GR = I − αA with α a well-chosen parameter. Remark 1.3.4. While Jacobi updates all components simultaneously from the old iterate, Gauss– Seidel reuses fresh values within the same sweep. SOR goes one step further: each update is extrapolated by ω, effectively accelerating convergence if ω is chosen well. Lewis Fry Richardson also proposed the first numerical weather forecasting approach, via the solution of differential equations. His paper was published in 1922 and is the foundation of modern numerical weather prediction. He attempted to perform computations by hand while serving with a Quaker ambulance unit on the Western Front in 1917 (his results were not good, but the method was later shown to be sound https://www.americanscientist.org/article/theweatherman). The SOR method was proposed by David M. Young in 1950 while working on elliptic PDEs as a PhD student at Harvard. His method was so effective that it immediately became the standard iterative scheme of its era, especially in engineering applications. Theorem 1.3.5 (Convergence under diagonal dominance). If A is strictly diagonally dominant by rows, then both Jacobi and Gauss–Seidel converge. Sketch. Write A = L + D + U the for Jacobi method we have GJ = D−1 (L + U ). For the ∥ · ∥∞ norm, P X aij j̸=i |aij | ∥GJ ∥∞ = max = max . i i aii |aii | j̸=i

P Strict row diagonal dominance (SRDD) means |aii | > j̸=i |aij | for all i, hence ∥GJ ∥∞ < 1. Therefore ρ(GJ ) ≤ ∥GJ ∥∞ < 1, and Jacobi converges by the fixed-point criterion. For Gauss-Seidel similar arguments can be employed. An idea would be to use a weighted ∞-norm ∥x∥∞,w = maxi |xi |/wi . SRDD implies the existence of w > 0 with ∥GGS ∥∞,w < 1. Theorem 1.3.6 (Gauss–Seidel for SPD matrices, Theorem 11.2.3 [12]). If A is symmetric positive definite, then Gauss–Seidel converges. Remark 1.3.7. Jacobi and GS are not comparable in general: examples exist where one converges and the other does not.

11

1.3. STATIONARY ITERATIVE METHODS Theorem 1.3.8 (Convergence SOR (Kahan)). For any A and ω ∈ R, ρ(Gω ) ≥ | ω − 1 |. Hence a necessary condition for convergence is 0 < ω < 2. Proof. See [3, Theorem 6.32]

Theorem 1.3.9 (Convergence SOR (Ostrowski [28])). Let A a SPD matrix with the standard splitting A = D + L + L⊤ such that the diagonal D is also positive. Then SOR converges for all 0 < ω < 2. Theorem 1.3.10 (Convergence rate SOR (Young [42])). Assume that the Jacobi eigenvalues µ(GJ ) are real and that ρ(GJ ) < 1. Then the optimal relaxation parameter in SOR method is ω⋆ =

2 1+

p

1 − ρ(GJ

)2

(0 < ω ⋆ < 2),

and the optimized SOR factor is ρ(Gω⋆ ) =

ρ(GJ ) p 1 + 1 − ρ(GJ )2

!2 .

Theorem 1.3.11 (Convergence of Richardson’s method for SPD matrices [41]). Let A a SPD matrix. Then: • Convergence iff 0 < α <

2 2 = . ρ(A) λmax (A)

• The optimal step is α⋆ =

2 . λmax (A) + λmin (A)

• The optimized factor is ρ(I − α⋆ A) =

κ(A) − 1 with κ(A) = λmax /λmin . κ(A) + 1

Example 1.3.12 (Diagonally dominant system). The tridiagonal matrix   4 −1 0 · · ·   −1 4 −1 . . .    A=   0 −1 4 . . .    .. .. .. .. . . . . is strictly diagonally dominant. Jacobi and Gauss–Seidel are guaranteed to converge. Example 1.3.13 (Ill-behaved system). The Hilbert matrix Hij = 1/(i + j − 1) is symmetric positive definite, but extremely ill-conditioned. Gauss–Seidel converges in theory, but in finite precision it stagnates quickly due to roundoff. It is now widely recognised that Jacobi and Gauss–Seidel should not be judged solely by their speed as solvers. Their true importance lies in their role as building blocks for the preconditioned methods widely used in modern scientific computing.

12

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION

Takeaways • Existence of a good splitting and ρ(G) < 1 are the universal convergence keys. • Jacobi/GS: guaranteed under SRDD; GS also for all SPD A. • SOR: must take 0 < ω < 2; for SPD, convergence holds for every such ω; if Jacobi’s ρ(GJ ) is known/estimated, Young’s formula gives a near-optimal ω. • Richardson: simple baseline; optimal α⋆ gives factor (κ − 1)/(κ + 1). • Though obsolete as standalone solvers for large systems, Jacobi, GS, and SOR remain indispensable as smoothers in multigrid and as preconditioners in Krylov methods.

13

1.4. EXERCISES

1.4

Exercises

Problem 1 (Condition number comparison). Compute the condition number of • the identity matrix I, • a diagonal matrix with entries ranging from 1 to 105 , • a 5 × 5 Hilbert matrix. Discuss what these condition numbers tell you. Problem 2 (Vector norm equivalence). Prove that for any vector x ∈ Rn , √ ∥x∥∞ ≤ ∥x∥2 ≤ n∥x∥∞ , in other words that these norms are equivalent. Problem 3 (Properties of the Frobenius norm, exam 2024/2025). Let A ∈ Rn×n and let U, V ∈ Rn×n be unitary (orthogonal) matrices, i.e., U ⊤ U = U U ⊤ = I and V ⊤ V = V V ⊤ = I. The Frobenius norm is !1/2 n X n X a2ij . ∥A∥F := i=1 j=1

The trace of A is trace(A) =

Pn

i=1 aii .

(a) Show that the Frobenius norm is a matrix norm, but it is not induced by any vector norm. (b) Show that ∥A∥2F = trace(A⊤ A). (c) Show that unitary transformations leave the Frobenius norm unchanged: ∥U AV ∥F = ∥A∥F . Problem 4 (Induced 2-norm, resit exam 2024/2025). (a) Let A be an n × n matrix. Define the induced 2-norm as ∥A∥2 = sup ∥Ax∥2 . ∥x∥2 =1

Prove that the induced 2-norm satisfies the properties of a matrix norm. (b) Show that the induced 2-norm of a diagonal matrix is equal to the maximum absolute value of its diagonal entries. (c) Let A ∈ Rn×n and U, V be unitary matrices of appropriate dimensions. Prove that ∥U AV ∥2 = ∥A∥2 . Problem 5 (Catastrophic cancellation). Let a = 1.0000000001, b = 1.0000000 and assume we work in double precision. • Compute a − b and compare to the expected result. • Explain what happens if these numbers are rounded to fewer significant digits (single precision). Problem 6 (Floating-point summation order). Design an experiment to show that the order of operations in floating-point summation matters. Sum a list of numbers in increasing and decreasing order and compare the results.

14

CHAPTER 1. SOLVING LINEAR SYSTEMS: A UNIFIED FOUNDATION

Problem 7 (Effect of conditioning). Solve the systems Ax = b for two matrices:       1 2 1 1 5 A1 = , A2 = , b= 3 4 1 1.0001 11 Compare the computed solutions and relate the error to the condition number of A. Problem 8 (Backward error estimation). For a given A and x̂ that solves Ax ≈ b, compute r = b − Ax̂ and interpret r as a backward error. What does this tell you about the quality of x̂? Problem 9 (Approximating the spectral radius by an induced norm). For any A ∈ Cn×n and any ε > 0, show there exists a (vector) norm on Cn such that the induced matrix norm satisfies ρ(A) ≤ ∥A∥ ≤ ρ(A) + ε. Problem 10 (Gelfand’s formula for induced norms). For any induced matrix norm and any G ∈ Rn×n , prove lim ∥Gk ∥1/k = ρ(G). k→∞

Problem 11 (Neumann series). Let B ∈ Cn×n with ρ(B) < 1. Prove that I − B is invertible and ∞ X −1 (I − B) = Bj . j=0

Problem 12 (Basic iterative methods, resit exam 2024/2025). Consider the linear system Ax = b, where     4 −1 0 2    A = −1 4 −1 , b = 6 . 0 −1 4 2 (a) Derive the iteration formula for the Jacobi method for solving this system. (b) Perform two iterations of the Jacobi method starting from x(0) = [0, 0, 0]T . (c) Prove that the iteration matrix G for the Jacobi method satisfies ρ(G) < 1, ensuring convergence of the method. Problem 13 (Stationary methods applied to a SOAR correlation matrix). A secondorder autoregressive (SOAR) correlation matrix, as used for spatial correlations in data assimilation, is defined in Chapter 2. Implement Jacobi, Gauss–Seidel, and SOR for a SOAR matrix with n = 100, L = 0.4, a = 1 (see the course GitHub repository for a Python notebook). How do the numerical results connect to the theory?

Chapter 2

Discretization of 1D and 2D PDEs (and Beyond) Overview Chapter 1 built the numerical foundations (norms, conditioning, direct and stationary iterative solvers) for Ax = b. In this chapter we will see where the matrix A actually comes from, and why it is almost always sparse: the discretization of a PDE by finite differences, the Laplacian of a graph, and the normal equations of a machine learning model all produce matrices whose nonzero pattern is dictated by local structure - a stencil, an edge list, or a feature’s support. This chapter covers: (1) from Taylor expansions to finite-difference formulas and the 1D tridiagonal Poisson system; (2) the 2D Poisson matrix, graph Laplacians, and sparsity in machine learning; (3) a toolbox for working with sparse matrices in practice - storage formats, bandwidth/fill-in, reordering, and the spectral properties that govern solver convergence. Learning objectives By the end of this chapter, you should be able to: • Derive finite-difference formulas from Taylor expansions and state their consistency order (truncation error). • Assemble the tridiagonal (1D) and block-tridiagonal/Kronecker (2D) linear systems that arise from discretizing the Poisson equation. • Define the graph Laplacian L = D − A and relate its sparsity pattern and spectrum to those of PDE and machine-learning matrices. • Explain why sparsity appears in ML (normal equations, Hessians) and data assimilation (localized correlation matrices). • Choose an appropriate sparse storage format (CSR/CSC, banded), and explain how bandwidth, profile, and reordering affect fill-in and solver cost. • Relate the spectral properties (eigenvalues, condition number) of a sparse matrix to the convergence of iterative solvers and the need for preconditioning.

2.1

Finite difference approximation

In Chapter 1 we discussed the importance of numerical solvers for differential equations. Some differential equations can be solved analytically (e.g. separable partial differential equations, or using the method of characteristics), giving solution in closed form. However, many PDEs 15

16

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

u(x) x−h x x+h

x

local Taylor expansion

Figure 2.1: Zoom-in of u(x) around a point with neighbors x ± h

of interest for applications cannot be solved analytically. In this case, we solve the problem approximately by discretising the continuous problem, which can then be solved on a computer. Three of the most common approaches include the finite element method, finite volume method, and finite difference method. In this course, we focus on finite difference approximations and consider how the properties of the resulting matrices play a role in the convergence of iterative methods. Long before computers, finite differences were already a powerful idea. Euler used them to tabulate functions, and a century later Richardson boldly extended the same spirit to the whole atmosphere. In his 1922 Weather Prediction by Numerical Process, he replaced derivatives by local difference operators to turn the equations of fluid dynamics into linear systems, attempting a six-hour forecast by hand. The calculation took him weeks and produced unstable and incorrect results, but his dream came back to life again with the advent of parallel computers [41]. Let u : R → R be smooth and h > 0. A Taylor expansion around x gives the following result. 2

3

4

u(x ± h) = u(x) ± hu′ (x) + h2 u′′ (x) ± h6 u(3) (x) + h24 u(4) (ξ± ), for some ξ± between x and x ± h. From this we obtain the following. Definition 2.1.1 (Finite difference operators). Forward: Backward: Central: Second derivative:

u(x + h) − u(x) = u′ (x) + h2 u′′ (x) + O(h2 ), h u(x) − u(x − h) = u′ (x) − h2 u′′ (x) + O(h2 ), h u(x + h) − u(x − h) = u′ (x) + O(h2 ), 2h u(x + h) − 2u(x) + u(x − h) = u′′ (x) + O(h2 ). h2

Definition 2.1.2 (Consistency and truncation error). A finite difference operator Lh is consistent with a differential operator L if (Lh u)(x) = (Lu)(x) + τh (x) with ∥τh ∥ → 0 as h → 0. The term τh is the local truncation error ; its order (O(h), O(h2 ),. . . ) is the accuracy. Proposition 2.1.3 (Consistency and order of the three-point stencil). Let u ∈ C 4 in a neighborhood of x ∈ R. Then h Dxx u(x) :=

u(x + h) − 2u(x) + u(x − h) h2 (4) ′′ = u (x) + u (ξ) h2 12

for some ξ ∈ (x − h, x + h).

h u(x) is consistent with u′′ (x) and has truncation In particular, the finite difference operator Dxx error O(h2 ) as h → 0.

17

2.1. FINITE DIFFERENCE APPROXIMATION [−1, 2, −1]/h2

α

u1

u2

···

un−1

un

+2

−1

β

h

i−1

i

−1 i+1

x

x

Figure 2.2: Left: Function values on discretised mesh. Right: Stencil for central differences in 1D

Proof. By Taylor’s theorem with Lagrange remainder, for some ξ+ ∈ (x, x+h) and ξ− ∈ (x−h, x), u(x ± h) = u(x) ± hu′ (x) +

h2 ′′ h3 h4 u (x) ± u(3) (x) + u(4) (ξ± ). 2 6 24

Adding the + and − expansions and subtracting 2u(x) yields u(x + h) − 2u(x) + u(x − h) = h2 u′′ (x) +

h4 (4) u (ξ+ ) + u(4) (ξ− ) . 24

Dividing by h2 gives
u(x + h) − 2u(x) + u(x − h) h2 (4) ′′ = u (x) + u (ξ+ ) + u(4) (ξ− ) . 2 h 24
By the intermediate value theorem, there exists ξ ∈ (x−h, x+h) such that 12 u(4) (ξ+ )+u(4) (ξ− ) = u(4) (ξ), hence u(x + h) − 2u(x) + u(x − h) h2 (4) = u′′ (x) + u (ξ). 2 h 12 Therefore h2 h Dxx u(x) − u′′ (x) ≤ max |u(4) (y)| = O(h2 ), 12 |y−x|≤h which proves consistency with second-order accuracy. Remark 2.1.4. Central differences are second order accurate for u′ and u′′ , while forward/backward are first order for u′ . Higher order formulas follow by eliminating higher derivatives from longer Taylor stencils. Consider now the boundary value problem (BVP) − u′′ (x) = f (x) for x ∈ (0, 1),

u(0) = α, u(1) = β.

Let xi = ih with h = 1/(n + 1) and unknowns ui ≈ u(xi ) for i = 1, . . . , n. Approximating u′′ (xi ) with the central difference gives −

ui−1 − 2ui + ui+1 = f (xi ), h2

i = 1, . . . , n,

where u0 = α and un+1 = β are known from the Dirichlet boundary conditions. Proposition 2.1.5 (Tridiagonal linear system). Collecting the n interior equations yields A u = g,

A=

1 tridiag(−1, 2, −1) ∈ Rn×n , h2

with g = (f (xi ))ni=1 modified at the endpoints: g1 ← g1 + α/h2 , gn ← gn + β/h2 .

18

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

Explicit form of the matrix. For n = 5, the system has the form      g1 u1 2 −1 0 0 0  u2  g2  −1 2 −1 0 0     1   0 −1 2 −1 0  u3  = g3  .      2 h  0 0 −1 2 −1 u4  g4  g5 u5 0 0 0 −1 2 The sparsity pattern corresponds exactly to the stencil [−1, 2, −1]/h2 . Takeaways • Finite difference discretisation converts differential operators into sparse matrices whose nonzero pattern is determined by the stencil. • The stencil is therefore the object to look at: it fixes the bandwidth, and with it the cost of every direct method applied to the resulting system.

2.2

Sparse Matrices: PDEs, Graphs, and Machine Learning

In 1928, Courant, Friedrichs, and Lewy introduced the finite difference method for PDEs in their seminal work. At the same time, scientists like Kirchhoff (1847) were studying electrical networks using adjacency matrices and Laplacians. Nearly a century later, sparse matrices from both PDEs and graphs lie at the core of modern scientific computing and machine learning.

2.2.1

The 2D Poisson equation: structure, stencils and blocks

Consider the Poisson problem on the unit square −∆u(x, y) = f (x, y),

(x, y) ∈ (0, 1)2 ,

u|∂Ω = 0.

1 Let xi = ih, yj = jh, with h = n+1 and interior indices i, j = 1, . . . , n. Denote the grid unknowns by ui,j ≈ u(xi , yj ). We discretize the two second derivatives separately:

uxx (xi , yj ) ≈

ui−1,j − 2ui,j + ui+1,j , h2

Thus −∆u(xi , yj ) ≈

uyy (xi , yj ) ≈

ui,j−1 − 2ui,j + ui,j+1 . h2

4ui,j − ui−1,j − ui+1,j − ui,j−1 − ui,j+1 . h2

Definition 2.2.1 (Kronecker product). For matrices A ∈ Rm×n and B ∈ Rp×q , their Kronecker product is the block matrix   a11 B a12 B · · · a1n B  a21 B a22 B · · · a2n B    mp×nq A⊗B =  . . .. ..  ∈ R . . .  . . . .  am1 B am2 B · · ·

amn B

From grid to linear system. Collect the n2 interior equations into A2D u = g. Using row-wise (lexicographic) ordering p = (j − 1)n + i, we obtain the block tridiagonal matrix  1 A2D = 2 I ⊗ T + T ⊗ In , T = tridiag(−1, 2, −1) ∈ Rn×n . h

2.2. SPARSE MATRICES: PDES, GRAPHS, AND MACHINE LEARNING

19

(i, j +1)

−1

(i, j)

(i−1, j) (i, j) (i+1, j)

−1

+4

−1

(i, j −1)

−1 Figure 2.3: Left: Indexing of grid. Right: stencil to update u(i, j) Written blockwise, the same matrix has the structure 

B

−I

 −I 1   A2D = 2  0 h   ..  .

B −I .. .

0

···

··· .. . −I .. . B .. .. . . 0 −I 0

 0 ..  .    . 0   −I  B

Each block B is itself tridiagonal with stencil [−1, 4, −1], and each −I couples neighboring rows of the grid. Remark 2.2.2. A2D is sparse (each row has at most 5 nonzeros), symmetric, and positive definite (SPD) under Dirichlet data. Boundary conditions. For homogeneous Dirichlet (u = 0 on ∂Ω) the boundary values are known and do not appear as unknowns. For inhomogeneous Dirichlet (u = g), add boundary contributions to the right hand side near the boundary; sparsity is preserved. One sided differences handle Neumann or Robin conditions, modifying only the first/last block rows.

2.2.2

Graph Laplacians

Graphs give rise to another fundamental family of sparse matrices. Let G = (V, E) be an undirected graph with |V | = n vertices and edges E. The adjacency matrix A is defined by ( 1 if (i, j) ∈ E, Aij = 0 otherwise. The degree matrix D is diagonal with entries Dii = deg(i) =

j Aij .

P

Graph Laplacian The graph Laplacian of G is L = D − A. It is symmetric, positive semidefinite, and sparse: nnz(L) = 2|E| + n.

20

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

Example. For the path graph on 5 vertices: 1

2

3

4

5

the adjacency matrix is  0 1  A= 0 0 0 so

1 0 1 0 0

0 1 0 1 0

0 0 1 0 1

 0 0  0 , 1 0

D = diag(1, 2, 2, 2, 1),



 1 −1 0 0 0 −1 2 −1 0 0   . 0 −1 2 −1 0 L=D−A=   0 0 −1 2 −1 0 0 0 −1 1

Remark 2.2.3. This L is a rank-2 perturbation of the tridiagonal matrix arising from the 1D finite difference discretisation of −u′′ (x). Topology changes the pattern. The path graph is only one extreme; the cycle Cn (periodic path) and the star K1,n−1 (one hub, n − 1 leaves) illustrate how differently L can look for the same number of vertices. 2 3

2

4

6 1

5

6

Cycle graph C6

3

1 5

4

Star graph K1,5

For the cycle Cn , L is a circulant tridiagonal matrix with wrap-around corner entries; its eigenvalues are λk = 2 − 2 cos(2πk/n) (see Section 2.3). For the star K1,n−1 , the hub has degree n − 1 and every leaf has degree 1, giving the sparse but highly non-tridiagonal spectrum σ(L) = {0, 1, . . . , 1, n}: a single dominant eigenvalue reflects the hub’s central role, in sharp | {z } n−2

contrast to the path and cycle. On a path graph, the Laplacian matrix “almost” coincides with the finite difference stencil of the PDE Laplacian.a a More rigorously, this would correspond to the discrete version of the Laplacian with homogeneous Neumann boundary conditions

Graph Laplacians are put to work in several directions. • Spectral clustering: eigenvectors of L reveal community structure in networks.

2.2. SPARSE MATRICES: PDES, GRAPHS, AND MACHINE LEARNING

21

• Graph neural networks: message passing relies on L or its normalized variants. • Electrical networks: resistor networks are governed by Kirchhoff’s laws, which are encoded in L. • Data assimilation: L is used to construct diffusion operators, which are commonly used as spatial covariance matrices in meteorological and oceanographic data assimilation. Takeaways Graph Laplacians are sparse and symmetric matrices that extend the notion of the Laplace operator from continuous domains to arbitrary discrete networks. Their role as “discrete PDEs on graphs” explains their central place in modern data science and network analysis.

2.2.3

Sparse Matrices in Machine Learning

Machine learning problems generate some of the largest sparse matrices encountered today. The source of sparsity is locality: each data point, feature, or parameter interacts only with a small subset of the whole system. In least squares regression, for instance, min ∥Ax − b∥22 , x

the optimality condition yields the normal equations A⊤ Ax = A⊤ b, where A ∈ Rm×n is the feature matrix and b ∈ Rm the labels. Sparsity of A and A⊤ A If each row of A contains only k ≪ n nonzeros (each sample depends on k features), then nnz(A) = O(mk) and A is extremely sparse. The matrix A⊤ A is SPSD and encodes feature correlations. It is usually denser than A, but it still preserves significant sparsity compared to a dense n × n matrix. Remark 2.2.4. In text classification, A is the document-term matrix : rows correspond to documents, columns to words, and entries record word frequencies. Each document contains only a few words compared to the dictionary size ⇒ A is sparse. Then A⊤ A counts word co-occurrences across documents, producing a sparse graph of word relations.

A A⊤ A Sparse design matrix A (left) vs. denser normal matrix A⊤ A (right). Second-order optimisation methods require the Hessian Hij =

∂2 F (x). ∂xi ∂xj

When the loss F couples only nearby features, H inherits sparsity. This occurs frequently in PDE-constrained optimization, and deep learning layers with local connectivity.

22

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

Example: Logistic regression For logistic regression with data matrix A (columns ak ) and labels bk ∈ {−1, 1}, F (x) =

m X

log 1 + exp(−bk a⊤ k x) ,

k=1

the Hessian takes the form H = A⊤ W A where W is a diagonal matrix of weights depending exp b a⊤ x on x (wk = (1+exp bk ak⊤ x)2 ). If A is sparse, so is H, and computations can be done efficiently. k k

Takeaways In machine learning, sparsity comes from locality in data or models. Normal equations, Hessians, and Jacobians inherit this structure. Harnessing it is essential: without sparse techniques, modern large-scale ML training would be infeasible.

2.2.4

Sparsity in data assimilation

In the final part of the course we will see that covariance matrices are important for weighting the contribution of different components to a data assimilation system. The sparsity structure of these covariance matrices has important computational and physical consequences. In ensemble data assimilation methods, underlying ‘physical’ covariances are obtained using sample estimates. This can lead to spurious long-range correlations which destroy sparsity. In a later part of the course we will consider localisation methods, which are used to prescribe a given local sparsity structure on estimated spatial covariance matrices.

Sometimes we wish to estimate the value of the covariance matrix from measurement data. One method is to use the measurements to fit appropriate parameter values for a standard correlation function (e.g. Gaussian, auto-regressive functions) see e.g. [25]. These distributions all have long tails and therefore yield dense matrices. In practice, correlation distributions are often truncated when they fall below a threshold value (typically chosen as 0.2) to ensure sparsity/fixed bandwidth of the resulting matrices.

2.3

Sparse Matrix Toolbox

We have now seen where sparse matrices come from: discretised PDEs, graphs, and machine learning. Knowing where they arise is not enough to compute with them. We also need to know how to store them, how to operate on them efficiently, and how their spectral properties affect the algorithms of the following chapters. That is the subject of this section.

2.3.1

Sparsity Patterns and storage formats

Definition The sparsity pattern of a matrix A is the set of positions (i, j) where aij ̸= 0. Example. For the tridiagonal finite difference Poisson matrix: nnz(A) ≈ 3n,

density(A) = nnz(A) ≈ 1/n. n2

23

2.3. SPARSE MATRIX TOOLBOX

1D Poisson (tridiagonal). For n interior points, A = h12 tridiag(−1, 2, −1) has nonzeros only on the main diagonal and the first sub-/super-diagonals.

Spy pattern for 1D Poisson, n = 9. Bandwidth = 1, nnz/row ≈ 3 (except at ends). 2D Poisson with lexicographic ordering. Let the grid be n × n (N = n2 unknowns) and use row-wise (lexicographic) ordering. Then A has nonzeros at indices j = i (main), j = i ± 1 within the same row block, and j = i ± n (coupling to the row above/below). That is exactly the 5-point stencil.

Spy pattern for 2D Poisson, n = 4. Bandwidth = 4, nnz/row ≈ 5 (except at end blocks). Graph Laplacians (path vs. star). Two small graphs with n = 8 nodes to highlight how topology changes the pattern. Path graph L

Star graph L

Left: path L ⇒ tridiagonal. Right: star L ⇒ dense first row/column + diagonal.

24

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

Matrix structure and storage. Sparse matrices are not stored entry by entry; only the nonzero pattern is kept. Common formats include: • CSC/CSR (Compressed Sparse Column/Row): general-purpose, store indices and values compactly. • Diagonal/banded storage: efficient when nonzeros cluster near the diagonal. Example 2.3.1 (CSR and CSC formats). Consider the 4 × 4 sparse matrix   1 0 0 2 0 3 0 4  A= 5 6 7 0 . 0 0 8 9 CSR (Compressed Sparse Row). val = [1, 2, 3, 4, 5, 6, 7, 8, 9], col_ind = [1, 4, 2, 4, 1, 2, 3, 3, 4], row_ptr = [1, 3, 5, 8, 10].

CSC (Compressed Sparse Column). val = [1, 5, 3, 6, 7, 8, 2, 4, 9], row_ind = [1, 3, 2, 3, 3, 4, 1, 2, 4], col_ptr = [1, 3, 5, 7, 10].

Here, indices are 1-based for readability (note that Python/NumPy uses 0-based indexing). - In CSR, row_ptr[i] points to the start of row i in val. - In CSC, col_ptr[j] points to the start of column j in val. Bandwidth vs. profile. Two measures quantify “how far” nonzeros lie from the diagonal: • Bandwidth: b = max{|i − j| : aij ̸= 0}. The worst-case distance from the diagonal. Narrow bands allow cheap banded LU; wide bands inflate fill and storage. • Profile (envelope): sum of distances of each nonzero to the diagonal. This reflects actual memory usage more sensitively than bandwidth alone. Remark 2.3.2. In the 2D Poisson problem with N = n2 , each row has only 5 nonzeros, yet the bandwidth grows like n. Naive banded LU then costs O(N b2 ) = O(N 2 ), whereas good reorderings (e.g. nested dissection) reduce the cost to O(N 3/2 ). Fill-in. During LU factorization, zeros in A can become nonzeros in L or U . • For banded matrices, fill grows with bandwidth. • In general, the amount of fill depends strongly on the sparsity pattern.

A LU Illustration: fill-in introduces new nonzeros (red).

25

2.3. SPARSE MATRIX TOOLBOX

0

0

2

2

4

4

6

6

8

8

10

10 0

5

10

0

nz = 20

0

0

2

2

4

4

6

6

8

8

10

10 0

5

nz = 19

5

10

nz = 91

10

0

5

10

nz = 18

Figure 2.4: Visualisation of fill-in with LU vs inverse. Top left: sparsity pattern original matrix, A. Top right: sparsity pattern of A−1 . Lower left: sparsity pattern of L. Lower right: sparsity pattern of U . Remark 2.3.3. By permuting rows/columns before factorization, one can reduce bandwidth and hence fill-in. Effective reordering can save orders of magnitude in memory and runtime—critical for large-scale finite element systems. The invention of compressed storage and reordering heuristics in the 1970s was pivotal: without them, sparse PDE systems of size 105 could not be solved on mainframes. The inverse of a sparse matrix is not guaranteed to be sparse, and can suffer from more fill-in than an LU or Cholesky decomposition. One example is the given in Figure 2.4, where A is a lower triangular matrix with ones on the diagonal and first subdiagonal, plus a single entry above the diagonal. The LU factorization of A is sparse, with some fill-in seen in U (lower right panel). This aligns with the element above the diagonal. On the other hand, A−1 is a dense matrix, with almost all entries being non-zero. In general we aim to avoid using inverse matrices, but in the final section of the course we will see that inverse covariance matrices play a very important role in data assimilation and cannot be completely avoided.

2.3.2

Spectral Properties

For some particular sparsity patterns, we can compute the eigenvalues explicitly: Tridiagonal matrices If a tridiagonal matrix is Toeplitz i.e. with stencil [c, a, b] then the eigenvalues can be computed explicitly as: √ jπ λj = a − 2 bc cos( ) j = 1, . . . , n. n+1

26

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

Definition A circulant matrix is a particular form of Toeplitz matrix which has rows composed of all the same elements, with each row being shifted one element to the right relative to the preceding row, i.e.   c0 c1 . . . cn−2 cn−1 cn−1 c0 c1 cn−2     .. ..  ..  . cn−1 c0 .  C= . .   . . . .  c2 . . c1  c1 c2 . . . cn−1 c0 Circulant matrices are diagonalized by a discrete Fourier transform, i.e. the eigenvalues are given by λj = c0 + cn−1 ω −j + cn−2 ω −2j + · · · + c1 ω −(n−1)j

j = 0, . . . , n − 1,

where ω = exp( 2πi n ) is an n − th root of unity. The eigenvectors of C are given by 1 vj = √ (1, ω j , ω 2j , . . . , ω (n−1)j )T , n

j = 0, 1, 2, . . . , n − 1.

Eigenvalues of model problems. For the 1D Poisson matrix,   4 jπ 2 λj = 2 sin , j = 1, . . . , n. h 2(n + 1) Hence κ(A) =

λmax = O(n2 ). λmin

Similarly for the 2D diffusion matrix, A = I − hν2 ∈ Rn×n , with Dirichlet boundary conditions on the unit square, the n eigenvalues of A are given by      4ν i j 2 2 µi,j = 1 + 2 sin π + sin π , i, j = 1, 2, . . . , n. h 2(n + 1) 2(n + 1) Implication These condition numbers grow as h−2 : iterative solvers (CG, GMRES, etc.) slow down as the mesh refines. Spectral clustering of eigenvalues governs convergence rates: if most eigenvalues lie in a tight interval, Krylov methods converge quickly; wide spreads or outliers cause stagnation. Preconditioners are therefore essential to compress the spectrum. Spectral viewpoint on sparsity. Even though A is sparse, its eigenvalues may cluster near zero, leading to ill-conditioning. Sparsity alone does not guarantee fast convergence: what matters is the distribution of eigenvalues. For example, reordering does not change eigenvalues (it is a similarity transformation), but preconditioning actively reshapes the spectrum to improve solver efficiency. Graph Laplacians. Eigenvalues encode combinatorial structure: • Smallest eigenvalue = 0 (multiplicity = number of connected components). • The second smallest eigenvalue (Fiedler value) measures algebraic connectivity and governs diffusion/mixing rates.

2.3. SPARSE MATRIX TOOLBOX

27

Machine learning. Normal equations A⊤ A have eigenvalues σi2 , the squares of the singular values of A. Ill-conditioning corresponds to nearly collinear features, where small singular values slow down gradient descent and amplify noise. Regularization (ridge regression, Tikhonov) shifts the spectrum, lifting small eigenvalues away from zero. Takeaways • Sparsity patterns (bandwidth, profile) determine storage and fill-in. • Compressed storage formats (CSR/CSC) make large problems feasible. • Fill-in can destroy sparsity; reordering strategies mitigate this. • Spectral properties are decisive: eigenvalues dictate conditioning, convergence, and the design of effective preconditioners. Efficient sparse linear algebra requires both algebraic insight (spectrum, conditioning) and practical storage/reordering techniques.

28

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

2.4

Exercises

Problem 1 (FD accuracy at a point). For u(x) = sin(πx), approximate u′′ (0.5) with the 3-point stencil u(0.5 + h) − 2u(0.5) + u(0.5 − h) h Dxx u(0.5) = h2 using h = 0.25 and h = 0.125. Compare with the exact u′′ (0.5) and the error ratio. Problem 2 (1D Poisson: assemble & solve by hand). Discretize −u′′ = 1 on (0, 1) with u(0) = u(1) = 0 and n = 4 interior points (h = 1/5). 1. Write the linear system Au = g. 2. Solve it explicitly for u ∈ R4 . Problem 3 (Discretisation of a BVP and application of a basic iterative method, exam 2024/2025). Consider the one-dimensional boundary value problem ( − u′′ (x) + c u(x) = 0, x ∈ (0, 1), u(0) = 1, u(1) = 0, 1 where c > 0 is a constant. Use a uniform mesh with grid points xj = jh, h = n+1 , j = 0, . . . , n + 1, and denote uj ≈ u(xj ).

(a) Give the finite difference scheme for the discretisation of the BVP. (b) Let u = (u1 , . . . , un )⊤ be the unknown vector solving Au = b. Determine A and b. (c) We solve the linear system by the Jacobi method based on the splitting A = D − E − F , where D is the diagonal part and E (resp. F ) the strict lower (resp. upper) triangular part. Give the exact Jacobi iteration. (d) If ek = u − uk denotes the error at iteration k, then ek = Gk e0 with iteration matrix G. Give G explicitly. Problem 4 (Finite Differences for Boundary Value Problems, resit exam 2024/2025). Consider the boundary value problem: −u′′ (x) + cu(x) = f (x),

x ∈ (0, 1),

u(0) = α,

u(1) = β,

where c is a positive constant, and f (x) is a given source term. (a) Discretize the equation using the finite difference method with a uniform grid xj = jh, j = 0, 1, . . . , n + 1, where h = 1/(n + 1). Derive the corresponding difference scheme. (b) Write the scheme in matrix form. Give the expression of the matrix and right hand side. (c) Compute the spectral condition number of the matrix A. (A) Hint: For a symmetric positive definite matrix we have: κ(A) = λλmax , where λmax (A) min (A) and λmin (A) are the largest and smallest eigenvalues of A, respectively. For a tridiagonal matrix of the form A = tridiag(a, b, a), the eigenvalues are given by:   kπ λk = b + 2a cos , k = 1, 2, . . . , n. n+1 This formula can be used to compute λmax and λmin . (d) Analyze how the spectral condition number κ(A) behaves when c decreases and n increases and interpret the result.

29

2.4. EXERCISES

Problem 5 (Bandwidth and naive LU cost). For a 2D n × n grid with natural ordering, take n = 20 (N = n2 ). Estimate the semi-bandwidth b and the work of banded LU via O(N b2 ). Problem 6 (Graph spectra, small and exact). For the path graph on 4 vertices, compute the eigenvalues of its Laplacian L. Then do the same for the complete graph K4 . Compare the spectral gaps. Problem 7 (2D Laplacian on an n × n grid: spectrum and conditioning). Let Tn = tridiag(−1, 2, −1) (Dirichlet 1D Laplacian without 1/h2 scaling) and 2

2

A2D = Tn ⊗ In + In ⊗ Tn ∈ Rn ×n . (i) Show that the eigenvalues of A2D are  pπ   qπ  λp,q = λp (Tn ) + λq (Tn ) = 2 − 2 cos + 2 − 2 cos , n+1 n+1

p, q = 1, . . . , n.

b2D = 12 A2D . Prove that (ii) With the standard spacing h = 1/(n + 1), the stiffness matrix is A h b2D ) ∼ 8 , λmax (A h2

b2D ) → 2π 2 λmin (A

(n → ∞),

b2D ) = O(n2 ). and deduce the asymptotic condition number κ2 (A (iii) (Optional, dD generalization.) For d spatial dimensions with nd unknowns, show that bdD ) = O(n2 ) still holds. κ2 (A Problem 8 (Graph Laplacians: spectra and expansion, general n). Let G = (V, E) be a simple undirected graph on n vertices with (combinatorial) Laplacian L = D − A. (i) Show that 0 is an eigenvalue of L with multiplicity equal to the number of connected components. (ii) For the cycle path graph Cn , provide the eigenvalues and eigenvectors of Cn in terms of ω = exp( 2πi n . (You may use the properties of circulant matrices.) Give an interval in which these eigenvalues lie. (iii) For the star graph Sn , show that σ(L) = {0, 1 (mult n−2), n}. Problem 9 (1D Poisson eigenvalues and κ(A), small n). For n = 5 interior points, use  jπ  4 2 λj = 2 sin , h 2(n + 1)

h=

1 , n+1

to list λ1 , . . . , λ5 numerically and compute κ(A). Repeat this for      4ν i j 2 2 µi,j = 1 + 2 sin π + sin π , h 2(n + 1) 2(n + 1)

i, j = 1, 2, . . . , n.

with ν = 0.0025. What differences do you notice between the two sets of eigenvalues? What will happen to the largest/smallest eigenvalues as n increases?

30

CHAPTER 2. DISCRETIZATION OF 1D AND 2D PDES (AND BEYOND)

Problem 10 (Data assimilation spectrum example). We can construct a second order autoregressive (SOAR) correlation function on the unit circle for N equally distributed gridpoints, via ! ! θi,j θ −|2a sin( |2a sin( i,j )| )| 2 2 exp , 1 ≤ i, j ≤ N C(i, j) = 1 + L L where L > 0 is the correlation lengthscale, θi,j denotes the angle between grid points i and j, and a is the radius of the domain [15]. This yields a dense matrix, which can be ill-conditioned. What happens to the spectrum if we replace all values less than some threshold (e.g. 0.2) by zero? Investigate how this changes as the lengthscale, threshold value and problem dimension increase. Problem 11 (Concrete ML sparsity example). Four documents over a 6-word vocabulary with bags: d1 = {1, 2}, d2 = {2, 3}, d3 = {3, 4}, d4 = {1, 4}. Build the sparse document-term matrix A ∈ R4×6 (entries 0/1). Compute A⊤ A explicitly and report its sparsity (which off-diagonal entries are nonzero?). Interpret the pattern.

Chapter 3

Applications of Basic Iterative Methods Overview Chapter 1 introduced Jacobi, Gauss–Seidel, SOR, and Richardson as instances of a single fixedpoint iteration x(k+1) = Gx(k) + c, convergent whenever ρ(G) < 1. This chapter shows how this works on three model problems where the spectrum of G can be computed explicitly, which turns the abstract convergence criterion into concrete, comparable convergence rates. This chapter covers: (1) the 1D Poisson matrix, where explicit eigenpairs give exact convergence factors for Jacobi/GS/SOR; (2) graph Laplacians, where the same iterations become diffusion/random-walk processes and connect to spectral clustering; (3) machine learning, where Richardson iteration is exactly gradient descent and the spectrum of X ⊤ X governs both convergence speed and spectral bias. Learning objectives By the end of this chapter, you should be able to: • Derive the explicit eigenvalues of the 1D Poisson matrix and use them to compute the exact spectral radius (and hence convergence rate) of Jacobi, Gauss–Seidel, and SOR. • Explain why Gauss–Seidel converges roughly twice as fast as Jacobi, and why SOR with an optimal ω improves the rate further. • Interpret Jacobi iteration on a graph Laplacian as a random walk / diffusion process, and relate its slow modes to the Fiedler vector and spectral clustering. • Recognize Richardson iteration as gradient descent on a quadratic loss, and relate the spectrum of A (or X ⊤ X) to the choice of step size and the resulting convergence rate. • Explain the phenomenon of spectral bias in machine learning as the same spectralfiltering mechanism that governs smoothing in PDE solvers, with the roles of low- and high-frequency modes swapped.

3.1

Basic Iterative Methods for PDE discretisations

We now go beyond the general fixed point framework of Chapter 1 and analyze in detail the convergence of classical iterative methods for the model problem of the one-dimensional Poisson equation with Dirichlet boundary conditions introduced in Chapter 2. This case is simple enough to allow explicit eigenvalue analysis, yet rich enough to illustrate the key phenomena that will reappear in higher–dimensional PDEs, graphs, and machine learning applications. Recall the discretization of −u′′ = f on (0, 1) with homogeneous Dirichlet boundary conditions 31

32

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

1 using n interior points, mesh width h = n+1 . The stiffness matrix is



2 −1 −1 2 −1 1  1  .. .. A = 2 tridiag(−1, 2, −1) = 2  . . h h   −1

    ..  ∈ Rn×n . .  2 −1 −1 2

(3.1)

This matrix is symmetric, positive definite, and tridiagonal. Remark 3.1.1. The eigenvalues of A are explicitly known:
λk (A) = h12 2 − 2 cos θk ,

kπ , k = 1, . . . , n. n+1

θk =

(3.2)

or alternatively λk (A) =

    kπ 4 4 2 2 kπh = , sin sin h2 2(n + 1) h2 2

k = 1, . . . , n.

Hence λmin ∼ π 2 , λmax ∼ 4/h2 , and κ(A) = O(n2 ). low k j high k (k)

Low-frequency (k small) vs. high-frequency (k large) eigenvectors vj = sin(jkπ/(n + 1)) of A: high-frequency modes oscillate rapidly and are damped fast by Jacobi/GS/SOR, while low-frequency modes decay slowly — the smoothing property used later in multigrid [7].

3.1.1

Convergence of the Jacobi iteration

We split A = D − L − U with D = h22 I. The Jacobi iteration matrix is   2 GJ = I − D−1 A = I − h2 A = I − 21 tridiag(−1, 2, −1) = tridiag 12 , 0, 12 . 2

Since GJ = I − h2 A and given (3.2) the Jacobi eigenvalues are  
2 2 kπ µk = 1 − h2 λk (A) = 1 − h2 · h12 2 − 2 cos θk = cos θk = cos n+1 . Hence the spectral radius is   π ρ(GJ ) = max |µk | = cos n+1 =1− k

π2 + O n−4 , 2 2(n + 1)

so Jacobi converges very slowly as n grows.

3.1.2

Convergence of the Gauss–Seidel iteration

For Gauss–Seidel we split A = (D − L) − U , with D = h22 I,

L = h12 tridiag(1, 0, 0),

U = h12 tridiag(0, 0, 1).

Then the iteration matrix is (see also the first exercise on the list for the exact formula) GGS = (D − L)−1 U.

3.1. BASIC ITERATIVE METHODS FOR PDE DISCRETISATIONS

33

Theorem 3.1.2 (Relation between eigenvalues of Gauss–Seidel and Jacobi). Let A = D−L−U with D diagonal, −L strictly lower, −U strictly upper, and all diagonal entries of D nonzero. Define GJ = D−1 (L + U ), GGS = (D − L)−1 U. If µ is an eigenvalue of GJ , then λ = µ2 is an eigenvalue of GGS . Proof. By definition,

GGS = (I − D−1 L)−1 D−1 U.

Hence λ ∈ σ(GGS ) if and only if
det (I − D−1 L)−1 D−1 U − λI = 0, which is equivalent to det D Thus

√

−1

U − λ(I − D

λ is an eigenvalue of

−1

   √ √
1 −1 L) = 0 ⇐⇒ det = 0. λI − D λL + √ U λ

(1)

 √  M (λ) := D−1 √1λ U + λ L .

Now observe that M (λ) is similar to D−1 (L + U ). 1 Therefore σ(M (λ)) = σ(D−1 (L + U )) = σ(GJ ). √ √ Hence λ ∈ σ(GJ ). Same holds for − λ. Conclusion. The Gauss–Seidel eigenvalues are   kπ µk = cos2 n+1 , and the spectral radius ρ(GGS ) = ρ(G2J ) hence Gauss-Seidel converges twice as fast as Jacobi. Remark 3.1.3 (Spectrum and non-diagonalisability of Gauss–Seidel). For the one-dimensional Poisson problem with n interior points and natural ordering, the Jacobi iteration matrix GJ is symmetric and diagonalisable. By contrast, the Gauss–Seidel iteration matrix is non-normal and not diagonalisable. Its spectrum is n   o kπ σ(GGS ) = cos2 n+1 : k = 1, . . . , ⌊n/2⌋ ∪ {0 (with algebraic multiplicity ⌈n/2⌉)}. Thus GGS has only ⌊n/2⌋ distinct nonzero eigenvalues, each simple, together with the eigenvalue 0. The eigenvalue 0 has algebraic multiplicity ⌈n/2⌉ but geometric multiplicity one, so Jordan blocks of size > 1 occur. This shows that GGS is not diagonalisable.

3.1.3

Convergence of the SOR iteration

For
SOR with relaxation parameter ω ∈ (0, 2) we use the splitting A = ω1 (D − ωL) − (1 − ω)D + ωU , so that the iteration matrix is
Gω = (D − ωL)−1 (1 − ω)D + ωU . 1

Here we use a bit a shortcut and we skip some steps. In fact this is true for matrices which have property A (i.e. can be reduced to a block structure via a permutation of rows and columns. In this proof we assume that this permutation has already been done).

34

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

Young’s theorem [42, 37]2 gives for certain matrices a relation between the SOR convergence factor and that of Jacobi and the optimal parameter ω: p 1 − 1 − ρ(GJ )2 2 p p ωopt = , ρ(Gωopt ) = . 1 + 1 − ρ(GJ )2 1 + 1 − ρ(GJ )2
π For our 1D Poisson problem, ρ(GJ ) = cos n+1 , which yields
π 1 − sin n+1 2
, ρ(Gωopt ) =
. ωopt = π π 1 + sin n+1 1 + sin n+1 Proposition 3.1.4 (Asymptotic Convergence and Limitations). Based on the previous results we can derive the asymptotic convergence behaviour of the three classical iterative methods. For large n, the spectral radii expand as 2

−4 π ), ρ(GJ ) = 1 − 2(n+1) 2 + O(n 2

−4 π ρ(GGS ) = 1 − (n+1) ), 2 + O(n 2π ρ(Gωopt ) = 1 − n+1 + O(n−2 ),

Example 3.1.5 (Numerical comparison, n = 49). Take n = 49 interior points (so h = 1/50). The three spectral radii evaluate to   π ρ(GJ ) = cos 50 ≈ 0.99803, ρ(GGS ) = ρ(GJ )2 ≈ 0.99606, ρ(Gωopt ) ≈ 0.88180. To reduce the initial error by a factor 10−6 , the number of iterations k needed satisfies ρk ≤ 10−6 , i.e. k ≥ −6/ log10 ρ. This gives kJ ≈ 7000,

kGS ≈ 3500,

kωopt ≈ 110.

So for this modest grid, optimally relaxed SOR needs roughly 64× fewer iterations than Jacobi to reach the same accuracy which is a direct, numerical illustration of the asymptotic rates above. Even so, 110 iterations of a method that costs O(n) per sweep is still far more expensive than the direct O(n) Thomas algorithm available for this tridiagonal system; the real payoff of iterative methods appears in 2D/3D, where direct factorization becomes infeasible but the same spectral machinery still applies. Takeaways • Jacobi converges very slowly: ρ(GJ ) → 1 as n → ∞, with rate O(1/n2 ). • Gauss–Seidel roughly doubles the effective rate, but is still O(1/n2 ). • SOR with optimal relaxation improves to O(1/n) per iteration, but this is still poor compared to modern solvers (e.g. multigrid or Krylov methods). Thus none of these methods are suitable as standalone solvers for large PDE systems.

3.2

Basic Iterative Methods for Graph Laplacians

Graph Laplacians defined in Chapter 2 provide another fundamental class of sparse matrices. We recall that for any graph, the smallest eigenvalue of L is given by 0, with corresponding 2

we accept this result without proof

35

3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS

eigenvector e = [1, 1, 1, . . . , 1]T i.e. the vector with all unit entries. This can be seen as by definition, each row of the graphc Laplacian matrix sums to 0. Below a few examples with corresponding spectral properties

P6 C6

S6

Path graph Pn . Here L is a rank−2 perturbation of the 1D Poisson matrix. We can bound the eigenvalues of L using the following result Theorem 3.2.1 ([40]). Consider two symmetric matrices S1 , S2 ∈ Rn×n . Let λ1 ≥ λ2 ≥ · · · ≥ λn . Then the kth eigenvalue of the matrix S1 + S2 satisfies λk (S1 ) + λn (S2 ) ≤ λk (S1 + S2 ) ≤ λk (S1 ) + λ1 (S2 ) Corollary 3.2.2 (Bound on eigenvalues of L for Pn ). Let ALap be defined as in (3.1). We write L = ALap + D where D is a diagonal matrix with D(1, 1) = −1, D(n, n) = −1 and all other entries zero. Following Theorem 3.2.1 the kth eigenvalue of L is bounded by max(λk (ALap ) − 1, λn (ALap + λk (D)) ≤ λk (L) ≤ min(λk (ALap ), λ1 (A) − λk (D)). In particular the extreme eigenvalues are bounded by: λmin (ALap ) − 1 ≤ λmin (L) ≤ min(λmin (ALap ), λmax (ALap ) − 1), max(λmax (ALap ) − 1, λmin (ALap )) ≤ λmax (L) ≤ λmax (ALap )

Cycle graph Cn . Eigenvectors are Fourier modes; eigenvalues are λk = 2 − 2 cos(2πk/n). Jacobi eigenvalues are cos(2πk/n). Star graph Sn . The spectrum of the graph Laplacian is {0, 1(×n − 2), n}. For real graph algorithms one typically turns to more advanced methods. Still, the classical iterations remain a useful “microscope” for understanding graph structure. They also have clear interpretations from the perspective of classical graph theory as we shall see later. Example 3.2.3 (Jacobi = random walk on a path graph). Consider the path P4 with vertices {1, 2, 3, 4}. The adjacency matrix and degree matrix are  0 1 A= 0 0

1 0 1 0

0 1 0 1

 0 0 , 1 0

D = diag(1, 2, 2, 1).

36

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

The Jacobi iteration matrix is 

0

1 2 GJ = D−1 A =  0 0

1 0

0

1 2

0 1

1 2

0

 0 0  1 . 2 0

Interpretation. Each row of GJ gives the transition probabilities of a random walk: • From vertex 1, you move to vertex 2 with probability 1. • From vertex 2, you move to 1 or 3 with equal probability 1/2. • From vertex 3, you move to 2 or 4 with probability 1/2. • From vertex 4, you move to vertex 3 with probability 1. Starting from an initial vector x(0) = [1, 0, 0, 0]⊤ (all mass at vertex 1), the Jacobi updates are x(1) = GJ x(0) = [0, 12 , 0, 0]⊤ ,

x(2) = GJ x(1) = 12 [1, 0, 12 , 0]⊤ ,

and so on. The values spread along the path exactly like the probability distribution of a random walk or like heat diffusion on the network. Thus Jacobi iteration coincides with a synchronous random walk (or diffusion) on the graph. This interpretation is standard in spectral graph theory. Takeaways Jacobi iteration on Lx = b is equivalent to the evolution of a random walk on the underlying graph, with eigenvalues describing diffusion rates. Thus Jacobi corresponds to a synchronous diffusion process or a “lazy random walk” on the graph. Relation to spectral clustering. The slowest-decaying Jacobi modes are precisely the eigenvectors of the Laplacian associated with the smallest nonzero eigenvalues. These eigenvectors are smooth signals on the graph: they vary little across densely connected regions and change more abruptly across weak connections. In spectral graph theory this property is exploited to detect communities [38]. For instance, the Fiedler vector [10] (the eigenvector associated with the second-smallest Laplacian eigenvalue) typically has nearly constant values on different clusters and a sharp transition across the cut between them. Thresholding this vector provides a natural bipartition of the graph. Example 3.2.4 (Jacobi iteration and community detection). Consider the graph consisting of two triangles connected by a single edge: vertices {1, 2, 3} fully connected to each other, vertices {4, 5, 6} fully connected to each other, and one edge (3, 4) connecting the two groups.

2

1

5

4

3

6

37

3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS Graph Laplacian construction. The adjacency matrix is   0 1 1 0 0 0 1 0 1 0 0 0   1 1 0 1 0 0  A= 0 0 1 0 1 1 ,   0 0 0 1 0 1 0 0 0 1 1 0 and the degree matrix is D = diag(2, 2, 3, 3, 2, 2). Thus the Laplacian is  2 −1 −1 0 0 0 −1 2 −1 0 0 0   −1 −1 3 −1 0 0 . L=D−A= 0 0 −1 3 −1 −1   0 0 0 −1 2 −1 0 0 0 −1 −1 2 

Spectrum and Fiedler vector. For the L above, λ(L) ≈ { 0, 0.438447, 3, 3, 3, 4.561553 }, so λ2 ≈ 0.438447 is smallest nonzero. An associated Fiedler vector is vFiedler ≈ [−0.465, −0.465, −0.261, 0.261, 0.465, 0.465]⊤ . It has a sign change change across the bridge (3, 4). Thresholding at 0 (or performing a sweep cut) yields the two clusters.

Jacobi as both power iteration and low-pass filter. Let G = (V, E) be an undirected graph with adjacency A, degree matrix D, and random walk matrix P := D−1 A. Define the symmetric normalized adjacency S := D−1/2 AD−1/2 and the normalized Laplacians Lrw := I − P and Lsym := I − S. We have the similarity P = D−1/2 SD1/2 , hence P and S share eigenvalues 1 = µ1 ≥ µ2 ≥ · · · ≥ µn ≥ −1. Let {vi }ni=1 be an orthonormal basis of eigenvectors of S 1/2 1 (Svi = µi vi ), so v1 = ∥DD1/2 1∥ and λsym = 1 − µi are the eigenvalues of Lsym . i 2

Low–pass view. Jacobi with b = 0 is x(k+1) = P x(k) , so with y (k) := D1/2 x(k) , y

(k)

k (0)

= S y

=

n X i=1

µik (y (0) , vi ) vi

=⇒

x

(k)

= D

−1/2

n X

µik (D1/2 x(0) , vi ) vi .

i=1

Thus each spectral component is multiplied by |µi | k = |1 − λsym | k : modes with large λsym (high i i “graph frequency”, i.e., oscillatory across edges) are attenuated rapidly, while smooth modes (small λsym ) persist. Jacobi is therefore a low-pass graph filter (i.e. a filter which removes high-frequency i components and permits low-frequency components to pass through). Power method view The low-pass picture has a classical counterpart. Let M ∈ Rn×n have a dominant eigenvalue, that is, let its eigenvalues be ordered as |λ1 | > |λ2 | ≥ · · · ≥ |λn |, so that ρ(M ) = |λ1 |. Multiplying a vector repeatedly by M amplifies its component along the dominant eigenvector v1 more than any other, so that after rescaling at each step the iterate lines up with v1 . This is the power method, the oldest of all eigenvalue iterations, and the simplest predecessor of the Lanczos method of Chapter 4.

38

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

Algorithm 1 Power method Require: Matrix M ∈ Rn×n with a dominant eigenvalue, starting vector y (0) with ∥y (0) ∥2 = 1 and a nonzero component along v1 , tolerance ε > 0. 1: for k = 0, 1, 2, . . . do 2: z ← M y (k) {one matrix–vector product per step} 3: µk ← (y (k) )⊤ z {Rayleigh quotient, since ∥y (k) ∥2 = 1} 4: if ∥z − µk y (k) ∥2 ≤ ε then 5: stop 6: end if 7: y (k+1) ← z/∥z∥2 8: end for Ensure: Approximate dominant eigenpair (µk , y (k) ).

Provided y (0) has a nonzero component along v1 , the iterates converge (up to sign) to the dominant eigenvector, and the Rayleigh quotients µk converge to λ1 . The rate is governed by the ratio of the two largest eigenvalue moduli: the eigenvector error decays like |λ2 /λ1 |k and the eigenvalue error, for symmetric M , twice as fast. Slow convergence is therefore expected when the two largest eigenvalues are not well separated. See Problem 8 for details. Deflating the trivial mode to reveal v2 . For the random walk, M = P has dominant eigenvalue 1 with right eigenvector 1 (since P 1 = 1). Because of this and in order to find the first eigenpair different from (1, 1) (and identify the Fiedler vector), we project it out after each step and normalize in the D-inner product ⟨u, v⟩D := u⊤ Dv:

x̃(k+1) = P x(k) ,

x(k+1) = x̃(k+1) −

⟨x̃(k+1) , 1⟩D 1, ⟨1, 1⟩D

∥x(k+1) ∥D = 1.

P P (k) = k Setting y (k) := D1/2 x(k) and expanding y (0) = i≥2 αi vi gives y i≥2 αi µi vi , hence x(k) → ±D−1/2 v2 at rate |µ3 /µ2 | k . Thus Jacobi = power iteration on P (with deflation) converging to the Fiedler vector. Classical point of view This viewpoint is standard in spectral graph theory and network science: P = D−1 A is the random-walk operator; S diagonalizes it; Jacobi is power iteration on P ; and the prominence of the Fiedler vector underlies spectral clustering. Takeaways • Classical iterations such as Jacobi act as diffusion or averaging processes on graphs. • Their eigenmode analysis shows that high–frequency oscillations on the graph are quickly damped, while smooth global modes persist. • The slow modes correspond to eigenvectors of the Laplacian with small eigenvalues; these capture global connectivity patterns and reveal community structure. • This perspective connects basic iterative methods to spectral clustering: Jacobi iterations naturally enhance low–frequency information that underlies graph partitioning. • In practice their convergence is too slow to serve as efficient graph solvers, but they remain conceptually valuable as “microscopes” for understanding diffusion and spectral structure.

39

3.3. BASIC ITERATIVE METHODS IN MACHINE LEARNING

3.3

Basic Iterative Methods in Machine Learning

The linear algebra perspective on Jacobi, Gauss–Seidel, and SOR extends naturally to optimization and machine learning. At their core, many ML algorithms amount to solving large linear or least-squares problems via iterative schemes. This section shows how the classical iterations appear as special cases of gradient descent, and how their spectral behaviour connects to convergence properties of training.

3.3.1

Richardson iteration and gradient descent

Recall the Richardson scheme for Ax = b with step size τ > 0: x(k+1) = x(k) + τ (b − Ax(k) ). This is exactly gradient descent for the quadratic loss min f (x) = 12 x⊤ Ax − b⊤ x, x

since ∇f (x) = Ax − b. Thus gradient descent on quadratic objectives is nothing but Richardson iteration. Spectral behaviour. If A is symmetric positive definite with eigenpairs (λj , vj ), then the error evolves as X e(k) = (1 − τ λj )k αj vj . j

Hence convergence requires 0 < τ < 2/λmax , and the rate is controlled by the condition number κ(A) = λmax /λmin [27]: the closer eigenvalues are to zero, the slower their modes decay. This is the same “bottleneck” phenomenon seen for smooth modes in PDEs and graph Laplacians. mode amplitude large λj (fast)

small λj (slow)

iteration k

Choosing the step size τ . Each mode j is damped at rate gj (τ ) = |1 − τ λj |, so the overall (worst-case) rate is q(τ ) = maxj gj (τ ), which for λj ∈ [λmin , λmax ] reduces to q(τ ) = max{|1−τ λmin |, |1−τ λmax |}. This is exactly the Richardson-for-SPD-matrices result of Chapter 1: the minimax-optimal step τ⋆ =

2 λmin + λmax

equalizes the two extreme mode gains and gives q(τ ⋆ ) =

κ−1 λmax ,κ = . κ+1 λmin

Example 3.3.1 (Optimal step size for a two-eigenvalue toy problem). Suppose A (or X ⊤ X) has only two distinct eigenvalues, λmin = 1 and λmax = 9, so κ = 9. The optimal step is τ ⋆ = 2/(1 + 9) = 0.2, giving mode gains |1 − τ ⋆ λmin | = |1 − 0.2| = 0.8 and |1 − τ ⋆ λmax | = |1 − 1.8| = 0.8: both extreme modes decay at exactly the same rate q ⋆ = (κ − 1)/(κ + 1) = 8/10 = 0.8, as predicted. A poorly chosen step, e.g. τ = 2/λmax = 2/9 ≈ 0.222, leaves the λmax mode exactly on the boundary of divergence (g = |1 − 2| = 1) while the λmin mode still decays comfortably (g = |1 − 0.222| ≈ 0.778) – a reminder that the largest eigenvalue alone does not determine the best step size; both extremes of the spectrum must be balanced.

40

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

Connection to Jacobi. Jacobi can be viewed as Richardson with a preconditioner D−1 : x(k+1) = x(k) + D−1 (b − Ax(k) ). This interpretation carries over directly to ML: preconditioning accelerates gradient descent by scaling the update direction.

3.3.2

Applications in machine learning

(1) Linear regression. Given data X ∈ Rm×n and labels y ∈ Rm , the least squares problem min 12 ∥Xw − y∥2 w

leads to the normal equations X ⊤ Xw = X ⊤ y. Richardson/gradient descent updates are
w(k+1) = w(k) − τ X ⊤ Xw(k) − X ⊤ y . The spectrum of X ⊤ X (squared singular values of X) dictates the convergence speed. If X is ill-conditioned, convergence is slow, mirroring the PDE case where low frequencies persist. This motivates preconditioning and methods like stochastic gradient descent (SGD). (2) Logistic regression and nonlinear losses. For convex but nonlinear losses, e.g. f (w) =

m X

log 1 + exp(−yi x⊤ i w) ,

i=1

gradient descent is no longer exactly Richardson, but locally (near the solution) the Hessian H = ∇2 f (w∗ ) plays the role of A. Convergence is again governed by the spectrum of H: flat directions (small eigenvalues) cause slow learning, steep directions (large eigenvalues) require small step sizes for stability.

3.3.3

Smoothing and spectral bias

Just as Jacobi smooths high-frequency error modes in PDEs, gradient descent smooths highfrequency parameter components. In training neural networks this is visible as spectral bias [31]: early iterations capture low-frequency/global patterns in the data, while fine details emerge later. Example 3.3.2 (Smoothing vs. spectral bias). Consider two settings: (a) Jacobi for the 1D Poisson problem. Let A be the Laplacian matrix on n grid points. The Jacobi error iteration reads   (k) (0) πj . e(k+1) = GJ e(k) , αj = µkj αj , µj = cos n+1 The eigenvectors are sine waves with frequency j. Large j corresponds to oscillatory modes, for which |µj | is small, so they are quickly damped. Smooth low–frequency modes (small j) have |µj | ≈ 1 and decay slowly. Thus Jacobi acts as a low–pass smoother : it removes oscillations first and leaves global structure last. (b) Gradient descent for least squares regression. Suppose we fit Xw ≈ y with normal equations X ⊤ Xw = X ⊤ y. Let X ⊤ X = V ΛV ⊤ with eigenpairs (λj , vj ). Gradient descent with step size τ yields (k) (0) αj = (1 − τ λj )k αj . For many data distributions (e.g. Fourier features, random features, kernel regression), it turns out that low-frequency patterns correspond to larger eigenvalues of the kernel or data covariance,

41

3.3. BASIC ITERATIVE METHODS IN MACHINE LEARNING

while high-frequency patterns correspond to smaller eigenvalues. This phenomenon is known as spectral bias. (k)

(0)

Duality. Both Jacobi and gradient descent share the same filtering mechanism: αj = g(λj )k αj , where λj indexes spectral modes. In PDE solvers, large λj corresponds to high–frequency error and is damped quickly (smoothing). In learning, large λj corresponds to low–frequency patterns in the data, which are learned quickly. Thus the roles of low vs. high frequency are swapped, but the underlying principle is the same: iterative methods act as frequency filters, with the spectrum of A or X ⊤ X dictating which modes persist. Takeaways • Richardson iteration is gradient descent on quadratic problems; Jacobi is Richardson with a diagonal preconditioner. • Convergence speed depends on the eigenvalue spectrum: directions corresponding to small eigenvalues decay slowly, leading to bottlenecks. • In ML problems this explains the role of conditioning, the need for preconditioning/normalization, and the phenomenon of spectral bias. • Classical iterative methods thus provide a unifying framework: from PDE solvers to optimization and machine learning, the spectral mechanisms are the same.

42

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

3.4

Exercises

Problem 1 (Closed form of the 1D Gauss–Seidel iteration matrix). Let 2 tridiag(−1, 2, −1) ∈ Rn×n , A = D − L − U, h2 where the Gauss-Seidel splitting uses nonnegative strict parts A=

D=

4 I, h2

L=

2 tridiag(1, 0, 0), h2

U=

2 tridiag(0, 0, 1). h2

The Gauss–Seidel iteration matrix is GGS = (D − L)−1 U . 1. Show that GGS is independent of h and can be written as GGS = B −1 S+ ,

B := 2I − S− ,

where S− has ones on the subdiagonal and S+ has ones on the superdiagonal. 2. Compute B −1 explicitly and deduce a closed form for the entries of GGS . 3. Write the first few rows of GGS to display its “plain matrix” pattern. Problem 2 (Error components for Jacobi). Let A be the 1D Poisson matrix with n interior points and GJ the Jacobi iteration matrix. (a) Show that the error e(k) = x(k) − x∗ can be expanded in the eigenbasis (k)

e

=

n X

(k)

(k+1)

αj v (j) ,

αj

(k)

= µj αj ,

j=1

  jπ with eigenvalues µj = cos n+1 . (b) Which modes (low j or high j) have eigenvalues closest to 1? (c) Conclude that Jacobi damps oscillatory (high–frequency) error components quickly, but smooth (low–frequency) error components persist. Preamble (for Problems 3 and 4, exam 2021/2022). For a0 , a−1 , a1 ∈ R with a−1 · a1 > 0 the Toeplitz tridiagonal matrix   a0 a1  a−1 a0 a1     .. .. .. A=  ∈ Rn×n . . .    a−1 a0 a1  a−1 a0 has eigenvalues √ λp = a0 + 2sgn(a−1 ) a−1 a1 cos and associated eigenvectors sp =

q

a−1 a1

i



pπ n+1

 ,

p = 1, . . . , n,

 n ipπ sin n+1 . i=1

If A is symmetric (a−1 = a1 ), then the eigenvectors form an orthogonal basis. Moreover, for any matrix B, lim B k = 0 ⇐⇒ ρ(B) < 1, k→∞

where ρ(B) denotes the spectral radius.

43

3.4. EXERCISES

Problem 3 (Tridiagonal Laplace matrix, exam 2021/2022). Consider the symmetric tridiagonal matrix   2 −1 −1 2 −1      .. .. .. A=  ∈ Rn×n . . . .    −1 2 −1 −1 2 Assume n is odd. (a) Show that the vector b = [1, 0, −1, 0, 1, 0, −1, . . . , 0, 1]⊤ ∈ Rn is an eigenvector of A. (b) Show that all eigenvalues of A lie in the interval (0, 4). Problem 4 (Jacobi’s method and convergence, exam 2021/2022). Let A ∈ Rn×n be nonsingular and let D = diag(A). Jacobi’s iteration for solving Ax = b is x(k+1) = (I − D−1 A)x(k) + D−1 b,

k = 0, 1, 2, . . .

(3.3)

Let B = I − D−1 A. (a) Show that if x̂ = limk→∞ x(k) exists, then Jacobi’s method converges if and only if A has no zeros on its diagonal and ρ(B) < 1. (b) Suppose A is the Toeplitz tridiagonal matrix above. Give necessary and sufficient conditions on (a−1 , a0 , a1 ) independent of n for Jacobi’s method to converge. (c) Consider the special case  a0 a0   A=  

a1 a0 .. .

 a1 .. . a0

   ..  ∈ Rn×n . .  a0 a1  a0 a0

What are the necessary and sufficient conditions on a0 , a1 (independent of n) for ρ(B) < 1? (d) Let A be as in (c) with a0 = 1 and a1 = ̸ 0, b = en , and initial guess x(0) = 0. Does there (k) exist a Jacobi iterate x that solves Ax = b exactly? Problem 5 (Numerical verification of convergence rates). For the 1D Poisson matrix with n = 19 interior points: (a) Compute ρ(GJ ), ρ(GGS ), and ωopt , ρ(Gωopt ) using the closed-form expressions derived in this chapter. (b) Estimate the number of iterations each method needs to reduce the initial error by a factor 10−8 . (c) Implement Jacobi, Gauss–Seidel, and SOR (with ω = ωopt ) for this system and confirm numerically that the observed asymptotic convergence factors match your predictions from (a). (See the course GitHub repository for a Python starting point.)

44

CHAPTER 3. APPLICATIONS OF BASIC ITERATIVE METHODS

Problem 6 (Spectrum of a small path graph). Let P5 be the path graph on 5 vertices with graph Laplacian L (as constructed in Chapter 2). (a) Write down L explicitly and compute its eigenvalues and eigenvectors directly (by hand or numerically). (b) Verify that they match the general formula λk (L) = 2 − 2 cos(πk/n) for the path graph. (c) Identify the Fiedler vector (the eigenvector for the smallest nonzero eigenvalue) and describe the bipartition it induces on the 5 vertices. Problem 7 (Gradient descent convergence rate from the data spectrum). Consider a least-squares problem minw 12 ∥Xw − y∥2 where X ∈ Rm×n has singular values σ1 = 4, σ2 = 2, σ3 = 1 (so X ⊤ X has eigenvalues 16, 4, 1). (a) Compute the condition number κ(X ⊤ X) and the largest admissible step size τmax = 2/λmax for gradient descent on X ⊤ Xw = X ⊤ y. (b) Compute the minimax-optimal step size τ ⋆ and the resulting worst-case convergence factor q ⋆ = (κ − 1)/(κ + 1). (c) Using τ ⋆ , estimate how many iterations are needed to reduce the worst-case mode’s error by a factor of 100. (d) Which mode (associated with σ1 , σ2 , or σ3 ) is learned fastest, and which slowest? Relate your answer to the notion of spectral bias. Problem 8 (Convergence of the power method). Let M ∈ Rn×n be symmetric, with orthonormal eigenvectors v1 , . P . . , vn and eigenvalues ordered by modulus, |λ1 | > |λ2 | ≥ · · · ≥ |λn |. Run Algorithm 1 from y (0) = ni=1 αi vi with ∥y (0) ∥2 = 1 and α1 ̸= 0, and let θk ∈ [0, π/2] be the angle between y (k) and the line spanned by v1 . (a) Show that rescaling at every step is the same as rescaling once at the end, i.e. y

(k)

M k y (0) = , ∥M k y (0) ∥2

so

k (0)

M y

=

n X

αi λik vi .

i=1

(b) Splitting this sum into its component along v1 and the rest, deduce that
P 2 1/2 λ2 k i≥2 αi tan θk ≤ , |α1 | λ1 so that y (k) → ±v1 geometrically with ratio |λ2 /λ1 |. P (c) Writing ci = (y (k) )⊤ vi , show that µk = i λi c2i and hence that |µk − λ1 | ≤ (λmax − λmin ) sin2 θk . Conclude that the eigenvalue converges at rate |λ2 /λ1 |2k , twice as fast as the eigenvector. (d) Take the 1D Poisson   matrix A = tridiag(−1, 2, −1) with n = 19, whose eigenvalues are λj = 2 − 2 cos jπ 20 . Identify the two largest in modulus, compute the ratio |λ2 /λ1 |, and estimate how many power iterations are needed to reduce the eigenvector error by a factor 10−2 . Relate the answer to the discussion of smoothing in Section 3.1. (e) The random-walk matrix P = D−1 A of a connected graph has dominant eigenvalue 1 with eigenvector 1. Explain what the power method converges to in that case, and why projecting out 1 at every step (as in this chapter) makes it converge to the Fiedler vector instead.

Chapter 4

Krylov Methods for Symmetric Systems Overview Chapters 1 and 3 built iterative solvers from a fixed matrix splitting A = M − N : every step applies the same operator G = M −1 N , so the error can be reduced by at best a fixed factor ρ(G) per iteration. For the PDE model matrices of Chapter 3 that factor behaves like 1 − O(h2 ), which makes stationary methods unusually slow as the mesh is refined. This chapter replaces that fixed recipe by an adaptive one. At step k, a Krylov subspace method relies on the subspace Kk (A, r(0) ) = span{r(0) , Ar(0) , . . . , Ak−1 r(0) } and then selects the best approximation in the affine space x(0) + Kk (A, r(0) ). For symmetric positive definite (SPD) A, this optimal choice can be computed with short recurrences and constant p storage: this is the Conjugate Gradient (CG) method [21], and its convergence is governed by κ(A) instead of κ(A). This chapter covers: (1) CG, its optimality properties and its Chebyshev convergence bound; (2) the Lanczos algorithm, including the fact that CG and Lanczos are one recurrence seen from two sides, and its use to compute the Fiedler vector of a graph; (3) machine learning, where we show that ridge regression, truncated SVD and CG are all spectral filters of the same data, so that stopping CG early is itself a regularisation. Learning objectives By the end of this chapter, you should be able to: • Derive the Conjugate Gradient method from the minimisation of the quadratic functional J(x) = 12 x⊤ Ax − b⊤ x, and explain why the steepest descent is not enough. • State and use the optimality property of CG: the iterate x(k) minimises the A-norm of the error over x(0) + Kk (A, r(0) ). p • Use the Chebyshev bound to relate the CG convergence rate to κ(A), and contrast this with the κ(A) dependence of steepest descent and stationary methods. • Describe Lanczos iteration and explain why Ritz values converge at the ends of the spectrum first, and use deflation to compute the Fiedler vector of a graph Laplacian. • Explain the CG-Lanczos equivalence, and interpret ridge regression and truncated CG as spectral filters - so that early stopping is itself a regularisation. • Explain what semi-convergence is, and why some stopping rule is still needed when the data are noisy.

45

46

4.1

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

The Conjugate Gradient method

We now move from stationary iterations such as Jacobi, Gauss–Seidel, and SOR to Krylov subspace methods. These methods generate successively richer approximation spaces and the most famous example is the Conjugate Gradient (CG) method. For symmetric positive definite (SPD) matrices, such as those arising from elliptic PDEs, CG achieves convergence far superior to classical schemes. We wish to solve Ax = b with A ∈ Rn×n a SPD matrix. In this case, it makes sense to consider the quadratic functional J(x) = 12 x⊤ Ax − b⊤ x. Minimising J(x) is equivalent to solving ∇J = Ax − b = 0 hence the linear system Ax = b. Definition 4.1.1 (Steepest descent method). Starting from an initial guess x0 , the steepest descent method generates iterates x(k+1) = x(k) + τk r(k) ,

r(k) = b − Ax(k) ,

where the step size τk is chosen to minimize J along the search direction opposite to the gradient of J i.e. r(k) . Let’s compute the formula of the step
length in the steepest descent method. With the notations above, let ϕ(τ ) = J x(k) + τ r(k) . Then  ⊤  ⊤  ⊤
ϕ′ (τ ) = r(k) Ax(k) − b + τ r(k) Ar(k) = − ∥r(k) ∥2 + τ r(k) Ar(k) . Setting ϕ′ (τ ) = 0 gives τk =

r(k)⊤ r(k) ∥r(k) ∥2 = , r(k)⊤ Ar(k) ∥r(k) ∥2A

which is equivalent to minimising J along the direction r(k) . Proposition 4.1.2 (Convergence of steepest descent). Let x∗ be the exact solution of Ax = b and let κ(A) = λmax (A)/λmin (A) denote the condition number of A. Then the steepest descent iterates satisfy   κ(A) − 1 k (0) (k) ∗ ∥x − x ∥A ≤ ∥x − x∗ ∥A , κ(A) + 1 where ∥v∥2A = v ⊤ Av (norm induced by A). Sketch of proof. Let e(k) = x(k) − x∗ and r(k) = b − Ax(k) = −Ae(k) . A steepest descent step gives e(k+1) = e(k) + τk r(k) . Expanding the A-norm of the error yields ∥e(k+1) ∥2A = ∥e(k) ∥2A − τk ∥r(k) ∥2 , where we used the optimal choice of τk . Dividing by ∥e(k) ∥2A and substituting the formula of gives ∥e(k+1) ∥2A ∥r(k) ∥4 = 1 − . ∥e(k) ∥2A (r(k)⊤ Ar(k) )(r(k)⊤ A−1 r(k) )

47

4.1. THE CONJUGATE GRADIENT METHOD

The denominator is bounded below via the Kantorovich inequality 1 , which relates (v ⊤ Av)(v ⊤ A−1 v) to the condition number κ(A). Inserting this bound yields κ(A) − 1 ∥e(k+1) ∥A ≤ . (k) κ(A) + 1 ∥e ∥A By induction over k, the stated result follows. Remark 4.1.3. In particular, for PDE discretizations of the Poisson equation introduced in Chapter 2 −u′′ = f,

u(0) = u(1) = 0,

which yields the SPD stiffness matrix A ∈ Rn×n . Here, the condition number scales as κ(A) = O(n2 ). Hence steepest descent requires O(n2 ) iterations, which becomes prohibitive for large systems. Steepest descent

Conjugate Gradient x(0)

x(0) x∗

x∗

p(1)

p(0)

two A-conjugate directions: exact solution after n = 2 steps

successive residuals are orthogonal: the iterates zig-zag down the valley

Figure 4.1: Steepest descent versus Conjugate Gradient on the level sets of J(x) = 12 x⊤ Ax − b⊤ x with A = 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∗ ; the number of steps grows like κ(A), so the picture gets dramatically worse as the contours become more elongated. CG replaces orthogonality of the residuals by A-conjugacy of the search directions: no direction is ever revisited, and in exact arithmetic n steps suffice.

When Hestenes and Stiefel introduced Conjugate Gradient in 1952, computing was still in its infancy: machines like the SWAC in Los Angeles (then fastest computer in the world) had memory measured in kilobytes and could barely store a few thousand numbers. Steepest descent was well known, but its zig-zagging inefficiency was frustrating. The breakthrough of CG was to design successive directions that “remember” past progress and never waste effort repeating the same zig-zag. This made it one of the first practical Krylov methods and a cornerstone of numerical linear algebra. The idea of Conjugate Gradient (CG), developed by Hestenes and Stiefel [21], is to construct search directions p(k) that are not only descent directions but are A-orthogonal. Using such conjugate directions ensures that each step eliminates error components independently, guaranteeing convergence in at most n steps in exact arithmetic. 1

⊤

⊤

A See the exercise list: 1 ≤ (v Av)(v ∥v∥4

−1

v)

≤ 14 (

p κ(A) + √ 1

κ(A)

)2

48

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

Definition 4.1.4 (Conjugate directions). Let A ∈ Rn×n be SPD. Two vectors p, q ∈ Rn are said to be A-conjugate if p⊤ Aq = 0. A family of vectors {p(0) , . . . , p(k) } is called mutually A-conjugate if every pair is A-conjugate. Such directions play the role of “orthogonal axes” with respect to the A-inner product ⟨u, v⟩A = u⊤ Av. Definition 4.1.5 (Krylov subspaces). Given an initial residual r(0) = b − Ax(0) , the k-th Krylov subspace is Kk (A, r(0) ) = span{r(0) , Ar(0) , . . . , Ak−1 r(0) }. Equivalently,

Kk (A, r(0) ) = { q(A)r(0) : q ∈ Πk−1 },

where Πk−1 denotes the set of polynomials of degree at most k − 1. Definition 4.1.6 (Conjugate Gradient method). Let Ax = b with A ∈ Rn×n be SPD and the right-hand side b ∈ Rn . Starting from an initial guess x(0) and the initial descent direction p(0) = r(0) , the Conjugate Gradient (CG) generates a sequence of approximations {x(k) } by choosing search directions p(k) ∈ Kk+1 (A, r(0) ) iteratively as follows: x(k+1) = x(k) + αk p(k) ,

αk =

r(k)⊤ r(k) , p(k)⊤ Ap(k)

r(k+1) = r(k) − αk Ap(k) ,

βk =

r(k+1)⊤ r(k+1) , r(k)⊤ r(k)

p(k+1) = r(k+1) + βk p(k) . Proposition 4.1.7 (Properties of CG). Let A be SPD. Then the iterates of CG satisfy: • Residuals are mutually orthogonal: r(i)⊤ r(j) = 0 for i ̸= j. • Search directions are A-orthogonal: p(i)⊤ Ap(j) = 0 for i ̸= j. • x(k) minimizes ∥x − x∗ ∥A over the affine space x(0) + Kk (A, r(0) ).

x(0) + Kk (A, r(0) ) x(k)

x∗

x(0)

 x(k) = arg min ∥x − x∗ ∥A : x ∈ x(0) + Kk (A, r(0) )

Figure 4.2: The optimality property of CG. The ellipses are level sets of the energy norm ∥x−x∗ ∥A . At step k the iterate is the point of the affine space x(0) + Kk (A, r(0) ) (dashed line) closest to x∗ in the A-norm, i.e. the A-orthogonal projection of x∗ onto it.

49

4.1. THE CONJUGATE GRADIENT METHOD

Proposition 4.1.8 (Polynomial interpretation and finite termination of CG). Let A ∈ Rn×n be SPD, x∗ denote the exact solution of Ax = b, and e(k) = x(k) − x∗ the CG error at step k (with x(0) arbitrary). 1. ( Polynomial representation) For each k ≥ 0 there exists a polynomial qk with deg qk ≤ k and qk (0) = 1 such that e(k) = qk (A) e(0) . The polynomial qk minimizes the A norm of the error at step k over the set  q(A)e(0) : q ∈ Πk , q(0) = 1 . 2. ( Finite termination) Let m be the minimal polynomial of A and let d = deg m (for symmetric A, d equals the number of distinct eigenvalues). Then CG attains the exact solution in at most d ≤ n steps e(d) = 0. Proof. (1) CG produces x(k) ∈ x(0) + Kk (A, r(0) ), hence e(k) = x(k) − x∗ = e(0) − vk ,

vk ∈ Kk (A, r(0) ).

Write vk = A sk−1 (A)e(0) with deg sk−1 ≤ k − 1 (since Kk (A, r(0) ) = A Kk (A, e(0) )). Then
e(k) = I − A sk−1 (A) e(0) = qk (A)e(0) , where qk (λ) = 1 − λsk−1 (λ) has degree ≤ k and satisfies qk (0) = 1. This also shows e(k) ∈ Kk+1 (A, e(0) ) = span{e(0) , Ae(0) , . . . , Ak e(0) }. The minimization follows from Prop. 4.1.7. (2) Let m(λ) be the minimal polynomial of A, so m(A) = 0 and deg m = d. Define qd (λ) =

m(λ) . m(0)

Then qd (0) = 1, qd (A) = 0, hence qd (A)e(0) = 0. By part (1), the set of attainable errors after k ≥ d steps includes 0, so the CG minimizer in the A-norm satisfies ∥e(d) ∥A = 0, i.e., e(d) = 0. Q For symmetric A with eigenvalues {λ1 , . . . , λn } and d distinct values, m(λ) = dj=1 (λ − λj ), so d ≤ n, with equality when all eigenvalues are distinct. Remark: All the above is true in exact arithmetic. In floating-point arithmetic, loss of orthogonality may delay (or prevent) exact termination; the classical analysis of this effect is due to Paige [30]. Proposition 4.1.9 (Convergence of Conjugate Gradient). Let x∗ = A−1 b denote the exact solution and κ(A) = λmax (A)/λmin (A) the condition number of A. Then the CG iterates satisfy the error bound !k p κ(A) − 1 ∥x(k) − x∗ ∥A ≤2 p , ∥x(0) − x∗ ∥A κ(A) + 1 where ∥v∥2A = v ⊤ Av. In exact arithmetic the method terminates in at most n steps.

Example. For the 1D Poisson matrix of size n, the Steepest descent needs on the order of n2 iterations, while CG requires only n iterations. This dramatic acceleration is what makes CG the standard choice for large SPD PDE systems. A full treatment of CG and its convergence theory can be found in [13, 33].

50

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS ∥e(k) ∥A /∥e(0) ∥A

SD: √

CG: 2 √κ−1 κ+1

κ−1 k κ+1

k

k

Figure 4.3: Error decay of steepest descent and CG for the same p SPD matrix. Both are linear (geometric) rates, but CG’s contraction factor is governed by κ(A) rather than κ(A). For the 1D Poisson matrix, where κ(A) = O(n2 ), this is the difference between O(n2 ) and O(n) iterations. The bound is only an upper estimate: when the eigenvalues of A are clustered, CG typically converges considerably faster than the curve shown. Takeaways • Comparison with classical iterations. – Jacobi, Gauss–Seidel, and SOR are simple stationary schemes with iteration counts scaling like O(n2 ) for PDEs. – Steepest descent improves the picture but still requires O(κ(A)) iterations, which is prohibitive when κ(A) scales like n2 . p – Conjugate Gradient reduces the required iterations to O( κ(A)), leading to dramatic savings in practice. • Finite termination property. In exact arithmetic CG produces the exact solution in at most n steps, unlike stationary methods which only converge asymptotically. • Energy minimization viewpoint. Each CG iterate minimizes the quadratic energy functional over a growing Krylov subspace, making the method “optimal” among all polynomial approximations of degree k.

Method Jacobi Gauss–Seidel SOR (optimal ω) Steepest descent Conjugate Gradient

Convergence
k 1 − cn−2
k 1 − cn−2
k 1 − cn−1
κ−1 k 2

κ+1 √ k √κ−1 κ+1

Iterations (PDE)

Remarks

O(n2 )

Very slow; mainly used as smoother

O(n2 )

Roughly twice as fast as Jacobi

O(n)

Improved, but still asymptotically slow

O(κ) √ O( κ)

Gradient method; zig-zag convergence Optimal Krylov solver

Table 4.1: Comparison of classical stationary methods, steepest descent, and Conjugate Gradient for SPD PDE systems.

4.2

The Lanczos Method and the Symmetric Eigenproblem

So far, the Krylov subspace Km (A, r(0) ) has been a useful tool for solving Ax = b. The same subspace answers a different question as well: what does the spectrum of A look like? There are two different aspects to this question:

4.2. THE LANCZOS METHOD AND THE SYMMETRIC EIGENPROBLEM

51

• We need κ(A). Every bound in Section 4.1 is stated in terms of λmin (A) and λmax (A), yet for a large sparse A these are exactly the quantities we cannot compute directly. Estimating them is not a side task: it is what tells us whether CG will take ten iterations or ten thousand, and, from Chapter 6, whether a preconditioner is doing its job. • The spectrum is needed but also the eigenvectors. In Chapter 3 we saw that for a graph Laplacian L = D − A the eigenvector belonging to the smallest nonzero eigenvalue, i.e. the Fiedler vector [10], encodes how the graph splits into communities. In both cases n is far too large for a full eigendecomposition, which costs O(n3 ) and destroys sparsity. What we can afford, as always, is the product Av. Instead, we rely on iterative Krylov subspace methods such as the Lanczos algorithm [23]. The Lanczos algorithm is an efficient iterative method for approximating eigenvalues and eigenvectors of a large symmetric matrix A. The key idea is: • Start with a normalized vector q1 . • Build the Krylov subspace Km (A, q1 ) = span{q1 , Aq1 , A2 q1 , . . . , Am−1 q1 }. • Iteratively construct an orthonormal basis Qm = [q1 , . . . , qm ] of Km (A, q1 ) and project: Tm = Q⊤ m AQm . • For a general matrix A, the projected matrix Tm is upper Hessenberg, i.e. upper triangular with one extra subdiagonal. This is the Arnoldi case and is the subject of Chapter 5. • When A is symmetric, however, ⊤ ⊤ ⊤ ⊤ ⊤ Tm = (Q⊤ m AQm ) = Qm A Qm = Qm AQm = Tm ,

so Tm is symmetric as well. Being simultaneously Hessenberg and symmetric, it must be tridiagonal :   α1 β1   ..  β1 α2  .  . Tm =   .. ..  . . βm−1  βm−1 αm • Only the scalars αj and βj therefore need to be computed, and each new basis vector is orthogonal to all but two of its predecessors automatically. This is the Lanczos algorithm, displayed in Algorithm 2. • Eigenvalues of Tm (Ritz values) approximate eigenvalues of A, and Qm y approximates an eigenvector of A if y is an eigenvector of Tm . Notation: two different α’s and β’s. The symbols α and β are unfortunately standard in both algorithms, with unrelated meanings. In CG (Definition 4.1.6) αk is a step length and βk a direction-update coefficient; in Lanczos (Algorithm 2) αj and βj are the diagonal and off-diagonal entries of Tm . Both conventions are too well established to be worth changing, so we keep them. In Section 4.2.1, where the two algorithms are compared side by side, we write αkcg , βkcg for the CG scalars and reserve the αj , βj for the entries of Tm .

52

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

Algorithm 2 Lanczos Algorithm Require: Symmetric matrix A ∈ Rn×n , starting vector q1 with ∥q1 ∥2 = 1, number of steps m. 1: Set β0 = 0, q0 = 0. 2: for j = 1, . . . , m do 3: w ← Aqj − βj−1 qj−1 4: αj ← qj⊤ w 5: w ← w − αj qj 6: βj ← ∥w∥2 7: if βj = 0 then 8: stop 9: else 10: qj+1 ← w/βj 11: end if 12: end for 13: Construct the tridiagonal matrix Tm .

A ∈ Rn×n

Krylov basis m matvecs

n huge and sparse

−→

Qm ∈ Rn×m

Q⊤ m AQm tridiagonal

Tm ∈ Rm×m

m ≪ n, dense but tiny: eigenvalues by any direct method

Figure 4.4: The Lanczos idea in one line. The large sparse symmetric matrix A is never factorised; it is only applied to vectors. After m steps the problem has been compressed into an m × m symmetric tridiagonal matrix Tm whose eigenvalues (the Ritz values) approximate the extreme eigenvalues of A, and whose eigenvectors lift back to approximate eigenvectors of A through Qm .

Interpretation. The matrix Tm ∈ Rm×m represents the orthogonal projection of A onto the Krylov subspace Km (A, q1 ), expressed in the orthonormal basis q1 , . . . , qm . Its eigenvalues, called Ritz values, provide approximations to the eigenvalues of A and, as the next proposition makes precise, to the extremal ones first. Since Tm is symmetric and tridiagonal, its eigenvalues can be computed efficiently and stably using standard algorithms [12]. Eigenvectors of Tm , lifted back via Qm = [q1 , . . . , qm ], yield approximations of the eigenvectors of A. Proposition 4.2.1 (Ritz values converge at the ends of the spectrum). Let A = A⊤ have eigenvalues λ1 ≤ · · · ≤ λn and let θ1 ≤ · · · ≤ θm be the eigenvalues of Tm = Q⊤ m AQm . Then λ1 ≤ θ 1 ,

θm ≤ λn ,

and more generally the θi interlace the λi (Problem 7). Moreover θ1 and θm are the smallest and largest values of the Rayleigh quotient RA (x) = x⊤ Ax/x⊤ x over Km (A, q1 ), so they improve monotonically with m. Two consequences matter in practice. First, the Ritz values approach λ1 and λn from the inside, so a handful of Lanczos steps already brackets the spectrum which is precisely the estimate of κ(A) that the CG bound of Section 4.1 needs. Second, interior eigenvalues are reached much more slowly. If the eigenvalue of interest is not extremal, one must deflate the eigenvectors that are in the way. The next example uses the first of these. Sharp convergence estimates for the Ritz values, which show that the rate again improves by a square root over the power method, are given in [32].

4.2. THE LANCZOS METHOD AND THE SYMMETRIC EIGENPROBLEM

4.2.1

53

Lanczos and CG are the same recurrence

It is no accident that Lanczos and CG appear in the same chapter. Both build an orthonormal basis of the same Krylov subspace by a three-term recurrence, and in exact arithmetic they build the same one. The following proposition makes this precise and, more usefully, shows that the entries of Tm can be read off from a CG run without any extra work. Proposition 4.2.2 (CG and Lanczos generate the same basis, [6, Section 3.1, Theorem 3]). Let A be SPD and b ̸= 0. Run CG (Definition 4.1.6) from x(0) = 0, so that r(0) = b, and assume r(j) ̸= 0 for j = 0, . . . , m−1; denote its scalars by αkcg , βkcg . Run Lanczos (Algorithm 2) on A with q1 = b/∥b∥2 , producing q1 , . . . , qm and the entries αj , βj of Tm . Then (i) the two methods produce the same orthonormal basis of Km (A, b), namely qj+1 = (−1)j

r(j) , ∥r(j) ∥2

j = 0, . . . , m − 1;

(ii) the entries of Tm are determined by the CG scalars,

α1 =

1 , α0cg

cg βj−2 1 αj = cg + cg αj−1 αj−2

(j ≥ 2),

βj =

q cg βj−1 cg αj−1

(j ≥ 1).

Proof. (i) By Proposition 4.1.7 the residuals r(0) , . . . , r(m−1) are mutually orthogonal, and r(j) ∈ Kj+1 (A, b) \ Kj (A, b). Once normalised they are therefore an orthonormal basis of Km (A, b) adapted to the nested sequence K1 ⊂ K2 ⊂ · · · . The Lanczos vectors are, by construction, the Gram–Schmidt basis of that same nested sequence, and such a basis is unique up to the sign of each vector. The signs are pinned down by the convention βj = ∥w∥2 > 0 in Algorithm 2, which produces the alternating factor (−1)j . (ii) Eliminate the search directions. From r(k+1) = r(k) − αkcg Ap(k) we get Ap(k) = r(k) −
cg r(k+1) /αkcg , and substituting this into Ar(k) = Ap(k) − βk−1 Ap(k−1) , which follows from p(k) = cg r(k) + βk−1 p(k−1) , gives the three-term recurrence Ar(k) = −

cg  cg βk−1 βk−1 r(k+1)  1 (k) (k−1) + + r − . cg cg r αkcg αkcg αk−1 αk−1

cg (k) (k) k Now q write ρk = ∥r ∥2 , so that r = (−1) ρk qk+1 by (i). The definition of βk gives ρk+1 /ρk = βkcg . Dividing the recurrence by (−1)k ρk turns every sign positive and yields

Aqk+1 =

q βkcg αkcg

qk+2 +

cg  βk−1 + cg qk+1 + αkcg αk−1

 1

q cg βk−1 cg αk−1

qk .

Comparing with the Lanczos recurrence Aqj = βj qj+1 + αj qj + βj−1 qj−1 and setting j = k + 1 gives the stated formulas. Remark 4.2.3. The alternating sign in (i) is a normalisation convention: replacing qj+1 by −qj+1 conjugates Tm by a diagonal matrix of signs, which changes the signs of the off-diagonal entries but leaves the Ritz values, and hence everything we use them for, unchanged. The practical payoff is that spectral information comes for free from a run one was performing anyway: a CG solve can report an estimate of λmin , λmax and hence κ(A) at the cost of a few scalar operations and no extra matrix–vector products. Proposition 4.2.2 also explains the finite-termination property of Section 4.1 from the other side: CG stops after d steps because

54

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

Lanczos exhausts the Krylov space after d steps, d being the number of distinct eigenvalues. A systematic account of this correspondence, together with the consequences it has for error estimation inside a running CG solve, is given in [6].

4.2.2

Application: computing the Fiedler vector

We return to the graph of Chapter 3: two triangles joined by a single edge,   0 1 1 0 0 0 1 0 1 0 0 0   1 1 0 1 0 0  , A= L = D − A, D = diag(2, 2, 3, 3, 2, 2).  0 0 1 0 1 1 0 0 0 1 0 1 0 0 0 1 1 0 There we ran Jacobi on this Laplacian and observed that the slowest-decaying error mode was the Fiedler vector. We can perform a similar computation with Lanczos. Deflating the known null vector. L is symmetric positive semidefinite and always has the eigenpair (0, 1), since every row sums to zero. The Fiedler value λ2 is therefore not extremal, but it sits just above a known eigenvalue. By the proposition above, Lanczos would converge to 0 and to λmax long before it resolved λ2 . The remedy is immediate: start from a vector orthogonal to 1 and keep the iteration in that subspace, q1 ⊥ 1,

w ← w − (1⊤ w) 1/n at every step.

On the deflated subspace λ2 is the smallest eigenvalue, and the proposition applies again. Example 4.2.4 (Lanczos on the two-triangle Laplacian). The exact spectrum is {0, 0.4384, 3, 3, 3, 4.5616}. Starting from a normalised q1 ⊥ 1 and deflating as above, the Ritz values are m 1 2 3

Ritz values θi 3.1323 2.4887, 3.7217 0.4384, 3.0000, 4.5616

After three steps the Ritz values are the exact nonzero eigenvalues, and the Ritz vector belonging to θ1 = 0.4384 is v = Q3 y1 = (0.4647, 0.4647, 0.2610, −0.2610, −0.4647, −0.4647)⊤ , the Fiedler vector to machine precision, and thresholding it is the simplest form of spectral clustering. Its sign pattern splits the vertices exactly as the two triangles, {1, 2, 3} against {4, 5, 6}, and the small entries ±0.2610 identify nodes 3 and 4 as the bottleneck across the cut. Takeaways • Lanczos is a Krylov projection: m matrix-vector products compress A into an m × m tridiagonal Tm , whose eigenvalues (Ritz values) approximate those of A. • Convergence happens at the ends of the spectrum first. Extremal eigenvalues come cheaply; a few Lanczos steps therefore bracket λmin , λmax hence κ(A), which the CG bound of Section 4.1 needs. • On graph Laplacians, deflating the constant vector makes the Fiedler value extremal, and Lanczos delivers the community structure in an optimal number of steps.

4.3. KRYLOV METHODS AND MACHINE LEARNING (OPTIONAL)

4.3

55

Krylov Methods and Machine Learning (optional)

Fitting a linear model to data by least squares is the simplest learning task there is, and one of the most widely used. The same problem appears well beyond machine learning: reconstructing a sharp image from a blurred one, recovering a physical parameter from indirect measurements, and combining a forecast with observations in the data assimilation systems of Chapter 2 all reduce to it. In each case the matrix is ill-conditioned and the data are noisy, so the point of view differs from that of the earlier sections: there the system had a unique solution and we wanted it accurately, whereas here solving to machine precision reproduces the noise rather than the underlying model. 2

4.3.1

Least squares and the SVD

We are given N observations, stored as the rows of a design matrix X ∈ RN ×d , together with the corresponding targets y ∈ RN , and we look for coefficients x ∈ Rd for which Xx is as close to y as possible. This is the linear least squares problem min ∥Xx − y∥22

X ⊤ X x = X ⊤ y,

⇐⇒

x∈Rd

(4.1)

the equivalence being obtained by setting the gradient of the objective to zero. The matrix X ⊤ X is symmetric positive semidefinite, and positive definite whenever X has full column rank, so we are back in the setting of this chapter. Proposition 4.3.1 (Least-squares solution via the SVD). Let X ∈ RN ×d have full column rank, with thin SVD X = U ΣV ⊤ and Σ = diag(σ1 , . . . , σd ), σ1 ≥ · · · ≥ σd > 0. Then the solution of (4.1) is unique and given by xls =

d X u⊤ y i

i=1

σi

vi .

(4.2)

Proof. Since U has orthonormal columns, X ⊤ X = V ΣU ⊤ U ΣV ⊤ = V Σ2 V ⊤ , which is invertible because every σi > 0; the normal equations in (4.1) therefore have the unique solution xls = (X ⊤ X)−1 X ⊤ y = V Σ−2 V ⊤ V ΣU ⊤ y = V Σ−1 U ⊤ y, which written columnwise is (4.2). (If rank X = r < d the solution is no longer unique, but replacing d by r in (4.2) gives the one of smallest norm.) Formula (4.2) also shows where the difficulty lies. The size of σi measures how much the data say about the direction vi . Suppose now the targets carry noise, y being replaced by y + e. −3 By (4.2) the coefficient i -th changes by u⊤ i e/σi : the same noise is divided by σi . If σi = 10 , a noise component of size 10−3 shifts that coefficient by 1. The coefficients the data determine the least well are thus the ones the noise affects the most, which is why xls is of little use when X is ill-conditioned.

4.3.2

Spectral filters and the truncated SVD

The remedy is modify the form of (4.2) while damping the terms responsible for amplification of errors. We therefore replace xls by x b=

d X i=1

2

fi

u⊤ i y vi , σi

fi ∈ [0, 1],

This section is not examinable but the ideas in it are likely to be useful for the project.

(4.3)

56

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

where the coefficients fi are called filter factors: fi = 1 retains the i-th direction exactly as least squares had it, fi = 0 discards it, and intermediate values damp it. Choosing the fi , keeping the directions the data determine well and suppressing the others, is what regularisation means [18]. The simplest choice is a sharp cutoff: for a fixed k ≤ d, set fi = 1 for i ≤ k and fi = 0 otherwise, keeping the k best-determined directions. This is the truncated SVD, the most transparent of the filters but also the most expensive: it requires the singular value decomposition of X, out of reach for large d, together with a cutoff k chosen by hand. The next two sections present filters that require neither.

4.3.3

Ridge regression: a smooth filter

Rather than truncating, one may penalise large coefficients. Ridge regression solves min 12 ∥Xx − y∥22 + λ2 ∥x∥22

⇐⇒

x

(X ⊤ X + λI) x = X ⊤ y,

where λ > 0 is a regularisation parameter. The matrix Aλ = X ⊤ X + λI is SPD, so CG applies directly and never requires Aλ itself: Aλ v = X ⊤ (Xv) + λv costs two passes over the data, and no d × d matrix is ever formed. This is what makes Krylov methods natural here. It remains to determine the effect of λ. Proposition 4.3.2 (Ridge regression is a spectral filter). Under the assumptions of Proposition 4.3.1 and for λ > 0, the ridge solution xλ = (X ⊤ X + λI)−1 X ⊤ y has the form (4.3) with fi = σi2 /(σi2 + λ) ∈ (0, 1). Proof. As above, X ⊤ X + λI = V (Σ2 + λI)V ⊤ and X ⊤ y = V ΣU ⊤ y, so xλ = V (Σ2 + λI)−1 ΣU ⊤ y =

X i

σi2 1 σi ⊤ (u y)v = · (u⊤ y)vi i i 2 2 σi + λ σi + λ σi i

which represents the solution in the form (4.3). Hence fi ≈ 1 when σi2 ≫ λ and fi ≈ 0 when σi2 ≪ λ: ridge regression is a smoothed truncated SVD, in which λ rather than k determines where the transition occurs (Figure 4.5, left). No singular value decomposition is needed, but a value of λ still has to be supplied. error

fi truncated SVD

1 ridge σi2

k⋆

k

Figure 4.5: Left: the truncated SVD cuts small-σ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 ⋆ .

4.3.4

Early stopping: the filter generated by CG

We now consider a third possibility, in which no filter is designed and no parameter selected: solve the unregularised normal equations in (4.1) by CG from x(0) = 0, and stop after k steps. The iteration then produces a filter of its own.

4.3. KRYLOV METHODS AND MACHINE LEARNING (OPTIONAL)

57

Proposition 4.3.3 (CG is a spectral filter). Apply CG to X ⊤ X x = X ⊤ y from x(0) = 0, and let qk be the CG residual polynomial of Proposition 4.1.8, so deg qk ≤ k and qk (0) = 1. Then x(k) has the form (4.3) with fi = 1 − qk (σi2 ), and among all filters built this way from a polynomial, that of CG minimises the A-norm of the error at every step. Proof. The matrix A = X ⊤ X is SPD with eigenpairs (σi2 , vi ) and Ax = X ⊤ y has solution xls . By Proposition 4.1.8, e(k) = qk (A)e(0) ; since x(0) = 0 we have e(0) = −xls , so x(k) = xls + e(k) = (I − qk (A))xls . Expanding xls in the basis {vi } and using qk (A)vi = qk (σi2 )vi gives the filter; the optimality is that of Proposition 4.1.8, read through this correspondence. Since qk (0) = 1 and qk must be small where the eigenvalues carry weight, the CG filter is close to 1 on the large σi2 while remaining close to 0 on the small ones after a few iterations. In other words, CG resolves the well-determined directions first. Stopping early is therefore a form of regularisation. The iteration count k plays the role that k plays for the truncated SVD and λ for ridge regression, but at no additional cost, since the iteration has to be carried out in any case. To be noted that the first two filters are fixed before the data are seen, whereas qk is built from Kk (A, X ⊤ y) and hence depends on y. CG selects a filter adapted to the problem at hand, which is why it outperforms a fixed filter at the same k. Remark 4.3.4 (When to stop). Filtering explains why early stopping helps, but not when to stop. Two effects compete, as shown on the right of Figure 4.5: • the part of the solution not yet resolved decreases with k; • the noise admitted through the filter grows with k, since the filter factors approach 1 on the small σi for which 1/σi is large. Their sum therefore attains a minimum at some k ⋆ and increases again beyond it, a phenomenon known as semi-convergence. A stopping rule is thus indispensable. If the noise level δ is known, one stops once ∥y − Xx(k) ∥2 ≤ τ δ with τ > 1, since fitting the data more closely than the noise amounts to fitting the noise; this makes CG a genuine regularisation method [17]. In machine learning δ is seldom known, and one stops when the error on a held-out validation set begins to rise. Takeaways • The least-squares solution divides by σi : the directions the data determine worst are exactly those noise corrupts most, so solving (4.1) exactly is the wrong thing to do. • Every remedy has the filtered form (4.3) and differs only in the shape of fi , each replacing the divergent 1/σi by something bounded. • Truncated SVD: fi = 1 for i ≤ k, else 0 is transparent, but it needs a full SVD and a cutoff k chosen by hand. • Ridge regression: fi = σi2 /(σi2 + λ), only matrix–vector products, one pass over the data each, so CG applies directly; but λ must be supplied. • CG, stopped early: fi = 1 − qk (σi2 ), no SVD and no parameter, optimal in the A-norm at every step, and adapted to the data through the Krylov space. • Early stopping is regularisation, with k in the role of λ. Semi-convergence means a stopping rule is still required.

58

4.4

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

Exercises

Problem 1 (A link between fixed point and Krylov iterative methods, exam 2023/2024). Let A, M ∈ Rn×n be non-singular, I := In , and b ∈ Rn . Let x0 ∈ Rn be arbitrary and define rk = b − Axk for k = 0, 1, 2, . . . . Abbreviate r := b − Ax0 . Assume one splits A = M − (M − A). (a) Show that Ax = b if and only if x = x + M −1 (b − Ax).

(4.4)

(b) Relation (4.4) leads to the fixed-point method xk+1 = xk + M −1 rk . Show that rk+1 = (I − AM −1 )rk . (c) Show that M −1 rk+1 = (I − M −1 A) M −1 rk . (d) Define the difference uk = xk+1 − xk . Show that uk = M −1 rk . (e) Define U k = [ u0 , u1 , . . . , uk−1 ] ∈ Rn×k . Show that   U k = M −1 r, (I − M −1 A)M −1 r, . . . , (I − M −1 A)k−1 M −1 r . Problem 2 (Krylov spaces, homework 2023/2024). Let n = 5 and 

 2 −1 0 0 0 −1 2 −1 0 0    A =  0 −1 2 −1 0  . 0 0 −1 2 −1 0 0 0 −1 2 One finds that its Schur decomposition is A = SDS −1 with  1 √

D = diag 2 −

√

√
3, 1, 2, 3, 2 + 3 ,

12

 1  2  S =  √13  1  2

√1 12

− 12 − 12 0 1 2 1 2

√1 3

1 2

0 − √13 0

− 21 0

√1 3

− 12

1 2

− √112



1  2  √1  . 3  1   2 − √112

Write S = [s1 , . . . , s5 ] as a matrix of five columns s1 , . . . , s5 . (a) Determine the minimal polynomial of A and the polynomial q for which A−1 = q(A). (b) Determine dim Kk (A, s2 ) and dim Kk (A, e1 ) for k = 1, 2, 3, . . . . (c) Let b = s2 + s4 and x0 = 0. Does there exist (in exact arithmetic) a k ∈ N for which the CG iterate xk equals the solution x̂ of Ax̂ = b? If yes, what is that value of k?

59

4.4. EXERCISES

Problem 3 (Kantorovich inequality). Let A be SPD with condition number κ(A). Show that for any v ̸= 0, !2 1 p 1 (v ⊤ Av)(v ⊤ A−1 v) ≤ κ(A) + p 1≤ . 4 (v ⊤ v)2 κ(A) Problem 4 (Orthogonality, A-conjugacy, best approximation). Let A be an SPD matrix. Prove the following fundamental CG properties: (a) (Residual orthogonality) r(i)⊤ r(j) = 0 for i ̸= j. (b) (A-orthogonality of search directions) p(i)⊤ Ap(j) = 0 for i ̸= j. (c) (Best A-norm approximation) x(k) minimizes ∥x − x∗ ∥A over x(0) + Kk (A, r(0) ). Problem 5 (Chebyshev bound and CG convergence rate). Let 0 < m = λmin (A) ≤ λ ≤ M = λmax (A) and κ(A) = M/m. (a) Show that ∥e(k) ∥A ≤ min p∈Πk ∥e(0) ∥A

max |p(λ)|, λ∈[m,M ]

p(0)=1

where Πk denotes polynomials of degree ≤ k. (b) We will assume as known the Chebyshev optimality property max |pk (λ)| = min

1

max |p(λ)| =

p∈Πk λ∈[m,M ] p(0)=1

λ∈[m,M ]

Tk



M +m M −m



where Tk is the Chebyshev polynomial having the complex representation √ √ (z + z 2 − 1)k + (z − z 2 − 1)k . Tk (z) = 2 Prove the bound !k p κ(A) − 1 , max |pk (λ)| ≤ 2 p λ∈[m,M ] κ(A) + 1 and conclude the classical CG estimate ∥e

(k)

∥A ≤ 2

!k p κ(A) − 1 p ∥e(0) ∥A . κ(A) + 1

Problem 6 (Orthogonality of Lanczos vectors). Let A = A⊤ and generate {qj } by the Lanczos steps w ← Aqj − βj−1 qj−1 , w ← w − αj qj ,

αj ← qj⊤ w,

βj ← ∥w∥2 ,

qj+1 ← w/βj

(βj ̸= 0),

with ∥q1 ∥2 = 1, q0 = 0, β0 = 0. Show that, in exact arithmetic, the vectors q1 , . . . , qm are orthonormal. Problem 7 (Networks: Lanczos Ritz values and interlacing). Let A ∈ Rn×n be symmetric, and let Tm = Q⊤ m AQm be the Lanczos tridiagonal matrix. Show that the eigenvalues of Tm interlace those of A, i.e. λi (A) ≤ θi ≤ λi+n−m (A), where θi are the eigenvalues of Tm .

60

CHAPTER 4. KRYLOV METHODS FOR SYMMETRIC SYSTEMS

Chapter 5

Krylov Methods for Non-Symmetric Systems Overview In Chapter 4 we exploited symmetry to the full but unfortunately none of these results will hold when A is non-symmetric. Convection–diffusion discretisations with upwinding, linearised Navier–Stokes Jacobians, transport and reaction networks, and the Google matrix are all nonsymmetric and often non-normal, meaning they have no orthogonal eigenbasis at all. This chapter develops the tools for that general setting. Replacing Lanczos by the Arnoldi process yields an upper Hessenberg projection H k and a full (long) recurrence; minimising ∥b − Ax∥2 over x(0) + Kk (A, r(0) ) then gives GMRES [34], at a storage and orthogonalisation cost that grows with k which is why restarting is needed in practice. The convergence theory changes character too: for non-normal A the eigenvalues alone no longer predict behaviour, and one turns to the field of values instead. This chapter applies the methodology to a wider class of problem: (1) Arnoldi and GMRES for Ax = b; (2) directed networks, where the same Arnoldi decomposition is used for eigenvalue computations; (3) rectangular systems from machine learning, where least squares poses the same minimisation but where the Krylov space is built from A and A⊤ . Learning objectives By the end of this chapter, you should be able to: • Derive the Arnoldi decomposition AVk = Vk+1 H k , explain why the projection is upper Hessenberg, and recover Lanczos as the symmetric special case. • Formulate GMRES as the least-squares problem miny ∥βe1 − H k y∥2 , and state its residual-minimisation, monotonicity, and finite-termination properties. • Express the GMRES residual as qk (A)r(0) with qk (0) = 1, and use eigenvalue and field-of-values bounds to explain why eigenvalues alone mispredict convergence when A is non-normal. • Formulate PageRank as an eigenproblem for the Google matrix, apply G without forming it, explain why |λ2 (G)| ≤ α bounds the power-iteration rate. • Explain why rectangular A forces a Krylov space built from A and A⊤ , and reduce LSQR to miny ∥γe1 − Bk y∥2 , mirroring the GMRES methodology.

61

62

5.1

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

GMRES and Arnoldi for non-normal matrices

We now wish to solve Ax = b when A is no longer symmetric. Two routes are possible in this case. CG requires A to be SPD, so it does not apply. Forming the normal equations A⊤ Ax = A⊤ b restores symmetry, but squares the condition number and destroys the sparsity of A; Section 5.3 will explain why this is a bad idea. What remains is the natural alternative: work directly with A, build a Krylov basis of Kk (A, r(0) ), and choose from that subspace the iterate whose residual is smallest in the 2-norm. Definition 5.1.1 (Krylov subspace and Arnoldi decomposition). Given A ∈ Rn×n (not necessarily symmetric) and r(0) = b − Ax(0) , Kk (A, r(0) ) = span{r(0) , Ar(0) , . . . , Ak−1 r(0) }. The Arnoldi process [2] (with v1 = r(0) /∥r(0) ∥2 ) builds an orthonormal basis Vk+1 = [v1 , . . . , vk+1 ] and an upper Hessenberg matrix H k ∈ R(k+1)×k such that   Hk AVk = Vk+1 H k , Hk = , hk+1,k e⊤ k where Hk ∈ Rk×k is upper Hessenberg with entries hij = (Avj , vi ) for 1 ≤ i ≤ j + 1 ≤ k. Hk

Vk+1

A Vk

nonzero k cols

=

one subdiagonal: upper Hessenberg

k + 1 cols 0

hk+1,k : the extra row

Figure 5.1: The Arnoldi decomposition AVk = Vk+1 H k , drawn for k = 4. Applying A to the orthonormal basis Vk of Kk never leaves Kk+1 , so the result factors through exactly one extra basis vector. Column j of H k is nonzero only in its first j + 1 entries, which is what upper Hessenberg means; in particular the first column has just two nonzeros. The matrix is (k + 1) × k, the last row containing the single entry hk+1,k that carries vk+1 . When A is symmetric, H k is additionally symmetric and therefore tridiagonal, and Arnoldi reduces to the Lanczos algorithm of Chapter 4. What is an upper Hessenberg matrix? A matrix H = [hij ] is upper Hessenberg if hij = 0 whenever i > j + 1, i.e. all entries strictly below the first subdiagonal are zero. Example (k = 6):   ⋆

⋆ 0 H= 0 0 0

⋆ ⋆ ⋆ 0 0 0

⋆ ⋆ ⋆ ⋆ 0 0

⋆ ⋆ ⋆ ⋆ ⋆ 0

⋆ ⋆ ⋆ . ⋆  ⋆ ⋆

Arnoldi produces such an Hk because Avj only has components along v1 , . . . , vj+1 .

5.1. GMRES AND ARNOLDI FOR NON-NORMAL MATRICES

63

Large non-symmetric systems were becoming routine in CFD and transport. Memory budgets were counted in megabytes, so one could not afford dense factorizations. Saad and Schultz’s 1986 idea [34] was elegantly pragmatic: let Arnoldi produce an orthonormal snapshot of the action of A, and within that small window solve the best possible least-squares.

5.1.1

GMRES: minimising the residual over the Krylov space

Arnoldi supplies a basis; it does not yet say which vector in the subspace to pick. GMRES makes the choice that CG made in Chapter 4, but in the 2-norm of the residual rather than the A-norm of the error, the latter no longer being available when A is not SPD. Definition 5.1.2 (GMRES — the principle (method definition)). Given A ∈ Rn×n , b, and an initial guess x(0) with residual r(0) = b − Ax(0) , the Generalized Minimal Residual method (GMRES) defines the k-th iterate as x(k) = arg

min x∈x(0) +Kk (A,r(0) )

Kk (A, r(0) ) = span{r(0) , Ar(0) , . . . , Ak−1 r(0) }.

∥b − Ax∥2 ,

Proposition 5.1.3 (Equivalence of the GMRES definition and the Arnoldi least-squares). Let r(0) = b − Ax(0) and Vk+1 = [v1 , . . . , vk+1 ] be an orthonormal Arnoldi basis of Kk (A, r(0) ) (0) with v1 = ∥rr(0) ∥ and set β = ∥r(0) ∥2 . Then the GMRES iterate x(k) is obtained by choosing 2

x(k) = x(0) + Vk y (k) , where y (k) solves the small least-squares y (k) = arg min

y∈Rk

βe1 − H k y

2

,

and the residual norm satisfies ∥r(k) ∥2 = ∥βe1 − H k y (k) ∥2 . Proof. Every x ∈ x(0) + Kk (A, r(0) ) can be written uniquely as x = x(0) + Vk y for some y ∈ Rk , because Vk has orthonormal columns spanning Kk . For such x, r(y) := b − Ax = b − A(x(0) + Vk y) = r(0) − AVk y. By Arnoldi, AVk = Vk+1 H k ; also r(0) = βv1 = Vk+1 (βe1 ). Hence
r(y) = Vk+1 βe1 − H k y . Because Vk+1 has orthonormal columns, it is an isometry on Rk+1 : ∥Vk+1 z∥2 = ∥z∥2 for all z ∈ Rk+1 . Therefore ∥r(y)∥2 = βe1 − H k y 2 . Minimizing ∥r(y)∥2 over x ∈ x(0) + Kk is thus exactly the (k+1) × k least-squares y (k) = arg min βe1 − H k y 2 , y

and the GMRES iterate is x(k) = x(0) + Vk y (k) with residual norm ∥r(k) ∥2 = ∥βe1 − H k y (k) ∥2 . If hk+1,k = 0, then H k has a zero last row, so the system is square (and upper Hessenberg) and the least-squares residual can be made zero; consequently r(k) = 0 and GMRES terminates in k steps. Remark 5.1.4 (Why orthonormality is the key feature). The proof used Vk+1 only once, but decisively: because its columns are orthonormal it preserves norms, so the n-dimensional residual and the (k + 1)-dimensional vector βe1 − H k y have exactly the same length. Had the basis merely ⊤ V been linearly independent, the reduced problem would have acquired a weighting by Vk+1 k+1 , and minimising it would no longer minimise the true residual. This is the same reason Lanczos orthonormalises in Chapter 4.

64

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

Proposition 5.1.3 turns each GMRES step into a small least-squares problem, but solving it from scratch at every step would be wasteful: H k differs from H k−1 by a single new column. Maintaining a QR factorisation of H k incrementally costs one rotation per step, and returns the residual norm as a by-product. Definition 5.1.5 (Givens rotations (notation used by GMRES)). For scalars a, b, a Givens rotation G(c, s) is the 2 × 2 orthogonal matrix   b a c s , s= √ , G(c, s) = , c= √ 2 2 2 −s c a +b a + b2   √  a a2 + b 2 7→ which maps . In GMRES we apply Gk to rows k and k+1 of the tall matrix b 0 H k to zero the subdiagonal entry hk+1,k , and we apply the same rotation to the right-hand side g = βe1 to keep g = Q⊤ (βe1 ) up to step k. Definition 5.1.6 (Full GMRES (Givens form) — construction). Given x(0) , set r(0) = b − Ax(0) , β = ∥r(0) ∥2 , v1 = r(0) /β. For k = 1, 2, . . . until convergence: 1. Arnoldi expansion. Compute w = Avk . For i = 1, . . . , k: hi,k = vi⊤ w,

w ← w − hi,k vi .

Set hk+1,k = ∥w∥2 . If hk+1,k = 0 (lucky breakdown), stop; else vk+1 = w/hk+1,k . 2. Update small LS (Givens). Apply the stored rotations G1 , . . . , Gk−1 to the new column of H k , then form Gk to annihilate hk+1,k and apply Gk to g = βe1 ; this maintains the QR of H k and yields ∥r(k) ∥2 = |gk+1 |. 3. Residual and test. The current residual norm is ∥r(k) ∥2 = |gk+1 |, available without forming x(k) . If ∥r(k) ∥2 /∥b∥2 ≤ tol, recover y (k) by back–substitution in the k × k upper triangular R and set x(k) = x(0) + Vk y (k) . Proposition 5.1.7 (GMRES basics in one shot). Let A ∈ Rn×n , x(0) be an initial guess, r(0) = b − Ax(0) , and let Vk+1 and H k come from k steps of Arnoldi with v1 = r(0) /∥r(0) ∥2 and AVk = Vk+1 H k ; set β = ∥r(0) ∥2 . The GMRES iterate x(k) = x(0) + Vk y (k) , where y (k) = arg miny ∥βe1 − H k y∥2 . This satisfies: (i) Orthogonality. r(k) ⊥ AKk (A, r(0) ). (ii) Monotonicity (full GMRES). ∥r(k+1) ∥2 ≤ ∥r(k) ∥2 . (iii) Residual polynomial and finite termination. There exists qk ∈ Πk with qk (0) = 1 such that r(k) = qk (A) r(0) , and qk minimizes ∥q(A)r(0) ∥2 over all such polynomials. If m is the minimal polynomial of A relative to r(0) with deg m = d, then GMRES terminates in at most d steps (i.e., r(d) = 0).

Proof. See Problem 3.

65

5.1. GMRES AND ARNOLDI FOR NON-NORMAL MATRICES normal A

non-normal A

ℑ

ℑ

F (A)

ℜ

0

ℜ

same eigenvalues, but 0 ∈ F (A) ⇒ GMRES may stagnate

eigenvalues clustered away from 0 ⇒ GMRES converges fast

Figure 5.2: Why eigenvalues are not enough. For a normal matrix the eigenvalues alone control ∥q(A)∥2 = maxj |q(λ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}, a convex set containing σ(A). If 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), no such bound is available, and GMRES can stagnate for n − 1 steps even with all eigenvalues equal to 1 (see the example at the end of this section).

5.1.2

Convergence: why eigenvalues are not enough

For CG, Chapter 4 gave a single convergence bound in terms of κ(A), and that was the end of the story. For GMRES there is no such clean result. The natural first attempt is to bound the residual using the eigenvalues of A, and it works whenever A is diagonalisable with a well-conditioned eigenvector matrix, but it becomes vacuous exactly when A is strongly non-normal, which is the case we introduced the method for. Proposition 5.1.8 (Diagonalizable-A bound). If A = V ΛV −1 , then ∥r(k) ∥2 ≤ κ2 (V ) min max |q(λi )|. i q∈Πk ∥r(0) ∥2 q(0)=1

Proof. r(k) = qk (A)r(0) = V qk (Λ)V −1 r(0) , so ∥r(k) ∥ ≤ ∥V ∥ ∥V −1 ∥ maxi |qk (λi )| ∥r(0) ∥. Definition 5.1.9 (Field of values (numerical range)).  ∗  x Ax n F (A) = : x ∈ C \ {0} . x∗ x In words, F (A) is the set of all Rayleigh quotients of A. It is compact and convex, and it contains every eigenvalue of A, as one sees by taking x to be an eigenvector. The point of replacing the spectrum by F (A) is that the latter is insensitive to the conditioning of the eigenvectors: its shape and position control ∥q(A)∥2 directly, and hence the GMRES residual through r(k) = qk (A)r(0) . Proposition 5.1.10 (Disk-in-F (A) bound). If F (A) ⊂ { |z − c| ≤ s } and 0 ∈ / { |z − c| ≤ s }, then √  s k (0) ∥r(k) ∥2 ≤ (1 + 2) ∥r ∥2 . |c|

66

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

√ Proof. Crouzeix-Palencia: ∥q(A)∥ ≤ (1 + 2) maxz∈F (A) |q(z)|. Test the degree-k polynomial qk (z) = (1 − z/c)k ; then maxF (A) |qk | ≤ (s/|c|)k . When eigenvalues mislead: a companion matrix. The following example shows how far the spectrum alone can be from the truth. For a monic polynomial p(z) = z n + cn−1 z n−1 + · · · + c0 , the companion matrix   0 1 0 ··· 0  0 0 1 ··· 0     .. .. ..  .. .. C(p) =  . . . . .     0 0 ··· 0 1  −c0 −c1 · · · −cn−2 −cn−1 has eigenvalues equal to the roots of p. These matrices are typically highly nonnormal. Take for example p(z) = (z − 1)3 = z 3 − 3z 2 + 3z − 1. Then   0 1 0   A = C(p) = 0 0 1 , Λ(A) = {1, 1, 1}. 1 −3 3 Run GMRES on Ax = f with x(0) = 0 and f = e1 . One observes ∥r(0) ∥2 = ∥r(1) ∥2 = ∥r(2) ∥2 = 1,

∥r(3) ∥2 = 0,

i.e., no decrease for two steps, then exact termination at k = 3. Moreover, e⊤ 1 Ae1 = 0, so 0 ∈ F (A) and the disk bound above does not apply here. This illustrates that tightly clustered eigenvalues (even all at 1) do not guarantee fast GMRES if the matrix is nonnormal; the field of values and nonnormality matter. In fact any nonincreasing residual curve whatsoever can be realised by some matrix with any prescribed spectrum [14]. A systematic account of non-normal behaviour is given in [36].

5.2

Networks: PageRank as an Eigenproblem

The original PageRank algorithm, developed by Larry Page and Sergey Brin in the late 1990s, was essentially an eigenvector computation on the web graph. At the time, the web already had tens of millions of pages, making full eigen-decompositions infeasible. Instead, iterative methods such as power iteration were the only scalable option. This numerical perspective i.e. treating web search as a giant linear algebra problem, was one of the key insights that allowed Google to leap ahead of earlier search engines. It remains a striking example of how classical numerical methods became transformative in an entirely new application domain. Chapter 4 used the Krylov subspace twice: to solve Ax = b and, through Lanczos, to compute eigenvalues of a symmetric matrix. Section 5.1 has just repeated the first of these for non-symmetric A; this section repeats the second. The Arnoldi decomposition AVm = Vm+1 H m is again all we need for the eigenvalues of the small matrix Hm rather than for a linear solve. The application we envision here is ranking the nodes of a directed graph, where the matrix is non-symmetric by construction and far too large to factorise. Definition 5.2.1 (Eigenvector centrality (directed graphs)). Let A ∈ Rn×n be the (possibly non-symmetric) adjacency matrix. Right eigenvector centrality assigns scores x ≥ 0 to nodes by solving Ax = λ1 x, ∥x∥1 = 1, while left eigenvector centrality assigns scores y ≥ 0 from A⊤ y = λ1 y,

∥y∥1 = 1.

67

5.2. NETWORKS: PAGERANK AS AN EIGENPROBLEM The Google matrix Let P ∈ Rn×n be the column-stochastic matrix where   1 , if page j links to page i, Pij = dj  0, otherwise, with dj the out-degree of page j. The Google matrix is defined as

1 G = αP + (1 − α) 11⊤ , 0 < α < 1. n By construction, G is stochastic1 and irreducible2 ; the Perron-Frobenius theorem guarantees a unique dominant eigenvector x with positive entries such that ∥x∥1 = 1, called the PageRank vector [24, 29]. Efficient application of the Google matrix We should notice that G is dense: every entry of 11⊤ is nonzero, so storing G for a web-scale graph is out of the question, but it doesn’t have to be stored. For any x with 1⊤ x = 1 the rank-one perturbation term acts as a constant vector, 1−α 1−α Gx = αP x + 1 (1⊤ x) = αP x + 1, n n so one application of G costs one sparse product with P plus a vector addition. The power iteration of Chapter 3 (Algorithm 1) applied to the dominant eigenvector is therefore simply 1−α 1, (5.1) x(k+1) = αP x(k) + n one sparse matrix–vector product per step, two vectors of storage, and no factorisation anywhere. This is the entire algorithm that ranked the early web. Example 5.2.2 (Google Page Rank). Consider a graph of three webpages with links and find which is the most central node 1

2

3 The adjacency matrix is

  0 0 1 A = 1 0 1 , 0 1 0 With damping factor α = 0.85,

  0 0 1/2 P = 1 0 1/2 . 0 1 0

  0.05 0.05 0.475 1−α ⊤ G = αP + 11 = 0.85 P + 0.05 11⊤ = 0.90 0.05 0.475 , n 0.05 0.90 0.05

(column-stochastic).

Start from the uniform vector x(0) = 31 (1, 1, 1)⊤ . Iterating with the power method gives   0.215 x⋆ ≈ 0.397 . 0.388 hence node 2 is the most central. 1

Every entry is a value from 0 to 1, and the numbers in each column (or row) add up to exactly 1 An irreducible matrix is a square matrix that cannot be rearranged into a block upper-triangular form using row and column permutations. It means you cannot split the matrix into smaller independent parts. 2

68

5.2.1

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

The role of the damping factor

First, what α means. The Google matrix models a random surfer who, at each step, follows one of the links on the current page with probability α, and with probability 1 − α abandons the page and jumps to a page chosen randomly. So α is a modelling choice. It is also the only thing that determines the cost of computing the ranking. Proposition 5.2.3 (Role of the damping factor [20]). Let P be column-stochastic, 0 < α < 1, ⊤ and G = αP + 1−α n 11 . Then every eigenvalue of G other than λ1 = 1 satisfies |λ| ≤ α. Consequently the power iteration (5.1) converges linearly with asymptotic factor at most α, independently of the size of the graph. Proof. Because G is column-stochastic, 1⊤ G = 1⊤ : the vector 1 is a left eigenvector for the eigenvalue 1. Let Gx = λx with λ ̸= 1. Left and right eigenvectors belonging to different eigenvalues are orthogonal, so 1⊤ x = 0. Then we have, Gx = αP x +

1−α 1 (1⊤ x) = αP x, n

so x is also an eigenvector of P , with eigenvalue λ/α. Since P is stochastic all its eigenvalues satisfy |µ| ≤ 1, hence |λ| ≤ α. Two consequences are resulting from this but they pull in opposite directions. • The method is usable at scale. The bound involves α and nothing else, in particular not n. Doubling the size of the web does not cost a single extra iteration, which is why (5.1) scaled from millions of pages to billions unchanged. • The model and the cost are linked. Taking α close to 1 makes the ranking faithful to the link structure and simultaneously makes the iteration slow: to reach a tolerance ε needs about log ε/ log α steps, which behaves like log(1/ε)/(1 − α) as α → 1 and therefore blows up. The three-node example makes this concrete. Its link matrix P has characteristic polynomial λ3 − 12 λ − 12 = (λ − 1)(λ2 + λ + 12 ), so besides λ = 1 its eigenvalues are (−1 ± i)/2, of modulus √ 1/ 2. By the√proof above the corresponding eigenvalues of G are exactly α times these, giving |λ2 (G)| = α/ 2 ≈ 0.707 α, which is what Table 5.1 reports. The bound |λ2 | ≤ α is thus attained here up to the factor |λ2 (P )|; on a real web graph |λ2 (P )| is very close to 1 and the bound is essentially sharp. α

|λ2 (G)|

power iterations to 10−8

0.50 0.85 0.95 0.99

0.354 0.601 0.672 0.700

17 36 45 52

√ Table 5.1: The damping factor on the three-node example of this section, where |λ2 (G)| = α/ 2 exactly. The last column is the number of iterations observed to gain eight digits, in agreement with the asymptotic estimate 8 log 10/ log(1/|λ2 (G)|). The value α = 0.85 used by Google is a compromise between fidelity to the link structure and cost, not a mathematical constant.

5.2.2

Beyond the power iteration: the Arnoldi process

Power iteration is attractive because (5.1) is one line of code, but it delivers a single eigenvector and its rate is fixed at α: it cannot be accelerated, and it has no way of certifying its answer.

69

5.3. KRYLOV METHODS FOR LARGE LEAST SQUARES (OPTIONAL)

Remark 5.2.4 (When is power iteration not enough?). Four situations call for Arnoldi instead. • Small spectral gap: if |λ1 | ≈ |λ2 |, power converges very slowly. • Multiple/clustered targets or several eigenpairs: power returns at most one direction; Arnoldi delivers many Ritz pairs at once. • Nonnormal matrices: eigenvectors may be poorly conditioned and power can behave erratically; Arnoldi tracks subspaces and provides residuals ∥Au − θu∥. • Stopping guarantees: Arnoldi gives rigorous a posteriori residuals |hm+1,m ||e⊤ m y|. The remedy is the Arnoldi process, which we run on A (here A = G, or a shifted problem), and from the decomposition AVm = Vm+1 H m take an eigenpair (θ, y) of the square part Hm . The Ritz pair is then (θ, u) with u = Vm y, and the associated residual is available for free, ∥Au − θu∥2 = |hm+1,m | |e⊤ m y|,

(5.2)

because Au − θu = (AVm − Vm Hm )y = hm+1,m vm+1 (e⊤ m y) and ∥vm+1 ∥2 = 1. Takeaways • PageRank is the dominant eigenvector of the Google matrix G; Perron–Frobenius guarantees it is unique and positive. • G is dense but never formed: one application costs one sparse product with P , so the power iteration (5.1) is one line of code. • The damping factor is the single difficult parameter; |λ2 (G)| ≤ α bounds the convergence rate independently of the size of the graph. • Power iteration gives one eigenpair at a fixed rate and no error guarantees; in turn, Arnoldi i.e. the same decomposition built in Section 5.1 gives several Ritz pairs at once, plus the residual |hm+1,m ||e⊤ m y| for free.

5.3

Krylov Methods for Large Least Squares (optional)

Section 5.1 solved Ax = b for square non-symmetric A by minimising ∥b − Ax∥2 over a Krylov subspace. Least squares asks for the minimiser of the same quantity when A is rectangular, and this section is about what changes. The answer is a single structural point: with A rectangular the sequence r(0) , Ar(0) , A2 r(0) , . . . does not even make sense, because Ar(0) lives in the wrong space. The Krylov space must be built by alternating A and A⊤ , and everything else follows from this single change. 3 A motivating example: sparse regression on text Consider a linear model fitted to a corpus of documents: row i of the design matrix X holds the word counts of document i over a vocabulary of n terms, and bi is a response, as for example a rating in recommender systems, a click-through rate, or a price. Fitting by least squares means solving min ∥Xx − b∥2 , x

3

X ∈ Rm×n ,

This section is not examinable but the ideas in it are likely to be useful for the project.

70

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

with m in the millions and n in the tens of thousands. Three features of this problem decide the choice of algorithm, and all of them are frequently encountered (recommender systems, genomics, and inverse problems). • X is rectangular, so “solve the linear system” is not even the right statement; there is no square matrix in sight and generally no exact solution. • X is extremely sparse. A document contains a few hundred of the n vocabulary terms, so X has a fraction of a percent of its entries nonzero. • X is ill-conditioned. Word frequencies are heavy-tailed and features are correlated (“machine” and “learning” co-occur), so the singular values of X decay over many orders of magnitude. In what follows will briefly explain why the obvious route fails. b r R(A) Ax∗

r ⊥ R(A) ⇐⇒ A⊤ (Ax∗ − b) = 0

Figure 5.3: Least squares as orthogonal projection. The minimiser of ∥Ax − b∥2 is the point of the column space R(A) closest to b, so the residual is orthogonal to R(A), which leads exactly to the normal equations A⊤ Ax = A⊤ b. LSQR realises the same projection without forming A⊤ A. As explained in Chapter 4, the first ideas would be to solve minx ∥Xx − b∥2 by differentiating: the minimiser satisfies the normal equations X ⊤ Xx = X ⊤ b, a square SPD system of size n, to which the CG of Chapter 4 applies immediately. Indeed Section 4.3 did precisely this: it read least squares as an SPD problem and ran CG on X ⊤ X in order to exhibit the filter factors 1 − qk (σi2 ). That was the right way to understand regularisation, because the filtering is visible only in the eigenbasis of X ⊤ X. It is not the right way to compute. At any serious scale one never forms X ⊤ X, and never runs CG on it, for two independent reasons. Sparsity is destroyed. Two terms contribute a nonzero to (X ⊤ X)jk as soon as they co-occur in a single document, so the Gram matrix is far denser than X. On a representative sparse design matrix, for example 20 000 × 5000 at 0.2% density, forming X ⊤ X produces a matrix that is 7.7% dense, nearly ten times as many nonzeros as X itself, and this ratio grows with the corpus. Conditioning is squared. Since the singular values of X ⊤ X are σi2 , κ2 (X ⊤ X) = κ2 (X)2 . ⊤ X (the method known as CGNR) therefore converges at a rate governed by Running CG on Xp κ2 (X) rather than κ2 (X), and in floating point the damage can be worse than slow: information present in X can be destroyed by the time X ⊤ X has been formed.

71

5.3. KRYLOV METHODS FOR LARGE LEAST SQUARES (OPTIONAL) A classical warning. Let ε > 0 and     1 1 1 + ε2 1 ⊤   X= ε 0 , X X= . 1 1 + ε2 0 ε

√ X has rank 2 for every ε = ̸ 0. But if ε < u, where u ≈ 2.2 · 10−16 is the unit roundoff, then 1 + ε2 rounds to 1 and the Gram matrix is [ 11 11 ], so exactly singular. In double precision this already happens at ε = 10−9 : the rank information is annihilated by forming X ⊤ X. What we want, then, is a method with the convergence of CG on X ⊤ X but the conditioning of X, that means one that uses X and X ⊤ only as operators and never multiplies them together. That is exactly what Golub–Kahan bidiagonalisation provides [11]. From here on we revert to the neutral notation A ∈ Rm×n for the rectangular matrix, to keep the parallel with Section 5.1 visible.

5.3.1

Golub–Kahan bidiagonalisation and LSQR

b Definition 5.3.1 (Golub–Kahan bidiagonalization (core recurrence)). Choose u1 = ∥b∥ and 2 set β0 = 0, v0 = 0. For k = 1, 2, . . . :

(i) r = A⊤ uk − βk−1 vk−1 , (ii) p = Avk − αk uk ,

αk = ∥r∥2 ,

βk = ∥p∥2 ,

vk = r/αk ,

uk+1 = p/βk .

Since Avk = αk uk + βk uk+1 , we obtain AVk = Uk+1 Bk with Bk ∈ R(k+1)×k lower bidiagonal, carrying α1 , . . . , αk on the diagonal and β1 , . . . , βk on the subdiagonal.

{vj } orthonormal in Rn (row space)

v1

v2 A⊤

A

u1

v3 A⊤

A

···

u2

{uj } orthonormal in Rm (column space)

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 A and one with A⊤ . The coefficients assemble into a lower bidiagonal Bk with AVk = Uk+1 Bk . This is Lanczos applied implicitly to A⊤ A but the squaring never happens numerically, which is precisely why LSQR is stable where CG on the normal equations is not. The reason this construction is worth the trouble is that it reduces the least-squares problem to a tiny one of exactly the shape met in Section 5.1.

72

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

Proposition 5.3.2 (LSQR as a small least-squares problem). Write γ = ∥b∥2 , so that u1 = b/γ, and let k steps of Golub–Kahan produce Vk , Uk+1 and the bidiagonal Bk ∈ R(k+1)×k with AVk = Uk+1 Bk . Then for any x = Vk y, ∥b − Ax∥2 = γe1 − Bk y 2 , and consequently the minimiser of ∥b − Ax∥2 over x ∈ R(Vk ) is x(k) = Vk y (k) with y (k) = arg min γe1 − Bk y 2 . y∈Rk

Proof. Since b = γu1 = Uk+1 (γe1 ) and AVk y = Uk+1 Bk y,
b − AVk y = Uk+1 γe1 − Bk y . The columns of Uk+1 are orthonormal, so Uk+1 is an isometry on Rk+1 and the norm is unchanged. Minimising over y gives the second claim. Compare this with Proposition 5.1.3: GMRES minimises ∥βe1 − H k y∥2 against an upper Hessenberg H k , LSQR minimises ∥γe1 − Bk y∥2 against a bidiagonal Bk . The algorithms differ only in how the small matrix is generated and, consequently, in how cheaply it is factorised: a bidiagonal Bk is reduced to triangular form by a single Givens rotation per step. This is what allows the LSQR iterate to be updated by short recurrences in O(m + n) work and storage, with no need to keep Vk . Mathematically it produces the same iterates as CG on A⊤ A, but the squaring never happens numerically, which is why it is stable where CGNR is not.

5.3.2

Early stopping and the truncated SVD σi

singular value index i captured in the first LSQR steps

reached late; these are the directions that amplify noise

Figure 5.5: LSQR as an iterative truncated SVD. Because the Krylov space is built from A⊤ b and grows by powers of A⊤ A, components associated with large singular values are resolved first. Stopping early therefore acts as a spectral filter, suppressing the small-σ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. Running ℓ = k +p steps of Golub–Kahan on a sparse term–document matrix M and computing the (tiny) SVD Bℓ = Û ΣV̂ ⊤ gives approximate singular triplets of M through U ≈ Uℓ+1 Û:,1:k and V ≈ Vℓ V̂:,1:k , with a small oversampling p to improve accuracy of the last few. This is latent semantic analysis: the leading singular directions of M define a low-dimensional semantic space in which related documents cluster, and it is obtained from a few dozen matrix–vector products rather than a dense O(mn min(m, n)) factorisation.

5.3. KRYLOV METHODS FOR LARGE LEAST SQUARES (OPTIONAL)

73

Vk doc-embedding (2D sketch) sports

finance

Figure 5.6: After truncated SVD, documents cluster in a low-dimensional semantic space. Takeaways • Rectangular A forces the Krylov space to be built from A and A⊤ ; Golub–Kahan bidiagonalisation is the resulting recurrence. • LSQR reduces the least-squares problem to miny ∥γe1 − Bk y∥2 , structurally identical to the GMRES problem of Section 5.1, but with a bidiagonal in place of a Hessenberg matrix. • Never form A⊤ A: the geometry of the normal equations is right, the arithmetic is not. • The same bidiagonalisation, stopped early, acts as a spectral filter (regularisation); kept and factorised, it yields a truncated SVD.

74

5.4

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

Exercises

Problem 1 (Krylov spaces, resit exam 2024/2025). Consider the matrix 

 2 −1 0 A = −1 2 −1 , 0 −1 2

  1 b = 0 . 1

(a) Define the Krylov space Kk (A, b) for k = 1, 2, . . . . (b) Compute the Krylov spaces K1 (A, b), K2 (A, b), and K3 (A, b). Determine their dimensions. (c) Explain the general principle of the GMRES method for solving Ax = b. How many iterations are needed for convergence? Problem 2 (Companion matrices and GMRES convergence, exam October 2024). Let   0 1 0   0 1     .. .. A= . . .    0 1  c0 c1 · · · cn−2 cn−1 P j and the polynomial p(z) = z n − n−1 j=0 cj z . Show that: (a) p(z) is its characteristic and minimal polynomial. (b) Let   0 1 0   A = 0 0 1 , 1 −3 3

p(z) = (z − 1)3 = z 3 − 3z 2 + 3z − 1,

  1 b = e1 = 0 . 0

Compute the Krylov spaces Kk (A, b), k = 1, 2, 3. (c) Compute the GMRES iterates uk for Au = b with u0 = 0 and conclude on how many iterations are needed for convergence. (d) Can this result be generalised to an n × n companion matrix? Explain your answer. Problem 3 (GMRES basics in one shot). Let A ∈ Rn×n , x(0) be an initial guess, r(0) = b − Ax(0) , and let Vk+1 and H k come from k steps of Arnoldi with v1 = r(0) /∥r(0) ∥2 and AVk = Vk+1 H k ; set β = ∥r(0) ∥2 . The GMRES iterate x(k) = x(0) + Vk y (k) , where y (k) = arg miny ∥βe1 − H k y∥2 . Prove that it satisfies: (i) Orthogonality. r(k) ⊥ AKk (A, r(0) ). (ii) Monotonicity (full GMRES). ∥r(k+1) ∥2 ≤ ∥r(k) ∥2 . (iii) Residual polynomial and finite termination. There exists qk ∈ Πk with qk (0) = 1 such that r(k) = qk (A) r(0) , and qk minimizes ∥q(A)r(0) ∥2 over all such polynomials. If m is the minimal polynomial of A relative to r(0) with deg m = d, then GMRES terminates in at most d steps (i.e., r(d) = 0).

75

5.4. EXERCISES

Problem 4 (Field of values and GMRES convergence). Let A ∈ Cn×n and suppose its field of values (numerical range) F (A) is contained in a disk {z : |z − c| ≤ s} with 0 ∈ / {z : (k) |z − c| ≤ s}. Let r denote the k-th GMRES residual for Ax = b with some initial guess x(0) and r(0) = b − Ax(0) . (a) Show that GMRES residuals admit a polynomial representation r(k) = qk (A)r(0) for some qk ∈ Πk with qk (0) = 1, and that GMRES chooses qk to minimize ∥q(A)r(0) ∥2 over all such polynomials. (b) Prove the “disk-in-F (A)” bound ∥r

(k)

 ∥2 ≤ C

s |c|

k

∥r(0) ∥2 ,

√ for a universal constant C. (Hint: Use Crouzeix–Palencia: ∥p(A)∥ ≤ (1+ 2) maxz∈F (A) |p(z)| and test p(z) = (1 − z/c)k .) (c) Explain briefly why this bound can be more descriptive than eigenvalue-only bounds for highly non-normal matrices. Give a concrete example of a matrix with tightly clustered eigenvalues but possibly slow GMRES and relate this to the geometry of F (A). Problem 5 (Arnoldi, Ritz values, and a posteriori eigen-residuals). Let A ∈ Cn×n and (0) suppose k steps of Arnoldi with v1 = ∥rr(0) ∥ produce 2

AVk = Vk+1 Hk ,

 Vk+1 = v1 · · ·

 vk+1 ,

 Hk =

bk H

hk+1,k e⊤ k

 ,

b k ∈ Ck×k is its where Vk+1 has orthonormal columns, Hk ∈ C(k+1)×k is upper Hessenberg, and H leading square part. b k with ∥y∥2 = 1, then (θ, v) with v := Vk y is a Ritz (a) Show that if (θ, y) is an eigenpair of H pair for A (i.e. θ is a Ritz value and v a Ritz vector). (b) Prove the a posteriori residual formula ∥Av − θv∥2 = |hk+1,k | |e⊤ k y|. (c) Interpret this bound: under what condition is θ a good eigenvalue approximation to A? How does this relate to (i) small subdiagonal element |hk+1,k | (a “lucky” or near-breakdown), b and (ii) the last component |e⊤ k y| of the eigenvector of Hk ? Problem 6 (PageRank on a four-page web). Four pages link as follows: page 1 → 2, 3; page 2 → 3; page 3 → 1; page 4 → 1, 3. (a) Write down the link matrix P (with Pij = 1/dj when page j links to page i) and check that it is column-stochastic. Page 4 receives no links: what does that say about row 4 of P , and what will make its PageRank nonzero nevertheless? (b) Take α = 0.85. Without writing out the dense matrix G, give the formula that computes Gx from P x, and say how many nonzeros are touched per application. (c) Starting from x(0) = 14 (1, 1, 1, 1)⊤ , perform two steps of the iteration (5.1) and rank the four pages.

76

CHAPTER 5. KRYLOV METHODS FOR NON-SYMMETRIC SYSTEMS

(d) Using Proposition 5.2.3, bound the number of iterations needed to reduce the error by 10−6 , first for α = 0.85 and then for α = 0.99. Comment on what the ranking gains and what the computation loses when α is increased. (e) Suppose a fifth page is added that links to nothing. Show that P is then no longer column-stochastic, and describe a way of repairing it that keeps the cost of applying G unchanged. Problem 7 (Golub–Kahan bidiagonalisation and LSQR). Let A ∈ Rm×n with m ≥ n, let b ̸= 0, and run the recurrence of Section 5.3.1 with γ = ∥b∥2 and u1 = b/γ. (a) Show that Avk = αk uk + βk uk+1 and A⊤ uk+1 = βk vk + αk+1 vk+1 . Deduce that AVk = Uk+1 Bk with Bk lower bidiagonal, and write B3 out explicitly. (b) Show by induction that span{v1 , . . . , vk } = Kk (A⊤ A, A⊤ b) and span{u1 , . . . , uk } = Kk (AA⊤ , b). This is the precise sense in which LSQR “is” CG on the normal equations. (c) Proposition 5.3.2 reduces step k to miny ∥γe1 − Bk y∥2 . Show that a single Givens rotation brings Bk to upper triangular form given the factorisation of Bk−1 , and contrast the work and the storage per step with GMRES, where step k costs k orthogonalisations against v1 , . . . , v k .

Chapter 6

Preconditioning and Accelerating Solvers Overview Chapters 4 and 5 showed that Krylov methods are optimal over the p subspace they build, but that their speed is dictated by the spectrum of A: CG needs O( κ(A)) iterations, and for the PDE matrices of Chapter 2 this means O(h−1 ) iterations, growing without bound as the mesh is refined. Optimality over a Krylov subspace cannot fix a badly conditioned problem. This is what preconditioning does. Instead of solving Ax = b we solve an equivalent system involving M −1 A, where M approximates A but is cheap to invert, and we choose M so that the transformed spectrum is clustered and the condition number is bounded independently of the mesh size; broad surveys of the subject are given in [4, 39]. This chapter covers: (1) the algebraic framework of left, right, and symmetric preconditioning, preconditioned CG, and simple choices such as Jacobi and incomplete factorisations; (2) Schwarz domain decomposition, where the preconditioner is built from independent subdomain solves and the overlap width controls the convergence rate; (3) two-level and multigrid ideas, where a coarse-grid correction removes the smooth error components that relaxation cannot touch, yielding preconditioners whose quality is independent of h. Learning objectives By the end of this chapter, you should be able to: • Explain why Krylov convergence rates make preconditioning necessary for mesh-refined PDE problems, and state the left, right, and symmetric preconditioning formulations. • Derive the preconditioned conjugate gradient algorithm and explain why it requires only solves with M , never M −1/2 , and why M must be SPD. • Construct and assess simple preconditioners (Jacobi, SSOR, incomplete Cholesky/LU) and estimate their effect on κ(M −1 A) for the discrete Laplacian. • Formulate the alternating and additive Schwarz methods, analyse the 1D convergence factor in terms of the overlap δ, and explain the role of a coarse space in restoring scalability. • Compute the smoothing factor of weighted Jacobi, explain why relaxation damps oscillatory but not smooth error, and describe how coarse-grid correction and the resulting two-grid/multigrid cycle give an h-independent preconditioner.

77

78

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS

6.1

Motivation and Basic Preconditioners

Preconditioning accelerates Krylov iterations by modifying the spectrum of the coefficient matrix. We aim to accelerate the solution of Ax = b, A ∈ Rn×n , using a matrix M that approximates A but is cheaper to invert. Definition 6.1.1 (Left, right, and symmetric preconditioning). We have three different options: • Left preconditioning:

ML−1 Ax = ML−1 b.

This changes the residual norm seen by the Krylov method. • Right preconditioning: AMR−1 y = b,

x = MR−1 y.

This preserves the original residual b − Ax but modifies the search space. • Symmetric preconditioning: for SPD problems, one uses an SPD matrix M and rewrites the system as M −1/2 AM −1/2 x̃ = M −1/2 b, x̃ = M 1/2 x. This preserves symmetry and positive definiteness. Key idea For SPD problems, the convergence of CG depends on the condition number of the operator. A good preconditioner replaces A by an operator whose eigenvalues are more tightly clustered, ideally near 1: M −1 A ≈ I,

κ(M −1 A) ≪ κ(A).

The gain can already be quantified exactly on simple model matrices.

M −1 A (clustered)

Unpreconditioned A 0.5 ℑ

ℑ 1

ℜ

ℜ 1

2

3

0.5

4

1

1.5

−1 −0.5

Figure 6.1: Effect of preconditioning on the spectrum: eigenvalues of A (left) are spread, whereas those of M −1 A (right) are tightly clustered around 1.

6.1.1

Preconditioning SPD problems

For SPD problems, preconditioning is most naturally understood through the symmetrically preconditioned operator e = M −1/2 AM −1/2 . A

79

6.1. MOTIVATION AND BASIC PRECONDITIONERS CG is then applied to

eb = M −1/2 b.

x e = M 1/2 x,

ex = eb, Ae

Equivalently, in the original variables, the search space is the preconditioned Krylov space x(0) + Kk (M −1 A, M −1 r(0) ). Proposition 6.1.2 (Preconditioned CG). Let A and M be SPD. Applying CG to M −1/2 AM −1/2 x e = M −1/2 b is equivalent to the following iteration in the original variables: z (0) = M −1 r(0) ,

r(0) = b − Ax(0) , q (k) = Ap(k) ,

αk =

⟨r(k) , z (k) ⟩ , ⟨p(k) , q (k) ⟩

x(k+1) = x(k) + αk p(k) , z (k+1) = M −1 r(k+1) ,

p(0) = z (0) ,

r(k+1) = r(k) − αk q (k) , βk+1 =

⟨r(k+1) , z (k+1) ⟩ , ⟨r(k) , z (k) ⟩

p(k+1) = z (k+1) + βk+1 p(k) . Moreover, spec(M −1/2 AM −1/2 ) = spec(M −1 A), so the usual CG convergence estimate holds with κ(M −1 A) replacing κ(A). Proof sketch. Set

x e = M 1/2 x,

re = M −1/2 r,

pe = M 1/2 p.

e = M −1/2 AM −1/2 gives Then standard CG on A αk =

⟨r(k) , M −1 r(k) ⟩ ⟨e r(k) , re(k) ⟩ = , ep(k) ⟩ ⟨p(k) , Ap(k) ⟩ ⟨e p(k) , Ae

which yields the stated recursion with z (k) = M −1 r(k) . Similarly, βk+1 = Finally,

⟨e r(k+1) , re(k+1) ⟩ ⟨r(k+1) , z (k+1) ⟩ = . ⟨e r(k) , re(k) ⟩ ⟨r(k) , z (k) ⟩

M −1 A = M −1/2 (M −1/2 AM 1/2 )M −1/2

is similar to M −1/2 AM −1/2 , so they have the same eigenvalues. Remark 6.1.3. Preconditioning does not alter the exact solution x∗ but changes the trajectory of residual minimization. The effectiveness depends on how well M approximates A in the subspace explored by the Krylov method.

6.1.2

A general framework for preconditioners

Jacobi preconditioning. Let A = D + L + U, where D is the diagonal of A. The Jacobi preconditioner is M = D. It is extremely cheap to apply, fully parallel, and often removes poor scaling between coordinates. Its main limitation is that it only uses diagonal information.

80

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS

Example: Jacobi for the 1D discrete Laplacian. Consider now the tridiagonal matrix arising from the standard finite-difference discretisation of −u′′ (x) = f (x),

x ∈ (0, 1),

u(0) = u(1) = 0.

With m interior grid points and mesh size h = 1/(m + 1), the linear system is Ah u = f,

Ah =

1 tridiag(−1, 2, −1) ∈ Rm×m . h2

The diagonal of Ah is constant: D = diag(Ah ) = h22 I. Therefore the Jacobi-preconditioned operator is simply h2 1 D−1 Ah = Ah = tridiag(−1, 2, −1). 2 2 Jacobi therefore rescales every eigenvalue by the same constant factor, and cannot improve the condition number κ(Ah ) ∼ h−2 , which grows quadratically as the mesh is refined. Takeaways Jacobi preconditioning is effective when the main difficulty is simple coordinate scaling, but it cannot improve problems such as the 1D discrete Laplacian where the diagonal is already uniform. In PDEs, the hard part is often not local scaling but global error propagation, which motivates Schwarz methods, coarse spaces, and multigrid. A more powerful algebraic alternative is the incomplete factorisation A ≈ L̃Ũ , in which fill-in arising outside a prescribed sparsity pattern is simply discarded [26].

6.1.3

One framework for the whole chapter

Before turning to better preconditioners it is worth noticing that many of them (and all those discussed in this lecture) are instances of a single construction. Let Ri ∈ Rni ×n be a restriction onto a subspace of the unknowns, let Ai = Ri ARi⊤ be the corresponding local matrix, and set X M −1 = Ri⊤ A−1 (6.1) i Ri . i

Each term solves the problem exactly on one subspace and ignores the rest; the preconditioner is their sum. Applying M −1 costs one small solve per subspace, and the solves are independent of one another which is what makes the construction parallel. Everything in this chapter is a choice of the subspaces in (6.1). Subspaces Ri

Preconditioner

What it cannot reach

single unknowns

Jacobi (§6.1)

anything non-local

overlapping subdomains

additive Schwarz (§6.2)

global, slowly varying error

subdomains + a coarse space

two-level, multigrid (§6.3)

—

The table above is the plan of the chapter, and each row fixes the disadvantages of the one above it. For example, a pure local preconditioner handles local error well and global error not at all (for example Jacobi in its simple and block variants or Schwarz preconditioners). To fix this, one should add a global subspace to catch what the local solves miss.

6.2

Schwarz Domain Decomposition

The algorithms of this section, their convergence theory, and their parallel implementation are treated at length in [9, 35].

81

6.2. SCHWARZ DOMAIN DECOMPOSITION

In the language of (6.1), this section makes the first non-trivial choice of subspaces: instead of correcting each unknown independently as Jacobi does, we take the Ri to restrict to overlapping subdomains of the mesh, so that each local solve of A−1 inverts the operator exactly on a whole i piece of the domain. The resulting preconditioner is additive Schwarz, and the question we must answer about it is the one posed in Section 6.1: how does its convergence behave? It is easiest to answer that question by first studying the corresponding stationary iteration, where the local solves are applied one after another. That iteration is of independent historical interest and its convergence factor exposes the two parameters that will govern the preconditioner as well: the overlap width, and the number of subdomains. Model problem: 1D Laplace equation. We consider −u′′ (x) = f (x),

x ∈ (0, 1),

u(0) = u(1) = 0.

We split the interval into two overlapping subdomains Ω1 = (0, β), Ω2 = (α, 1), 0 < α < β < 1. The overlap is the interval (α, β). The alternating Schwarz iteration is: overlap

Ω1 0

α

Ω2 1

β

Figure 6.2: 1D domain split into two overlapping subdomains.  (n) ′′   −(u1 ) = f in (0, β), (n) u1 (0) = 0,   (n) (n−1) u1 (β) = u2 (β),

 (n) ′′   −(u2 ) = f in (α, 1), (n) u2 (1) = 0,   (n) (n) u2 (α) = u1 (α).

Each subproblem is solved exactly, using boundary data from the latest iterate on the neighbouring subdomain. Proposition 6.2.1 (Contraction of alternating Schwarz in 1D). Let u be the exact solution (n) (n) of the Poisson BVP and let u1 , u2 be the alternating Schwarz iterates on Ω1 = (0, β) and (n) (n) (n) (n) Ω2 = (α, 1). Define the errors e1 = u|Ω1 − u1 , e2 = u|Ω2 − u2 . Then (n)

−(e1 )′′ = 0

in (0, β),

(n)

−(e2 )′′ = 0

in (α, 1),

and the interface errors satisfy (n)

(n−1)

|e2 (β)| ≤ ρ |e2

(β)|,

ρ=

α(1 − β) < 1. β(1 − α)

Thus the alternating Schwarz method converges geometrically, and the convergence improves as the overlap β − α increases. Proof. The iterates and the exact solution satisfy the same equation with different boundary data, (n) (n) (n) (n−1) so each error is harmonic on its subdomain: −(ei )′′ = 0, with e1 (0) = 0, e1 (β) = e2 (β) (n) (n) (n) and e2 (1) = 0, e2 (α) = e1 (α). In one dimension, a harmonic function is affine, so each error is determined by its two boundary values: (n)

e1 (x) =

x (n−1) e (β), β 2

(n)

e2 (x) =

1 − x (n) e (α). 1−α 1 (n)

(n−1)

Evaluating the first at x = α and substituting into the second at x = β gives e2 (β) = ρ e2 with ρ = α(1−β) β(1−α) , and 0 < α < β < 1 forces ρ < 1.

(β)

82

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS

Remark 6.2.2. The formula makes the role of the overlap explicit: as β − α → 0 we get ρ → 1 and convergence stalls, while a wider overlap decreases ρ. Overlap is the mechanism by which information is shared by the neighbouring subdomains.

interface error magnitude

small overlap (ρ ≈ 0.7) 100

medium overlap (ρ ≈ 0.4)

large overlap (ρ ≈ 0.2)

4 5 6 Schwarz iteration n

8

10−1 10−2 10−3 10−4

0

1

2

3

7

9

10

Figure 6.3: Alternating Schwarz converges geometrically in 1D. The contraction factor improves when the overlap increases.

Discrete formulation. Discretising with m = 9 interior points and taking β = 7h, α = 4h gives local matrices A1 ∈ R6×6 and A2 ∈ R5×5 , each the 1D Laplacian on its own subdomain. The iteration becomes (n)

A1 u1

(n−1)

= f1 + B̃12 u2

,

(n)

A2 u2

(n)

= f2 + B̃21 u1 ,

where B̃12 , B̃21 inject the interface values from one subdomain into the other (Dirichlet transmission). Note that only A1 and A2 are factorised: the global matrix is never formed. This makes the method efficient at a larger scale. The discrete contraction factor inherits the continuous one: the iteration converges geometrically at a rate set by the overlap and, crucially, independently of h such that refining the mesh does not slow it down, because the overlap is a fixed fraction of the domain. Takeaways The alternating Schwarz iteration quickly damps high-frequency errors and converges independently of the mesh size. It provides a prototype for overlapping domain decomposition preconditioners used in large-scale PDE solvers. From subdomain solves to a preconditioner. The discrete Schwarz viewpoint becomes clearer if we start from the global residual r(n) = f − Au(n) . Instead of trying to correct all unknowns at once, we restrict this residual to each subdomain and solve a local correction problem there. Let Ri be the restriction from the global vector to the unknowns belonging to subdomain Ωi , and let Ai = Ri ARiT be the local stiffness matrix on that subdomain. The local correction δui is defined by Ai δui = Ri r(n) ,

so that

(n) δui = A−1 . i Ri r

This means that on each subdomain we solve exactly the local error equation driven by the restricted residual. The correction is then prolongated back to the global space by RiT . If the

83

6.2. SCHWARZ DOMAIN DECOMPOSITION

subdomains overlap, several local corrections contribute to the same global node, so we blend them with diagonal weights Di satisfying the partition-of-unity condition X RiT Di Ri = I. i

This leads to the global correction X X (n) u(n+1) = u(n) + RiT Di δui = u(n) + RiT Di A−1 . i Ri r i

i

Additive Schwarz preconditioner. The operator X −1 MAS = RiT Di A−1 i Ri i

is called the additive Schwarz preconditioner. With this notation, one step of preconditioned Richardson iteration reads −1 u(n+1) = u(n) + MAS (f − Au(n) ). So the preconditioner acts as an approximate inverse of A obtained by: 1. restricting the residual to each subdomain, 2. solving local problems there, 3. prolongating the local corrections, 4. blending the overlap contributions. Remark 6.2.3. This additive preconditioner is closely related to the alternating Schwarz iteration, but it is not exactly the same algorithm. The alternating method is multiplicative: the second subdomain uses the updated information coming from the first. By contrast, additive Schwarz computes all local corrections from the same global residual and sums them together. The additive form is especially convenient as a preconditioner inside Krylov methods . Key idea Schwarz preconditioning improves on Jacobi by replacing pointwise updates with subdomain solves. Instead of using only diagonal information, it inverts the operator locally on each subdomain. This captures nearest-neighbour couplings much better and yields a significantly stronger preconditioner for discretised PDEs. Example of a decomposition. To make the link with the preconditioned formulation explicit, consider again m = 9 interior nodes and the decomposition Ω1 = (0, 7h),

Ω2 = (4h, 1),

which gives an overlap of three grid points (x5 , x6 , x7 ). The global unknown vector is u = (u1 , u2 , . . . , u9 )⊤ , and the local vectors on each subdomain are u1 = (u1 , . . . , u6 )⊤ ,

u2 = (u5 , . . . , u9 )⊤ .

The restriction matrices R1 ∈ R6×9 and R2 ∈ R5×9 are the corresponding rows of the identity: R1 = [ I6 0 ] picks out u1 , . . . , u6 , and R2 = [ 0 I5 ] picks out u5 , . . . , u9 . They are never stored as matrices since applying Ri is an indexing operation.

84

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS Ω1 D1

Ω2 1 2

1 x1

x2

x3

x4

x5

x6

x7

x8

1 2

D2

x

1 D1 = 1 D1 = 12

α = 4h

x9

D2 = 1 D2 = 12

β = 7h

Figure 6.4: Partition of unity in the overlap: D1 (top) and D2 (bottom). On the overlap, D1 = D2 = 12 and outside the overlap they sum to one: R1T D1 R1 + R2T D2 R2 = I. Each restriction operator extracts local unknowns, and its transpose RiT prolongates the local corrections into the global vector. In the overlap region (nodes x5 , x6 , x7 ), contributions from both subdomains are combined using a discrete partition of unity: D1 = diag(1, 1, 1, 1, 12 , 21 , 21 , 0, 0), These satisfy

D2 = diag(0, 0, 0, 0, 12 , 12 , 21 , 1, 1).

R1T D1 R1 + R2T D2 R2 = I,

ensuring that the overlapping corrections are weighted consistently and the global operator remains a proper approximation of A−1 . The additive Schwarz preconditioner then reads −1 T M −1 = R1T D1 A−1 1 R1 + R2 D2 A2 R2 , which applies local solves and blends them smoothly through the partition weights.

6.3

Two-Level Preconditioning and Multigrid Ideas

In order to understand the mechanism of error reduction or smoothing and then the need of a coarse space correction let us consider the discretisation of the 1D Poisson problem. The eigenvectors of Ah are the discrete sine modes  jπℓ  , ℓ = 1, . . . , m, j = 1, . . . , m, vj (ℓ) = sin m+1 with eigenvalues λj (Ah ) =

  jπ 4 2 sin . h2 2(m + 1)

Small values of j correspond to smooth, slowly varying modes, whereas large values of j correspond to oscillatory modes.

smooth
mode πj sin m+1 , small j

oscillatory
mode πj sin m+1 , large j

Figure 6.5: Low-frequency and high-frequency error modes on a 1D grid. Consider the Jacobi iteration u(k+1) = u(k) + D−1 (f − Ah u(k) ),

D=

2 I. h2

6.3. TWO-LEVEL PRECONDITIONING AND MULTIGRID IDEAS

85

If e(k) = u − u(k) denotes the error, then TJ = I − D−1 Ah .

e(k+1) = TJ e(k) ,

Proposition 6.3.1 (Mode-by-mode damping of Jacobi). For the 1D Laplacian, each sine mode vj is an eigenvector of the Jacobi error propagator TJ , with eigenvalue µj = 1 −

   jπ  λj (Ah ) jπ 2 = 1 − 2 sin = cos . 2/h2 2(m + 1) m+1

Hence, if e(0) =

m X

cj vj ,

j=1

then after k Jacobi steps e(k) =

m X

cj µkj vj .

j=1

Proof. Since the vj are eigenvectors of Ah , Ah vj = λj (Ah )vj . Therefore

   λj (Ah )  TJ vj = I − D−1 Ah vj = 1 − vj . 2/h2

Using λj (Ah ) = we obtain

  jπ 4 2 sin , h2 2(m + 1)

 µj = 1 − 2 sin2

 jπ . 2(m + 1)

The trigonometric identity 1 − 2 sin2 (θ) = cos(2θ) gives  jπ  µj = cos . m+1 The expansion of e(k) follows by diagonalising the iteration on the basis {vj }m j=1 . Remark 6.3.2. This formula shows explicitly why Jacobi behaves differently on different scales: j small whereas

j large

=⇒

=⇒

|µj | ≈ 1

(slow damping),

|µj | is much smaller

(fast damping).

So Jacobi is not a good solver by itself, but it is a good smoother.

6.3.1

Weighted Jacobi and the smoothing factor

A common variant is weighted Jacobi: u(k+1) = u(k) + ωD−1 (f − Ah u(k) ). Its error propagator is

TJ,ω = I − ωD−1 Ah .

86

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS

For the sine mode vj , (ω)

(ω)

TJ,ω vj = µj vj ,

µj

 = 1 − 2ω sin2

 jπ . 2(m + 1)

To study smoothing, we focus on the upper half of the spectrum, m+1 2 ≤ j ≤ m, which corresponds to oscillatory modes. For these indices,  sin2

 h1 i jπ ∈ ,1 , 2(m + 1) 2

(ω)

∈ [ 1 − 2ω, 1 − ω ].

hence µj

So the worst-case damping factor on the high-frequency range is
ρhigh (ω) = max |1 − 2ω|, |1 − ω| . Proposition 6.3.3 (Optimal weighted Jacobi parameter for smoothing). For the 1D Laplacian, the value of ω that minimizes the worst-case high-frequency damping factor is ω = 23 . For this choice, ρhigh = 13 . Thus every high-frequency mode is reduced by at least a factor 1/3 per step.

Proof. We minimize ρhigh (ω) = max(|1 − 2ω|, |1 − ω|). The optimum is obtained by balancing the two endpoint magnitudes: |1 − 2ω| = |1 − ω|. For ω ∈ (1/2, 1), this becomes 2ω − 1 = 1 − ω, hence ω = 23 . Then ρhigh = 13 . Takeaways Weighted Jacobi with ω = 2/3 is not a scalable solver, but it is an effective smoother: it damps oscillatory components uniformly while leaving smooth components for the coarse grid to remove. After several relaxation steps, the remaining error typically has the form e = esmooth + eosc ,

∥eosc ∥ ≪ ∥esmooth ∥.

That is, the oscillatory part has already been reduced, while the smooth part remains. The key multigrid observation is that a smooth error on the fine grid may look much more oscillatory on a coarser grid. This makes it easier to represent and eliminate there. Key idea A good solver for an elliptic problem must act on two scales at once: local smoothing

+

global coarse correction.

6.3. TWO-LEVEL PRECONDITIONING AND MULTIGRID IDEAS

87

after 5 smoothing sweeps initial error piecewise-linear interpolant of its coarse-node values

fine nodes

coarse nodes (H = 2h)

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 = 31; the thick blue curve is the same error after five weighted-Jacobi sweeps (ω = 23 ); the dashed orange curve is the piecewise-linear interpolant through its values at the 15 coarse nodes only.

6.3.2

Coarse-grid correction

Section 6.3.1 left us with an iterate whose error is smooth: the smoother has removed the oscillatory components and stalls on the rest. What remains satisfies the error equation Ae = r,

r = f − Au,

e = u∗ − u,

and solving that on the fine grid is of course just the original problem again. The way out rests on a single observation: a smooth vector is cheap to represent. It is determined, to good accuracy, by its values at half of the grid points, the rest being recovered by interpolation so the error equation may be solved in a space of half the dimension. Formally, let P : Rnc → Rn be the prolongation (linear interpolation from the coarse grid), let R : Rn → Rnc be the restriction, and define the coarse operator by AH = RAP. Taking R ∝ P ⊤ makes AH the Galerkin projection of A onto the coarse space; for the 1D Laplacian with full weighting R = 12 P ⊤ it coincides exactly with the direct discretisation on the coarse grid (Problem 4). The correction then reads: restrict the residual, solve AH eH = Rr, and prolongate, u ← u + P eH . Coarse-grid correction cannot be an iterative method on its own. Its error propagation operator is T = I − P A−1 H RA, and a short computation using AH = RAP shows that T 2 = T (Problem 5). A projection has eigenvalues 0 and 1 only, so ρ(T ) = 1: whatever lies in ker T is annihilated in one step, and whatever does not is left completely untouched, forever. Coarse-grid correction removes the part of the error that the coarse space can see, and by construction can do nothing about the rest. This is the precise sense in which the two ingredients are complementary, and it is why neither is optional. The smoother cannot touch smooth error; the coarse correction cannot touch anything outside its range. Each ingredient is useless on its own; the next subsection shows that together they are optimal.

88

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS two-grid S ν1

V-cycle S ν2

smooth

smooth

Ωh R

P

Ω2h exact solve

Ω4h coarsest solve

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) work and reduces the error by a factor independent of h. Used as a preconditioner for CG, this turns the O(h−1 ) iteration count of Chapter 4 into a constant.

6.3.3

From coarse correction to the two-grid method

The two ingredients are now in place, and Section 6.3.2 showed that each fails precisely where the other succeeds. Combining them in the obvious order gives the two-grid method : smooth to remove what the coarse space cannot see, correct on the coarse grid to remove what the smoother cannot touch, then smooth again to clean up the interpolation error introduced by P . Starting from an approximation u: 1. Pre-smoothing: u ← S ν1 (u), with S weighted Jacobi or Gauss–Seidel. 2. Residual restriction: rH = R(f − Au). 3. Coarse solve: AH eH = rH . 4. Coarse correction: u ← u + P eH . 5. Post-smoothing: u ← S ν2 (u). From two grids to multigrid. Step 3 still requires an exact solve with n/2 unknowns (n/4 in 2D), which may itself be far too large to solve directly. But the observation that motivated the coarse grid applies again one level down: once the coarse problem has itself been smoothed, its error is smooth relative to H and is representable on a grid of spacing 2H. So replace the exact coarse solve by one two-grid step on Ω2h , its coarse solve by a two-grid step on Ω4h , and so on until the grid is small enough to solve directly. That recursion is multigrid, and the grids visited on the way down and back up form the Vcycle of Figure 6.7. Two properties follow. The work is O(n), since the grids shrink geometrically; and the error reduction per cycle is independent of h, being governed by the smoothing factor of Section 6.3.1 e.g. 13 per step for weighted Jacobi in 1D rather than by κ(A). Takeaways Smoother and coarse correction are exactly complementary [7, 16]: the first removes oscillatory error and cannot touch smooth error, the second removes what lies in the coarse space and cannot touch anything else. Recursion turns the pair into an O(n) solver whose contraction factor does not depend on h.

6.3. TWO-LEVEL PRECONDITIONING AND MULTIGRID IDEAS

89

Multigrid as a preconditioner In practice, multigrid is often used not as a standalone solver, but as a preconditioner inside a Krylov method. One multigrid cycle is then viewed as an approximate inverse: −1 MMG ≈ A−1 . The linear system is solved in preconditioned form, for example by CG or GMRES: −1 −1 MMG Ax = MMG b.

This fits naturally into the general framework of preconditioning: multigrid constructs a strong approximate inverse by combining local relaxation with a global coarse-space correction. Remark 6.3.4. For elliptic problems, multigrid-preconditioned Krylov methods often have iteration counts that remain essentially bounded under mesh refinement. This is the sense in which multigrid is called mesh-independent or even optimal : the number of iterations no longer grows significantly as the grid is refined.

90

6.4

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS

Exercises

Problem 1 (Convergence of a parallel Schwarz method in 1D, past exam). Consider the one-dimensional Laplace problem − u′′ (x) = f (x),

x ∈ (0, 1),

u(0) = u(1) = 0,

and the overlapping decomposition Ω = Ω1 ∪ Ω2 with Ω1 = (0, α),

Ω2 = (β, 1),

0 < β < α < 1,

The parallel (additive) Schwarz iteration is ( (n+1)
′′ − u1 =f in (0, α), (n+1) u1 (0) = 0, (n)

where ui

δ = α − β > 0.

( (n+1)
′′ =f − u2

(n+1) (n) u1 (α) = u2 (α),

(n+1) (n) u2 (β) = u1 (β),

in (β, 1), (n+1) u2 (1) = 0,

is the nth iterate restricted to Ωi . (n)

(n)

(a) Let ei := ui − ui denote the errors on each subdomain (where u is the exact solution). (n+1) (n+1) Write the boundary value problems (BVPs) satisfied by e1 and e2 . (b) Solve these BVPs explicitly (they are linear functions on each subdomain). (c) Let (n)

(n)

ηn := e2 (α),

ζn := e1 (β).

Show that the interface errors satisfy the two-step recursions ηn+1 = and hence

(n)

e1 (α) (n−2) e1 (α)

1−α ζn , 1−β

ζn+1 =

(n)

=

e2 (β) (n−2) e2 (β)

= ρ,

ρ :=

β ηn , α β 1−α · . α 1−β

(d) Prove that 0 < ρ < 1 whenever β < α (i.e., when there is positive overlap), and conclude that the parallel Schwarz method is convergent. Problem 2 (Spectral equivalence ⇒ PCG convergence). Let A, B be SPD and suppose there exist c1 , c2 > 0 such that c1 x⊤ Bx ≤ x⊤ Ax ≤ c2 x⊤ Bx

∀x ∈ Rn .

(a) Show that spec(B −1/2 AB −1/2 ) ⊂ [c1 , c2 ], hence κ(B −1 A) ≤ c2 /c1 . (b) Conclude that PCG applied to Au = b with left preconditioner B satisfies !k p κ(B −1 A) − 1 ∥ek ∥H ≤ 2 p , ek = uk − u∗ . ∥e0 ∥H κ(B −1 A) + 1 with H = B −1/2 AB −1/2 . Problem 3 (Preconditioned conjugate gradient). Let A, M be SPD. Apply CG to Ãx̃ = b̃ with à = M −1/2 AM −1/2 and b̃ = M −1/2 b. Show that the mapped iterates x(k) = M −1/2 x̃(k) satisfy the iteration in the proposition; in particular p(0) = M −1 r(0) , and αk =

⟨r(k) , z (k) ⟩ , ⟨p(k) , Ap(k) ⟩

βk+1 =

⟨r(k+1) , z (k+1) ⟩ , ⟨r(k) , z (k) ⟩

z (k) = M −1 r(k) .

91

6.4. EXERCISES Problem 4 (Transfer operators and the Galerkin coarse matrix). Let Ah =

1 tridiag(−1, 2, −1) ∈ Rm×m h2

with m = 2mH + 1, h = 1/(m + 1), and let P ∈ Rm×mH be linear interpolation from the coarse grid H = 2h: (P eH )2j = 12 (eH )j , (P eH )2j+1 = (eH )j , (P eH )2j+2 = 21 (eH )j . (a) Write P explicitly for m = 7, mH = 3, and show that the full weighting restriction R = 12 P ⊤ applies the stencil 14 (1, 2, 1). (c) Compute AH = RAh P for m = 7 and verify that it coincides exactly with the direct discretisation H12 tridiag(−1, 2, −1) on the coarse grid. (d) What do you get from P ⊤ Ah P instead, and why must the scalings of R and P be consistent? Problem 5 (Coarse-grid correction alone cannot converge). With Ah , P, R and AH = RAh P as in Problem 5, the coarse-grid correction u ← u+P A−1 H R(f −Ah u) has error propagation operator T = I − P A−1 RA . h H (a) Show that T 2 = T (use AH = RAh P to simplify P A−1 H RAh P ), and deduce that every eigenvalue of T is 0 or 1, that ρ(T ) = 1 when mH < m, and that rank T = m − mH . (c) Conclude that coarse-grid correction alone is not convergent, identify which errors lie in ker T (compare Figure 6.6), and explain in one sentence why a smoother is therefore mandatory rather than merely helpful. Problem 6 (The two-grid contraction factor). Let S = I −ωD−1 Ah be the weighted Jacobi smoother, D = diag(Ah ), and T as in Problem 6; the two-grid operator with ν pre-smoothing steps is MTG = T S ν . jπ
(a) From the eigenvalues µj (ω) = 1 − 2ω sin2 2(m+1) of S, explain why ω = 1 is a poor smoother even though it minimises ρ(S). (b) Take ω = 23 and m = 7. Compute ρ(MTG ) numerically for ν = 1 and ν = 2, and compare with the smoothing factor 13 of Proposition 6.3.3. State the relationship you observe between ρ(MTG ) and ν. (c) Repeat for m = 15 and m = 31: is the contraction factor independent of h? Contrast with the O(h−2 ) iteration count of unpreconditioned Jacobi in Chapter 3. (d) A V-cycle visits grids of size m, m/2, m/4, . . . ; show its total cost is O(m), and explain why this together with (c) makes multigrid an optimal solver.

92

CHAPTER 6. PRECONDITIONING AND ACCELERATING SOLVERS

Bibliography [1] Ahmad Abdelfattah, Hartwig Anzt, Erik G Boman, Erin Carson, Terry Cojean, Jack Dongarra, Alyson Fox, Mark Gates, Nicholas J Higham, Xiaoye S Li, Jennifer Loe, Piotr Luszczek, Srikara Pranesh, Siva Rajamanickam, Tobias Ribizel, Barry F Smith, Kasia Swirydowicz, Stephen Thomas, Stanimire Tomov, Yaohung M Tsai, and Ulrike Meier Yang. A survey of numerical linear algebra methods utilizing mixed-precision arithmetic. The International Journal of High Performance Computing Applications, 35(4):344–369, July 2021. [2] Walter Edwin Arnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue problem. Quarterly of Applied Mathematics, 9(1):17–29, 1951. [3] Owe Axelsson. Iterative Solution Methods. Cambridge University Press, Cambridge, 1994. [4] Michele Benzi. Preconditioning techniques for large linear systems: A survey. Journal of Computational Physics, 182(2):418–477, 2002. [5] Ankit Bisain, Alan Edelman, and John Urschel. A new upper bound for the growth factor in Gaussian elimination with complete pivoting. Bulletin of the London Mathematical Society, 57(5):1369–1387, 2025. [6] R. Bouyouli, G. Meurant, L. Smoch, and H. Sadok. New results on the convergence of the conjugate gradient method. Numerical Linear Algebra with Applications, 16(3):223–236, 2009. [7] William L. Briggs, Van Emden Henson, and Steve F. McCormick. A Multigrid Tutorial. SIAM, Philadelphia, PA, 2nd edition, 2000. √ [8] Michel Crouzeix and César Palencia. The numerical range is a (1 + 2)-spectral set. SIAM Journal on Matrix Analysis and Applications, 38(2):649–655, 2017. [9] Victorita Dolean, Pierre Jolivet, and Frédéric Nataf. An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation. SIAM, Philadelphia, PA, 2015. [10] Miroslav Fiedler. Algebraic connectivity of graphs. Czechoslovak Mathematical Journal, 23(2):298–305, 1973. [11] Gene Golub and William Kahan. Calculating the singular values and pseudo-inverse of a matrix. Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, 2(2):205–224, 1965. [12] Gene H. Golub and Charles F. Van Loan. Matrix Computations - 4th Edition. Johns Hopkins University Press, Philadelphia, PA, 2013. [13] Anne Greenbaum. Iterative Methods for Solving Linear Systems. SIAM, Philadelphia, PA, 1997. 93

94

BIBLIOGRAPHY

[14] Anne Greenbaum, Vlastimil Pták, and Zdenve˘k Strakoš. Any Nonincreasing Convergence Curve is Possible for GMRES. SIAM Journal on Matrix Analysis and Applications, 17(3):465– 469, July 1996. [15] S. A. Haben. Conditioning and Preconditioning of the Minimisation Problem in Variational Data Assimilation. PhD thesis, Department of Mathematics and Statistics, University of Reading, 2011. [16] Wolfgang Hackbusch. Multi-Grid Methods and Applications, volume 4 of Springer Series in Computational Mathematics. Springer, Berlin, 1985. [17] Martin Hanke. Conjugate Gradient Type Methods for Ill-Posed Problems, volume 327 of Pitman Research Notes in Mathematics Series. Longman Scientific & Technical, Harlow, 1995. Reissued by Chapman & Hall/CRC, Boca Raton, 2018. [18] Per Christian Hansen. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. SIAM, Philadelphia, PA, 1998. [19] Sam Hatfield, Andrew McRae, Tim Palmer, and Peter Düben. Single-Precision in the Tangent-Linear and Adjoint Models of Incremental 4D-Var. Monthly Weather Review, 148(4):1541–1552, April 2020. [20] Taher H. Haveliwala and Sepandar D. Kamvar. The second eigenvalue of the Google matrix. Technical Report 2003-20, Stanford University, 2003. [21] Magnus R. Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6):409–436, 1952. [22] Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. Acta Numerica, 31:347–414, May 2022. [23] Cornelius Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Journal of Research of the National Bureau of Standards, 45(4):255–282, 1950. [24] Amy N. Langville and Carl D. Meyer. Google’s PageRank and Beyond: The Science of Search Engine Rankings. Princeton University Press, Princeton, NJ, 2006. [25] M. J. Martin, A. Hines, and M. J. Bell. Data assimilation in the FOAM operational shortrange ocean forecasting system: A description of the scheme and its impact. Quarterly Journal of the Royal Meteorological Society, 133(625):981–995, 2007. [26] J. A. Meijerink and H. A. van der Vorst. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Mathematics of Computation, 31(137):148–162, 1977. [27] Yurii Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Springer, Boston, MA, 2004. [28] Alexander M Ostrowski. On the linear iteration procedures for symmetric matrices. Rend. Mat. Appl., 14:140–163, 1954. [29] Lawrence Page, Sergey Brin, Rajeev Motwani, and Terry Winograd. The PageRank citation ranking: Bringing order to the web. Technical Report 1999-66, Stanford InfoLab, 1999. [30] Christopher C. Paige. Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem. Linear Algebra and its Applications, 34:235–258, 1980.

BIBLIOGRAPHY

95

[31] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the Spectral Bias of Neural Networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of PMLR, pages 5301– 5310, 2019. [32] Yousef Saad. On the rates of convergence of the Lanczos and the block-Lanczos methods. SIAM Journal on Numerical Analysis, 17(5):687–706, 1980. [33] Yousef Saad. Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia, PA, 2nd edition, 2003. [34] Yousef Saad and Martin H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986. [35] Andrea Toselli and Olof Widlund. Domain Decomposition Methods: Algorithms and Theory, volume 34 of Springer Series in Computational Mathematics. Springer, Berlin, 2005. [36] Lloyd N. Trefethen and Mark Embree. Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press, Princeton, NJ, 2005. [37] Richard S. Varga. Matrix Iterative Analysis. Springer, Berlin, Heidelberg, 2nd edition, 2000. [38] Ulrike von Luxburg. A tutorial on spectral clustering. Statistics and Computing, 17(4):395– 416, 2007. [39] Andrew J. Wathen. Preconditioning. Acta Numerica, 24:329–376, 2015. [40] J. H. Wilkinson. The Algebraic Eigenvalue Problem. Clarendon Press, 1965. [41] David Young. On Richardson’s Method for Solving Linear Systems with Positive Definite Matrices. Journal of Mathematics and Physics, 32(1-4):243–255, 1953. [42] David M. Young. Convergence properties of the symmetric and unsymmetric successive overrelaxation methods and related methods. Mathematics of Computation, 24(112):793–807, 1970.

Related documents

• Invertible continuous latent dynamic for long-term data assimilation in complex physical systems#1141476
• Erratum to: Computation of pseudospectral abscissa for large-scale nonlinear eigenvalue problems#955431
• Residual-Controlled Douglas--Rachford Splitting for Differentiable Solver Layers#609633
• Generative machine learning for multivariate angular simulation.#281019
• Structured Divide-and-Conquer for the Definite Generalized Eigenvalue Problem#625472