PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Victorita Dolean, Jemima Tabeart

arXiv:2608.21234v1math.NAcs.LG

TL;DR

Ill-conditioned, noisy least-squares problems arise in machine learning and other inverse tasks, where solving to machine precision can reproduce noise rather than the underlying model. These notes connect standard numerical linear algebra methods to such applications and show that early stopping conjugate gradient is itself a form of regularisation.

  • Problem

    Ill-conditioned, noisy least-squares problems arise in machine learning and other inverse tasks, where solving to machine precision can reproduce noise rather than the underlying model.

    Source: §4.3 · p. 61

  • Method

    The notes present classical numerical linear algebra algorithms and connect Krylov methods to machine learning applications.

    Source: Preface · p. 5, Overview · p. 51

  • Results

    Ridge regression, truncated SVD, and CG are spectral filters of the same data.

    Source: Overview · p. 51

  • Takeaways & Limitations

    Early stopping CG is itself a regularisation for machine learning.

    Source: Overview · p. 51

    As iteration increases, the filter admits more noise, so a stopping rule is indispensable.

    Source: §4.3 · p. 63

Abstract

from arXiv · show

These lecture notes form the first part of a master's-level course on advanced numerical linear algebra. Their aim is not only to present the classical algorithms, but to show why the subject has become considerably more central than it was a generation ago. Numerical linear algebra grew up alongside the numerical solution of partial differential equations, and for a long time that is where its large sparse systems came from. Ranking the nodes of a network, assimilating observations into a weather forecast, and fitting a model to a large noisy data set now lead to problems of the same kind: too large to factorise, structured, and accessible only through matrix-vector products. Strikingly few ideas are needed for all of them. Each chapter therefore develops a standard topic and then puts it to work outside its original setting. We treat norms, factorisations, conditioning and floating-point arithmetic; sparse matrices arising from finite differences, from graphs and from machine learning; stationary iterations and the smoothing property; the conjugate gradient and Lanczos methods, with spectral clustering and regularisation by early stopping; Arnoldi and GMRES, with PageRank and large least squares; and finally preconditioning, Schwarz domain decomposition and multigrid. We assume a first course in linear algebra. Every section closes with a summary of what should be retained and every chapter with exercises, several drawn from past examinations. Accompanying Python code reproduces the numerical illustrations.

Solving Linear Systems: A Unified Foundation · Overview

This section establishes numerical linear algebra foundations for solving systems, covering norms, conditioning, stability, finite-precision arithmetic, direct methods, and stationary iterations, including their costs and convergence properties. It also includes exercises on induced norms and stationary iterative methods, with a SOAR-matrix implementation task.

Overview

  • Linear algebra underlies scientific and engineering computations, and this chapter introduces norms, conditioning, stability, LU and Cholesky factorizations, and stationary iterative methods.Its motivating applications include structural analysis, machine learning, and weather data assimilation.Source: Overview · p. 7

1.1.1 Vector and Matrix Norms

  • Vector norms measure size, distance, convergence, and approximation error; all norms are equivalent in finite dimensions, though their values can differ.For example, ∥x∥∞ ≤ ∥x∥2 ≤ √n∥x∥∞.Source: §1.1.1 · p. 8

1.1.2 Eigenvalues and Matrix Decompositions

  • Eigenvalues and the spectral radius govern long-run iterative behavior, while Schur decomposition and SVD exist for every matrix.An eigendecomposition requires a diagonalizable matrix; the SVD uses nonnegative singular values.Source: §1.1.2 · p. 9

1.1.3 Matrix Norms

  • Induced matrix norms support analysis of algorithm and problem stability and conditioning; the spectral norm equals the largest singular value.Norms can agree on asymptotic convergence while differing numerically.Source: §1.1.3 · pp. 9–10

1.2. FLOATING-POINT ARITHMETIC, ERROR, CONDITIONING, AND FACTORIZATIONS5

  • Reduced-precision arithmetic can improve efficiency but degrade performance, motivating mixed-precision methods that seek computational gains while maintaining convergence and stability.The notes describe growing interest in reduced precision driven in part by machine-learning applications.Source: §1.2 · p. 11

Sources of Numerical Errors

  • Rounding, truncation, and cancellation cause numerical error; cancellation can erase significant digits, floating-point addition is non-associative, and lower precision raises overflow and rounding risks.In catastrophic cancellation, subtracting nearly equal stored values can yield an inaccurate result or even zero.Source: Sources of Numerical Errors · p. 11

1.2.2 Error, Conditioning and Stability of algorithms

  • Ill-conditioned systems can amplify small perturbations, and final accuracy depends jointly on problem conditioning and algorithm stability.A stable algorithm can still produce inaccurate results on an ill-conditioned problem.Source: §1.2.2 · p. 12, §1.2 · p. 13

1.2.3 Gaussian Elimination and LU Factorization

  • Gaussian elimination costs O(n^3) flops, while forward and backward substitution each cost O(n^2), motivating iterative solvers and preconditioners for large systems.The notes identify O(n^3) as infeasible for systems with n ∼ 10^6.Source: §1.2.3 · p. 13
  • GEPP’s row pivoting is backward stable, while elimination without pivoting can fail even for well-conditioned A; O(n^3) LU limits very large sparse systems.Iterative solvers are introduced to overcome that scaling barrier.Source: §1.2.3 · p. 14

1.2.4 Cholesky Factorization for SPD Matrices

  • For SPD A, Cholesky is backward stable without pivoting, uses half LU’s storage, and costs approximately n^3/3 flops; a nonpositive pivot calls SPD status or scaling into question.The notes recommend switching to LDL^T with symmetric pivoting or regularizing A with λI if a nonpositive pivot occurs.Source: §1.2.4 · p. 14

1.3 Stationary Iterative Methods

  • A stationary iteration converges for every initial guess exactly when ρ(G) < 1; iterative updates cost O(nnz(A)) and preserve sparsity.Factorization fill-in can expand storage from O(nnz(A)) by orders of magnitude.Source: §1.3 · p. 15

1.3.1 Examples of Stationary Methods and their convergence

  • Jacobi uses previous-iterate values, Gauss–Seidel reuses fresh values, and SOR extrapolates updates; Jacobi and Gauss–Seidel converge under strict row diagonal dominance, while Gauss–Seidel also converges for SPD A.Jacobi and Gauss–Seidel are not generally comparable: examples exist where one converges and the other does not.Source: §1.3.1 · p. 16

1.3. STATIONARY ITERATIVE METHODS

  • For SPD A, SOR converges when 0 < ω < 2, while Richardson converges iff 0 < α < 2/λmax(A), with optimal step 2/(λmax(A) + λmin(A)).Jacobi, Gauss–Seidel, and SOR remain useful as multigrid smoothers and Krylov preconditioners despite being obsolete as standalone large-system solvers.Source: §1.3 · pp. 17–18

1.4 Exercises

  • Exercises ask students to compare condition numbers and norms, explore cancellation and summation order, and assess conditioning and residual-based backward error.Source: §1.4 · pp. 19–20
  • The exercises ask students to relate induced norms to spectral radius, prove Gelfand’s formula, and establish the Neumann-series result when ρ(B) < 1.The prompts also ask students to show that I − B is invertible.Source: §1.4 · p. 20
  • A basic-iteration exercise covers deriving and performing Jacobi steps, proving convergence via ρ(G) < 1, and implementing Jacobi, Gauss–Seidel, and SOR on a SOAR correlation matrix.The SOAR task specifies n = 100, L = 0.4, and a = 1, and asks how numerical results connect to theory; no results are supplied.Source: §1.4 · p. 20

Discretization of 1D and 2D PDEs (and Beyond) · Overview

This section examines how discretizations and related models produce sparse matrices, and how their structure, storage, and reordering affect computational cost. It connects spectral properties to conditioning and iterative convergence, extending sparsity and spectral analysis across PDE matrices, graph Laplacians, data-assimilation correlations, and document-term matrices.

Overview

  • The chapter links PDE, graph, and machine-learning matrices through local structure, and introduces sparse storage, bandwidth, reordering, fill-in, and spectral properties relevant to iterative-solver convergence and preconditioning.It also covers CSR/CSC and banded storage formats.Source: Overview · p. 21

2.1 Finite difference approximation

  • The three-point central-difference approximation to u′′ is consistent with truncation error O(h^2), while central differences for u′ and u′′ are second order and forward or backward differences for u′ are first order.Higher-order formulas can be obtained by eliminating higher derivatives from longer Taylor stencils.Source: §2.1 · pp. 22–23

2.1 FINITE DIFFERENCE APPROXIMATION

  • For the 1D Poisson boundary-value problem, central differences at n interior grid points yield a tridiagonal system, with Dirichlet boundary values contributing to the endpoint right-hand sides.The matrix sparsity pattern corresponds to the stencil [−1, 2, −1]/h^2.Source: §2.1 · pp. 23–24
  • Finite-difference discretization turns differential operators into sparse matrices whose nonzero patterns follow the stencil, which also fixes bandwidth and direct-solver cost.For the 1D Poisson problem, the matrix has the stencil pattern [−1, 2, −1]/h^2.Source: §2.1 · p. 24

2.2 Sparse Matrices: PDEs, Graphs, and Machine Learning

  • Sparse matrices from PDEs and graphs are central to modern scientific computing and machine learning, with graph Laplacians extending the continuous Laplace operator to discrete networks.Their applications include graph neural-network message passing and electrical networks governed by Kirchhoff’s laws.Source: §2.2 · pp. 24, 27

2.2.1 The 2D Poisson equation: structure, stencils and blocks

  • The 2D Poisson discretization forms a block-tridiagonal system: each diagonal block is tridiagonal, each row has at most five nonzeros, and the matrix is SPD under Dirichlet data.Inhomogeneous Dirichlet data modify the right-hand side near the boundary without changing sparsity; one-sided differences for Neumann or Robin conditions modify only the first or last block rows.Source: §2.2.1 · p. 24, §2.2 · p. 25

2.2.2 Graph Laplacians

  • The graph Laplacian L = D − A is symmetric, positive semidefinite, and sparse, with nnz(L) = 2|E| + n; its eigenvectors support spectral clustering.Graph Laplacians also serve as discrete Laplace operators on networks, including diffusion operators used as spatial covariance matrices in data assimilation.Source: §2.2.2 · pp. 25–26, §2.2 · p. 27

2.2.3 Sparse Matrices in Machine Learning

  • Machine-learning sparsity arises from locality: when each feature-matrix row has k ≪ n nonzeros, nnz(A) = O(mk), and related normal equations and Hessians can retain sparsity.In logistic regression, H = A⊤WA; if A is sparse, H is sparse and computations can be performed efficiently.Source: §2.2.3 · pp. 27–28

2.2.4 Sparsity in data assimilation

  • Sample covariance estimates in ensemble data assimilation can create spurious long-range correlations; localisation imposes local sparsity, while fitted long-tail correlations can be truncated below a threshold, typically 0.2.Fitted correlation functions such as Gaussian or autoregressive forms yield dense matrices before truncation.Source: §2.2.4 · p. 28

2.3 Sparse Matrix Toolbox

  • Sparse matrices arise in discretized PDEs, graphs, and machine learning, and computing with them requires compact storage, efficient operations, and attention to spectral properties.Common formats include CSR/CSC and diagonal or banded storage for matrices whose nonzeros cluster near the diagonal.Source: §2.3 · p. 28, §2.3.1 · p. 30

2.3.1 Sparsity Patterns and storage formats

  • For 2D Poisson with N = n^2, bandwidth grows like n, so naive banded LU costs O(N^2), while nested dissection reduces the cost to O(N^3/2).Permuting rows and columns can reduce bandwidth and fill-in, saving orders of magnitude in memory and runtime; a sparse matrix inverse can nevertheless be dense.Source: §2.3.1 · p. 30
  • For sparse matrices, nonzero patterns determine storage and fill-in; reordering can reduce bandwidth and fill, while the inverse of a sparse matrix is not necessarily sparse.The notes contrast sparse LU with a dense inverse in an example and state that effective reordering can save orders of magnitude in memory and runtime.Source: §2.3.1 · p. 31, Figure 2.4 · p. 31

2.3.2 Spectral Properties

  • As the mesh refines, model-problem condition numbers grow like h^-2 and iterative solvers slow; eigenvalue distributions govern convergence, and preconditioning reshapes the spectrum.Sparsity alone does not guarantee fast convergence; clustering can help Krylov methods, while wide spreads or outliers can cause stagnation.Source: §2.3.2 · p. 32

2.4 Exercises

  • Exercises cover pointwise finite-difference accuracy for u(x)=sin(πx), BVP discretization with a basic iterative method, 2D-grid bandwidth/LU estimation at n=20, and exact Laplacian spectra for a four-vertex path and K4.The bandwidth task asks for semi-bandwidth b and banded-LU work O(N b2), with N=n2. The BVP exercises are labeled exam 2024/2025 and resit exam 2024/2025.Source: §2.4 · pp. 34–35
  • The exercises state the 2D stiffness relation bA2D = (1/h²)A2D with standard spacing h = 1/(n + 1), an optional d-dimensional case with nd unknowns, and numerical 1D Poisson eigenvalues and κ(A) for n = 5.Source: §2.4 · p. 35
  • Graph Laplacian exercises connect the multiplicity of eigenvalue 0 to connected components and ask for cycle and star graph spectra.For the star graph Sn, the stated spectrum is {0, 1 (mult n−2), n}; the cycle exercise also asks for an eigenvalue interval.Source: §2.4 · p. 35
  • A data-assimilation exercise uses a SOAR correlation matrix on a unit circle to examine how thresholding small entries affects its spectrum.The matrix may be dense and ill-conditioned; the exercise varies correlation lengthscale, threshold, and problem dimension.Source: §2.4 · p. 36
  • A document-term exercise builds a sparse 4 × 6 binary matrix from four documents, then asks for A⊤A, its nonzero off-diagonal entries, and their interpretation.The documents are bags of words drawn from a six-word vocabulary.Source: §2.4 · p. 36

Applications of Basic Iterative Methods · Overview · 3.1 Basic Iterative Methods for PDE discretisations

The chapter uses explicitly analyzable spectra to make fixed-point convergence concrete, beginning with the one-dimensional Dirichlet Poisson problem and connecting iterative methods to graph and machine-learning applications.

Overview

  • Gauss–Seidel converges roughly twice as fast as Jacobi, while SOR with an optimal ω improves the rate further.Source: Overview · p. 37
  • On graph Laplacians, Jacobi becomes a random-walk or diffusion process whose slow modes relate to the Fiedler vector and spectral clustering.Source: Overview · p. 37
  • Richardson iteration is gradient descent on a quadratic loss, with the spectrum of A or X⊤X governing step-size choice and convergence speed.Source: Overview · p. 37
  • Machine-learning spectral bias uses the same spectral-filtering mechanism as PDE smoothing, but with low- and high-frequency modes swapped.Source: Overview · p. 37

3.1 Basic Iterative Methods for PDE discretisations

  • For the 1D Poisson problem, the stiffness matrix is symmetric positive definite and tridiagonal, with κ(A) = O(n^2).Its extreme eigenvalues scale as λ_min ∼ π^2 and λ_max ∼ 4/h^2.Source: §3.1 · p. 38
  • High-frequency modes are damped quickly by Jacobi, Gauss–Seidel, and SOR, while low-frequency modes decay slowly, yielding the smoothing property used in multigrid.Source: §3.1 · p. 38

3.1.1 Convergence of the Jacobi iteration

  • Jacobi iteration converges very slowly as the number of interior points n grows.Source: §3.1.1 · p. 38

3.1.2 Convergence of the Gauss–Seidel iteration

  • The Gauss–Seidel analysis splits A into diagonal, lower, and upper parts and introduces the resulting iteration matrix.Source: §3.1.2 · p. 38

3.1. BASIC ITERATIVE METHODS FOR PDE DISCRETISATIONS

For PDE discretisations, Gauss–Seidel’s eigenvalues are related to Jacobi’s by squaring, while asymptotic rates show that Jacobi and Gauss–Seidel remain slow and even optimally relaxed SOR is unsuitable as a standalone solver for large systems.

3.1. BASIC ITERATIVE METHODS FOR PDE DISCRETISATIONS

  • If μ is an eigenvalue of Jacobi’s iteration matrix, μ^2 is an eigenvalue of Gauss–Seidel’s, yielding the stated conclusion that Gauss–Seidel converges twice as fast as Jacobi.Source: §3.1 · p. 39
  • For 1D Poisson with natural ordering, Jacobi is symmetric and diagonalisable, whereas Gauss–Seidel is non-normal and non-diagonalisable, with zero’s algebraic multiplicity ⌈n/2⌉ but geometric multiplicity one.Source: §3.1 · p. 39

3.1.3 Convergence of the SOR iteration

  • At n=49, kJ ≈7000, kGS ≈3500, and kωopt ≈110; optimal SOR uses roughly 64× fewer iterations than Jacobi, yet Thomas is cheaper for this tridiagonal system.The comparison targets a 10^-6 reduction in initial error; iterative sweeps cost O(n), while the Thomas algorithm costs O(n).Source: §3.1.3 · p. 40
  • Jacobi and Gauss–Seidel have O(1/n^2) rates, while optimally relaxed SOR improves to O(1/n); none is suitable as a standalone solver for large PDE systems.The notes contrast optimal SOR with modern solvers such as multigrid or Krylov methods.Source: §3.1.3 · p. 40

3.2 Basic Iterative Methods for Graph Laplacians · 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS

The section connects graph-Laplacian spectra with classical Jacobi iteration, interpreting it as diffusion on a graph and showing how slowly decaying modes can support community detection.

3.2 Basic Iterative Methods for Graph Laplacians

  • Every graph Laplacian has smallest eigenvalue 0 with the all-ones eigenvector, because each matrix row sums to 0.Source: §3.2 · pp. 40–41

3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS

  • For the path graph Pn, L is a rank-2 perturbation of the 1D Poisson matrix, and its eigenvalues are bounded using a symmetric-matrix perturbation result.The section states bounds for the kth and extreme eigenvalues but the supplied passages do not include their formulas.Source: §3.2 · p. 41
  • For the cycle graph Cn, Laplacian eigenvalues are λk = 2 − 2 cos(2πk/n), while Jacobi eigenvalues are cos(2πk/n).Source: §3.2 · p. 41
  • The star graph Sn has Laplacian spectrum {0, 1(×n − 2), n}.Source: §3.2 · p. 41
  • For P4, Jacobi encodes a walk that moves from each endpoint to its sole neighbor with probability 1 and splits equally between interior neighbors.Source: §3.2 · pp. 41–42
  • Jacobi iteration on a graph Laplacian evolves as a synchronous, or “lazy,” random walk whose eigenvalues describe diffusion rates.On the path example, the iterates spread like a probability distribution or heat diffusion.Source: §3.2 · p. 42
  • Slowly decaying Jacobi modes are Laplacian eigenvectors for small nonzero eigenvalues; thresholding the Fiedler vector can reveal graph communities [38] and bipartition the graph [10].These eigenvectors vary little within densely connected regions and more sharply across weak connections.Source: §3.2 · p. 42

3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS · 3.3. BASIC ITERATIVE METHODS IN MACHINE LEARNING

For graph Laplacians, Jacobi iterations act as diffusion and power iteration: they damp oscillatory graph modes while preserving smooth modes, linking basic iterative methods to Fiedler vectors and spectral clustering. The passages also characterize power-method convergence and note that these iterations are too slow to be efficient graph solvers.

3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS

  • For the example Laplacian, λ2 ≈ 0.438447 is the smallest nonzero eigenvalue, and thresholding its Fiedler vector at 0 yields two clusters.The vector changes sign across the bridge (3, 4).Source: §3.2 · p. 43
  • Jacobi with b = 0 is iteration by the random-walk matrix P; its spectral components show that high-frequency modes are attenuated while smooth modes persist.P is similar to the symmetric normalized adjacency S, which shares its eigenvalues.Source: §3.2 · p. 43
  • With a nonzero starting component along the dominant eigenvector, the power method converges to that eigenvector, with eigenvector error decaying like |λ2/λ1|^k.For symmetric matrices, eigenvalue error decays twice as fast; convergence is slow when the two largest eigenvalue moduli are close.Source: §3.2 · p. 44
  • After deflating the trivial eigenvector 1, Jacobi as power iteration on P converges to the Fiedler vector at rate |μ3/μ2|^k.The projection removes the dominant eigenpair (1, 1), using the D-inner product.Source: §3.2 · p. 44
  • Slow graph modes correspond to Laplacian eigenvectors with small eigenvalues, which capture global connectivity patterns and community structure.Jacobi’s enhancement of low-frequency information connects these iterations to spectral clustering.Source: §3.2 · p. 44
  • Basic graph iterations converge too slowly to serve as efficient graph solvers, but remain useful for understanding diffusion and spectral structure.Classical iterations act as diffusion or averaging processes on graphs.Source: §3.2 · p. 44

3.3 Basic Iterative Methods in Machine Learning

Classical iterative methods illuminate machine-learning optimization: gradient descent on quadratic objectives is Richardson iteration, and spectral properties govern convergence and learning speed. The section also connects Jacobi smoothing to spectral bias in neural-network training.

3.3.1 Richardson iteration and gradient descent

  • Gradient descent on a quadratic objective is exactly Richardson iteration, with the gradient ∇f(x) = Ax − b.This equivalence connects classical linear-system iterations directly to optimization.Source: §3.3.1 · p. 45
  • For SPD A, Richardson converges only when 0 < τ < 2/λmax; its rate depends on κ(A) = λmax/λmin, with small-eigenvalue modes decaying slowly.The slow modes are analogous to smooth modes that persist in PDEs and graph Laplacians.Source: §3.3.1 · p. 45
  • The minimax step balances the gains of the extreme eigenvalue modes, yielding q(τ⋆) = (κ − 1)/(κ + 1); for eigenvalues 1 and 9, τ⋆ = 0.2 and q⋆ = 0.8.A step based only on λmax can leave its mode at the divergence boundary while the λmin mode still decays.Source: §3.3.1 · p. 45
  • Jacobi is Richardson with preconditioner D−1, so preconditioning accelerates gradient descent by scaling its update direction.This interpretation carries over from classical iterations to machine-learning optimization.Source: §3.3.1 · p. 46

3.3.2 Applications in machine learning

  • In linear regression, the spectrum of X⊤X—the squared singular values of X—dictates convergence speed; ill-conditioning slows convergence and motivates preconditioning or SGD.The slow convergence mirrors the persistence of low frequencies in PDEs.Source: §3.3.2 · p. 46
  • For convex nonlinear losses, gradient descent is not exactly Richardson, but near the solution the Hessian spectrum governs learning: flat directions are slow and steep directions constrain stable step sizes.The local Hessian H = ∇2f(w∗) plays the role of A.Source: §3.3.2 · p. 46

3.3.3 Smoothing and spectral bias

  • Jacobi damps oscillatory high-frequency error while retaining smooth modes; analogously, neural-network training shows spectral bias, learning low-frequency/global patterns before fine details [31].In the Jacobi example, high-frequency modes are quickly damped, while smooth low-frequency modes decay slowly.Source: §3.3.3 · p. 46

3.3. BASIC ITERATIVE METHODS IN MACHINE LEARNING

Jacobi and gradient descent act as spectral filters, with the spectrum determining which modes persist. Their frequency roles differ between PDE solvers and learning, but the same spectral mechanism links iterative methods across both settings.

3.3. BASIC ITERATIVE METHODS IN MACHINE LEARNING

  • In learning, high-frequency patterns correspond to smaller eigenvalues, while low-frequency patterns are learned quickly, producing spectral bias.PDE solvers instead damp high-frequency error associated with large eigenvalues quickly.Source: §3.3 · p. 47
  • PDE solvers and learning use the same spectral filtering principle, although large eigenvalues correspond to high-frequency error in PDEs and low-frequency patterns in learning.The spectrum of A or X^⊤X dictates which modes persist, providing a framework spanning PDE solvers, optimization, and machine learning.Source: §3.3 · p. 47
  • Richardson iteration is gradient descent on quadratic problems, and Jacobi is Richardson iteration with a diagonal preconditioner.Source: §3.3 · p. 47
  • Directions corresponding to small eigenvalues decay slowly, creating convergence bottlenecks.The eigenvalue spectrum therefore determines convergence speed.Source: §3.3 · p. 47
  • In machine learning, the spectral account explains conditioning, the need for preconditioning or normalization, and spectral bias.Source: §3.3 · p. 47

3.4 Exercises

The exercises develop spectral and convergence analyses for classical iterations, then apply them to graph Laplacians, least squares, and power methods. They range from symbolic derivations to numerical verification and interpretation of error or spectral behavior.

Problem 1

  • Problem 1 asks for an h-independent Gauss–Seidel iteration matrix, its closed-form entries, and its plain-matrix pattern.The exercise directs students to compute B^-1 to obtain the entries.Source: §3.4 · p. 48

Problem 2

  • Problem 2 uses the Jacobi error eigenbasis to show that high-frequency error is damped quickly while smooth, low-frequency error persists.Source: §3.4 · p. 48

Problem 3

  • Problem 3 asks students to verify an eigenvector of the odd-sized tridiagonal Laplace matrix and show that all its eigenvalues lie in (0, 4).Source: §3.4 · p. 49

Problem 4

  • Problem 4 examines Jacobi convergence through the criterion of a nonzero diagonal and ρ(B) < 1, then asks for Toeplitz conditions and whether an iterate can solve a special case exactly.Source: §3.4 · p. 49

Problems 5–6

  • Problems 5 and 6 ask students to check predicted iteration convergence factors numerically and analyze the spectrum and Fiedler-vector bipartition of the five-vertex path graph.Source: §3.4 · pp. 49–50

Problem 7

  • Problem 7 uses the data spectrum of X^T X to ask about the condition number, admissible and optimal gradient-descent step sizes, convergence, and spectral bias.Source: §3.4 · p. 50

Problem 8

  • Problem 8 derives power-method eigenvalue convergence at rate |λ2/λ1|^2k, compares it with eigenvector convergence, and asks how projection affects convergence for a graph random-walk matrix.Source: §3.4 · p. 50

Krylov Methods for Symmetric Systems · Overview

This section develops Krylov methods for symmetric systems, connecting CG’s adaptive approximations to Lanczos projections for eigenvalue estimation and graph Fiedler vectors. It also examines early stopping and spectral filtering as regularization for noisy, ill-conditioned least-squares problems, alongside exercises on CG convergence and Lanczos properties.

Overview

  • Ridge regression, truncated SVD, and CG act as spectral filters of the same data, making early stopping of CG a form of regularization.For noisy data, a stopping rule remains necessary because of semi-convergence.Source: Overview · p. 51

4.1 The Conjugate Gradient method

  • For 1D Poisson systems, steepest descent takes on the order of n^2 iterations, while CG requires only n, motivating CG for large SPD PDE systems.Steepest descent's iteration count grows with the condition number; the Poisson stiffness matrix has condition number O(n^2).Source: §4.1 · pp. 53, 55

4.1. THE CONJUGATE GRADIENT METHOD

  • CG minimizes the A-norm error over each growing Krylov subspace and, in exact arithmetic, terminates within the number of distinct eigenvalues, at most n.In floating-point arithmetic, loss of orthogonality may delay or prevent exact termination.Source: §4.1 · pp. 54–55

4.2 The Lanczos Method and the Symmetric Eigenproblem

  • Lanczos supports computing a graph Laplacian's Fiedler vector, which encodes how the graph splits into communities [10].The chapter motivates iterative spectral methods because full eigendecomposition is too costly for large sparse matrices.Source: §4.2 · p. 57

4.2. THE LANCZOS METHOD AND THE SYMMETRIC EIGENPROBLEM

  • Lanczos compresses a large symmetric matrix into a tridiagonal projection using matrix-vector products; its Ritz values approximate eigenvalues, and lifted Ritz vectors approximate eigenvectors.Symmetry makes the projected matrix both symmetric and Hessenberg, hence tridiagonal; only scalar coefficients are needed at each step.Source: §4.2 · p. 57
  • Extremal Ritz values converge monotonically from inside the spectrum, interlace the eigenvalues, and can bracket the spectrum to estimate κ(A) for the CG bound.Interior eigenvalues converge more slowly and may require deflation.Source: §4.2 · p. 58

4.2.1 Lanczos and CG are the same recurrence

  • CG and Lanczos generate the same orthonormal Krylov basis in exact arithmetic, and CG scalars determine the tridiagonal projection entries.A CG run can estimate λmin, λmax, and κ(A) with a few scalar operations and no extra matrix-vector products; Lanczos exhausts the Krylov space after d steps, where d is the number of distinct eigenvalues.Source: §4.2.1 · pp. 59–60

4.2.2 Application: computing the Fiedler vector

  • After deflating the constant vector, Lanczos found the exact nonzero eigenvalues in three steps and recovered the Fiedler vector to machine precision, splitting the graph into its two triangles.The small entries ±0.2610 identify nodes 3 and 4 as the bottleneck across the cut.Source: §4.2.2 · p. 60

4.3.1 Least squares and the SVD

  • In noisy, ill-conditioned least squares, the solution divides noise by small singular values, so solving to machine precision can reproduce noise rather than the underlying model.When σi = 10^-3, a noise component of size 10^-3 shifts the corresponding coefficient by 1.Source: §4.3 · p. 61, §4.3.1 · p. 61

4.3.2 Spectral filters and the truncated SVD

  • Truncated SVD keeps the k best-determined directions and discards the rest, but requires a full SVD and a cutoff k chosen by hand.Source: §4.3.2 · p. 62

4.3.3 Ridge regression: a smooth filter

  • Ridge regression is a smooth spectral filter that damps small-singular-value directions; it avoids an SVD but requires choosing λ.CG can apply to X⊤X + λI using matrix-vector products without forming the d × d matrix.Source: §4.3.3 · p. 62

4.3.4 Early stopping: the filter generated by CG

  • Early-stopped CG generates a data-adapted spectral filter, resolves well-determined directions first, and minimizes the A-norm error among polynomial filters at every step.The iteration count k acts as a regularization parameter without requiring an SVD or a separately selected filter parameter.Source: §4.3 · p. 63

4.3. KRYLOV METHODS AND MACHINE LEARNING (OPTIONAL)

  • CG exhibits semi-convergence: unresolved solution components decrease with iteration while admitted noise grows, so error eventually rises and a stopping rule is indispensable.When noise level δ is known, stop when ∥y − Xx(k)∥2 ≤ τδ for τ > 1; otherwise, the passage recommends stopping when held-out validation error begins to rise.Source: Figure 4.5 · p. 62, §4.3 · p. 63

4.4 Exercises

  • The CG exercise uses Chebyshev optimality to derive the classical convergence estimate for a symmetric system with eigenvalue bounds m = λmin(A) and M = λmax(A).It defines κ(A) = M/m.Source: §4.4 · p. 65
  • In exact arithmetic, the Lanczos vectors q1, …, qm generated by the stated steps are to be shown orthonormal.The exercise assumes A = A⊤, ∥q1∥2 = 1, q0 = 0, and β0 = 0.Source: §4.4 · p. 65
  • For symmetric A, the exercise asks readers to show that the eigenvalues of the Lanczos tridiagonal matrix Tm interlace those of A.Here Tm = Qm⊤AQm, and θi denote its eigenvalues.Source: §4.4 · p. 65

Krylov Methods for Non-Symmetric Systems · Overview

The section develops Krylov methods for nonsymmetric systems, focusing on Arnoldi and GMRES, their convergence and residual properties, and applications to damped PageRank iteration and large least-squares problems. It also covers Golub–Kahan/LSQR connections to Krylov subspaces and CG on the normal equations, early stopping as a spectral filter, and related cost and per-step reduction exercises.

Overview

  • The chapter extends Krylov methods to nonsymmetric systems, using Arnoldi and GMRES for linear solves, Arnoldi for directed-network eigenproblems, and A and A⊤ for rectangular least squares.The overview identifies convection–diffusion, Navier–Stokes, transport and reaction networks, and the Google matrix as examples of nonsymmetric problems.Source: Overview · p. 67

5.1 GMRES and Arnoldi for non-normal matrices

  • Forming normal equations restores symmetry but squares the condition number and destroys sparsity, so the method instead works directly with A in K_k(A, r(0)).Source: §5.1 · p. 68
  • Arnoldi constructs the upper Hessenberg projection AV_k = V_{k+1}H_k; when A is symmetric, H_k is tridiagonal and Arnoldi reduces to Lanczos.Applying A to K_k stays within K_{k+1}, requiring one extra basis vector.Source: Figure 5.1 · p. 68, §5.1 · p. 68

5.1.1 GMRES: minimising the residual over the Krylov space

  • GMRES chooses the iterate with the smallest 2-norm residual over x(0) + K_k(A, r(0)), equivalently solving the reduced least-squares problem min_y ||βe1 − H_k y||2.The orthonormal Arnoldi basis preserves the residual norm in this reduced problem.Source: §5.1.1 · p. 69
  • GMRES residuals are orthogonal to A K_k, have the form q_k(A)r(0) with q_k(0) = 1, and terminate within the degree of the minimal polynomial relative to r(0).If h_{k+1,k} = 0, GMRES terminates in k steps.Source: Overview · p. 67, §5.1.1 · pp. 69–70

5.1.2 Convergence: why eigenvalues are not enough

  • For non-normal A, eigenvalues alone can mispredict GMRES convergence, whereas the field of values controls ||q(A)||2 and yields a rate when contained in a disk excluding zero.Even eigenvalues all equal to 1 can accompany stagnation for n − 1 steps when zero lies in the field of values.Source: §5.1.2 · pp. 71–72

5.2 Networks: PageRank as an Eigenproblem

  • For directed graphs, Arnoldi reduces eigenvalue computation to the small matrix H_m, while the Google matrix has a unique positive dominant eigenvector, the PageRank vector [24, 29], and can be applied using a sparse product plus vector addition.Source: §5.2 · pp. 72–73

5.2.1 The role of the damping factor

  • Every Google-matrix eigenvalue except λ1 = 1 satisfies |λ| ≤ α, so power iteration converges linearly with asymptotic factor at most α, independently of graph size.Source: §5.2.1 · p. 74
  • The damping factor makes PageRank scalable independently of graph size, but values near 1 improve link fidelity while slowing convergence.The iteration count behaves like log(1/ε)/(1 − α) as α approaches 1; on the three-node example, |λ2(G)| = α/√2, while on real web graphs the bound is essentially sharp.Source: §5.2.1 · p. 74

5.2.2 Beyond the power iteration: the Arnoldi process

  • Arnoldi can produce several Ritz pairs, tracks subspaces for nonnormal matrices, and supplies residuals for a posteriori stopping guarantees.Power iteration returns at most one eigenvector, has a fixed rate, and cannot certify its answer; Arnoldi's Ritz residual is |h_m+1,m||e⊤_m y|.Source: §5.2.2 · pp. 74–75

5.3 Krylov Methods for Large Least Squares (optional)

  • At 0.2% density, forming X⊤X produces a 7.7%-dense matrix, and forming the Gram matrix can erase rank information in floating point.For a 20,000 × 5,000 design matrix, X⊤X has nearly ten times as many nonzeros as X; with ε = 10^-9 in double precision, forming X⊤X annihilates the rank information.Source: §5.3 · pp. 76–77

5.3.1 Golub–Kahan bidiagonalisation and LSQR

  • Golub–Kahan bidiagonalisation alternates products with A and A⊤ to reduce least squares to a small problem, enabling LSQR to match CG-on-normal-equations iterates without numerical squaring.Short recurrences require O(m + n) work and storage and do not require retaining Vk; this makes LSQR stable where CGNR is not.Source: Figure 5.4 · p. 77, §5.3.1 · p. 78

5.3.2 Early stopping and the truncated SVD

  • Early stopping in LSQR acts as a spectral filter, suppressing small-singular-value directions that amplify noise; its iteration count serves as a regularisation parameter.For latent semantic analysis, a small SVD of the bidiagonal matrix approximates leading singular directions and yields a low-dimensional semantic space from a few dozen matrix-vector products.Source: Figure 5.5 · p. 78, §5.3.2 · p. 78

5.4 Exercises

  • The exercises cover Krylov spaces and GMRES convergence, residual polynomial and finite-termination properties, field-of-values bounds, Arnoldi Ritz residuals, and a PageRank example.Source: §5.4 · pp. 80–81
  • Problem 7 asks students to derive Golub–Kahan bidiagonalization from the recurrence, including the factorization AV_k = U_{k+1}B_k.It also asks them to write B_3 explicitly.Source: §5.4 · p. 82
  • The exercises ask students to show that the generated right and left vector spans are Krylov subspaces and establish the sense in which LSQR is CG on the normal equations.Source: §5.4 · p. 82
  • The exercises examine reducing each LSQR step to a bidiagonal least-squares problem and using a Givens rotation, contrasting its per-step work and storage with GMRES orthogonalization.Source: §5.4 · p. 82

Preconditioning and Accelerating Solvers · Overview

This section develops preconditioning to accelerate solvers, covering basic and SPD preconditioning, a subspace framework, additive Schwarz, and multigrid. It examines how overlap, local solves, smoothing, and coarse-grid correction affect convergence and solver efficiency.

Overview

  • Preconditioning targets mesh-dependent Krylov convergence by using a cheaper-to-invert approximation M to cluster the transformed spectrum and bound its condition number independently of mesh size.For PDE matrices, the unpreconditioned iteration count grows as the mesh is refined.Source: Overview · p. 83

6.1 Motivation and Basic Preconditioners

  • For SPD problems, a good preconditioner replaces A with an operator whose eigenvalues are more tightly clustered, ideally near 1.Preconditioning has left, right, and symmetric formulations; symmetric preconditioning preserves symmetry and positive definiteness.Source: §6.1 · p. 84

6.1.1 Preconditioning SPD problems

  • For SPD systems, preconditioned CG requires A and M to be SPD and uses solves with M in the original variables rather than requiring M^-1/2.The preconditioner preserves the exact solution while changing the Krylov method’s residual-minimization trajectory.Source: Overview · p. 83, §6.1 · p. 85

6.1.2 A general framework for preconditioners

  • For the 1D discrete Laplacian, Jacobi only rescales eigenvalues, leaving κ(Ah) ∼ h^-2 and unable to improve conditioning as the mesh is refined.Jacobi can help with coordinate scaling, but the uniform diagonal here means the difficulty is not local scaling.Source: §6.1.2 · p. 86
  • Incomplete factorization approximates A by sparse factors, discarding fill outside a prescribed sparsity pattern.Source: §6.1.2 · p. 86

6.1.3 One framework for the whole chapter

  • The chapter’s preconditioners share a subspace construction: independent local solves correct separate subspaces, and adding a global subspace addresses what local solves miss.Each local solve is cheap to apply and the independent solves can be parallelized.Source: §6.1.3 · p. 86

6.2 Schwarz Domain Decomposition

  • With overlap fixed as a fraction of the domain, discrete Schwarz convergence is geometric at a rate independent of h; increasing overlap improves the contraction factor.In the 1D model, vanishing overlap drives the factor toward 1 and stalls convergence.Source: §6.2 · p. 88

6.2. SCHWARZ DOMAIN DECOMPOSITION

  • Additive Schwarz restricts the residual to overlapping subdomains, solves local correction problems, prolongates and blends the corrections, and sums them from the same residual.Unlike additive Schwarz, alternating Schwarz is multiplicative because each subdomain uses updated information from the preceding solve.Source: §6.2 · p. 89
  • Additive Schwarz combines local subdomain solves and blends their corrections across the overlap using a discrete partition of unity.The restriction operators extract local unknowns, their transposes prolongate corrections, and the partition weights preserve a global approximation of A^-1.Source: §6.2 · p. 90

6.3 Two-Level Preconditioning and Multigrid Ideas

  • Jacobi is a smoother rather than a standalone solver: it damps oscillatory modes quickly but leaves smooth, slowly varying error comparatively untouched.The discrete sine modes distinguish these scales, and the Jacobi error propagator damps large-index modes faster.Source: §6.3 · pp. 90–91

6.3.1 Weighted Jacobi and the smoothing factor

  • ω = 2/3 minimizes weighted Jacobi’s worst-case high-frequency damping factor, reducing every such mode by at least 1/3 per step.Weighted Jacobi with this parameter is an effective smoother, not a scalable solver; it leaves smooth components for coarse correction.Source: §6.3.1 · p. 92

6.3.2 Coarse-grid correction

  • Coarse-grid correction represents smooth error using values at half the grid points and removes the error visible to the coarse space, but cannot converge alone.It annihilates error in its range while leaving error outside that range untouched; the smoother handles the complementary oscillatory component.Source: §6.3.2 · p. 93

6.3.3 From coarse correction to the two-grid method

  • The two-grid method pre-smooths, restricts the residual, solves and corrects on the coarse grid, then post-smooths.Its coarse solve can be replaced recursively by further two-grid steps until the grid is small enough for a direct solve.Source: §6.3.3 · p. 94
  • O(n) work and h-independent error reduction per V-cycle make recursive multigrid optimal; in 1D, weighted Jacobi reduces high-frequency modes by a factor of 1/3 per step.For elliptic problems, multigrid-preconditioned Krylov iteration counts remain essentially bounded under mesh refinement.Source: §6.3.3 · pp. 94–95

6.4 Exercises

  • Using weighted Jacobi S and the two-grid operator M_TG = TS^ν, the exercises ask why ω = 1 is a poor smoother despite minimizing ρ(S).Source: §6.4 · p. 97
  • The exercises ask readers to compute ρ(M_TG) for ν = 1 and 2, compare it with the smoothing factor, and state how it relates to ν.Source: §6.4 · p. 97
  • They ask whether the two-grid contraction factor is independent of h for m = 15 and 31, contrasting it with Jacobi’s O(h^-2) iteration count.Source: §6.4 · p. 97
  • The exercises ask readers to show that a V-cycle costs O(m) across successively halved grids and explain why this, together with h-independent contraction, makes multigrid optimal.Source: §6.4 · p. 97

Related papers

Loading 2608.21234v1…