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

These notes address structured linear-algebra problems arising in PDEs, networks, data assimilation, and noisy machine learning that are too large to factorise. They develop classical and Krylov-based algorithms and show that early-stopped CG acts as spectral-filter regularisation, while iterative methods support applications such as spectral clustering.

  • Problem

    Large structured systems in PDEs, networks, data assimilation, and noisy learning are often too large to factorise and accessible mainly through matrix-vector products.

  • Method

    The course develops classical factorizations, stationary iterations, Krylov methods, preconditioning, domain decomposition, and multigrid, applying them across these settings.

  • Results

    The notes show that ridge regression, truncated SVD, and CG are spectral filters of the same data, with early CG stopping providing regularisation and Fiedler-vector thresholding enabling spectral clustering.

  • Takeaways & Limitations

    A small set of numerical linear-algebra ideas transfers across sparse systems, graph ranking, data assimilation, spectral clustering, and large noisy least-squares problems.

  • Takeaways & Limitations

    In double precision, forming X^T X can destroy rank information when ε = 10^-9, making the Gram matrix exactly singular.

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

The chapter establishes a unified numerical foundation for solving Ax = b, covering norms, spectral decompositions, floating-point error, conditioning, stability, direct factorizations, and stationary iterations. It contrasts robust but poorly scaling direct methods with sparse, scalable iterative approaches and gives their convergence conditions and roles in later algorithms.

  • Overview: The chapter develops norms, conditioning, floating-point stability, LU and Cholesky factorizations, and stationary iterations as core tools for numerical linear algebra.These foundations support applications including structural analysis, machine learning, weather forecasting, and data assimilation.
  • 1.1 Norms and Matrix Decompositions · 1.1.1 Vector and Matrix Norms · 1.1.3 Matrix Norms: Vector and matrix norms measure size, distances, and approximation errors; all norms on finite-dimensional spaces are equivalent, although their numerical values can differ.For example, ∥x∥∞≤∥x∥2 ≤√n∥x∥∞ and ∥x∥2 ≤∥x∥1 ≤√n∥x∥2.
  • 1.1.2 Eigenvalues and Matrix Decompositions: Eigenvalues and the spectral radius govern iterative behavior, while the Schur decomposition and SVD exist for every matrix and underpin general-purpose numerical algorithms.The SVD represents A using unitary U and V and nonnegative singular values.
  • 1.2.1 Floating-Point Arithmetic · Sources of Numerical Errors: Floating-point representation introduces rounding, truncation, cancellation, overflow, and non-associativity, with lower precisions more exposed to overflow and rounding errors.Catastrophic cancellation can make subtracting nearly equal stored numbers inaccurate or even zero.
  • 1.2.2 Error, Conditioning and Stability of algorithms: Conditioning describes a problem’s sensitivity to perturbations, while stability describes error propagation; together they determine final accuracy.A system close to singularity can be highly sensitive to small perturbations in A or b, even when the algorithm is stable.
  • 1.2.3 Gaussian Elimination and LU Factorization · 1.3 Stationary Iterative Methods: Gaussian elimination computes A = LU in O(n3) flops, with O(n2) forward and backward substitutions, so its scaling and fill-in motivate iterative solvers for large sparse systems.For n ∼106, O(n3) is infeasible; partial pivoting yields PA = LU and backward stability.
  • 1.2.4 Cholesky Factorization for SPD Matrices: For SPD matrices, Cholesky A = LL⊤ requires no pivoting, is backward stable, halves LU storage, and costs ≈n3/3 flops.If a nonpositive diagonal factor occurs, the notes recommend LDL⊤ with symmetric pivoting or regularization via A ←A + λI.
  • 1.3 Stationary Iterative Methods · 1.3.1 Examples of Stationary Methods and their convergence: Stationary iterations update solutions through A = M −N and converge for every initial guess exactly when ρ(G) < 1; Jacobi, Gauss–Seidel, SOR, and Richardson have stated matrix-dependent conditions.Jacobi and Gauss–Seidel converge under strict row diagonal dominance, Gauss–Seidel converges for SPD matrices, SOR converges for 0 < ω < 2, and these methods remain useful as smoothers and preconditioners.

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

The chapter explains how finite-difference discretization produces sparse linear systems from PDEs, then connects their structure to graph, machine-learning, and data-assimilation matrices. It also shows how storage, ordering, fill-in, and spectral properties determine computational cost and iterative-solver convergence.

  • 2.1 Finite difference approximation: The three-point central-difference stencil is consistent with u′′(x) and has truncation error O(h2).Central differences are second order accurate for u′ and u′′, whereas forward/backward differences are first order for u′.
  • 2.1 Finite difference approximation: Finite-difference discretization converts differential operators into sparse matrices whose nonzero patterns are determined by local stencils.The 1D Poisson stencil is [−1, 2, −1]/h2, and its bandwidth controls direct-method cost.
  • 2.2.1 The 2D Poisson equation: structure, stencils and blocks: The 2D Poisson discretization yields a block-tridiagonal matrix with tridiagonal blocks, at most 5 nonzeros per row, symmetry, and positive definiteness under Dirichlet data.Lexicographic ordering couples each grid point to its same-row neighbors and neighboring grid rows.
  • 2.2.2 Graph Laplacians: Graph Laplacians L = D −A are sparse, symmetric, positive semidefinite matrices with nnz(L) = 2|E| + n.They extend the Laplace operator from continuous domains to arbitrary discrete networks and support spectral clustering, graph neural networks, and data assimilation.
  • 2.2.3 Sparse Matrices in Machine Learning: Machine-learning sparsity arises from locality, so normal equations, Hessians, and Jacobians inherit sparse structure that is essential for large-scale computation.If each feature row has k ≪n nonzeros, then nnz(A) = O(mk); for logistic regression, H = A⊤WA is sparse when A is sparse.
  • 2.2.4 Sparsity in data assimilation: Data-assimilation covariance estimates can create spurious long-range correlations or dense long-tailed matrices, so localisation and thresholding impose local sparsity or fixed bandwidth.Correlation distributions may be truncated below a typically chosen threshold of 0.2.
  • 2.3 Sparse Matrix Toolbox; 2.3.1 Sparsity Patterns and storage formats: Sparse computation requires suitable storage formats and attention to bandwidth, profile, reordering, and fill-in.CSR/CSC provide general-purpose compressed storage, while diagonal/banded formats suit near-diagonal nonzeros; effective reordering can save orders of magnitude in memory and runtime.
  • 2.3.1 Sparsity Patterns and storage formats: For 2D Poisson with N = n2, naive banded LU costs O(N2), whereas nested dissection reduces the cost to O(N 3/2).Although each row has only 5 nonzeros, the bandwidth grows like n; fill depends strongly on the sparsity pattern.

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

The chapter applies fixed-point iterations to one-dimensional Poisson discretisations, where explicit spectral analysis yields comparable convergence rates. It connects these mechanisms to smoothing, graphs, and machine learning applications.

  • Overview: Explicit eigenpairs for the 1D Poisson matrix make the convergence rates of Jacobi, Gauss–Seidel, and SOR exactly computable.The matrix is symmetric, positive definite, and tridiagonal, with κ(A) = O(n^2).
  • Overview: Gauss–Seidel converges roughly twice as fast as Jacobi, while optimally weighted SOR improves the rate further.These comparisons follow from the explicit spectral-radius analysis of the model problem.
  • Applications of Basic Iterative Methods: Jacobi iteration on graph Laplacians can be interpreted as random-walk or diffusion dynamics linked to the Fiedler vector and spectral clustering.Its slow modes provide the connection between iterative convergence and graph structure.
  • Applications of Basic Iterative Methods: Richardson iteration is gradient descent on a quadratic loss, with the spectrum of A or X^⊤X determining step size and convergence rate.The same spectral viewpoint links basic iterations to machine learning optimization.
  • Applications of Basic Iterative Methods: Spectral bias in machine learning is described as the same spectral-filtering mechanism as PDE smoothing, with low- and high-frequency roles reversed.The comparison transfers the iterative-method perspective across application domains.
  • 3.1 Basic Iterative Methods for PDE discretisations: The Poisson model is simple enough for explicit eigenvalue analysis yet representative of phenomena that recur in higher-dimensional PDEs, graphs, and machine learning.The discretisation uses n interior points and homogeneous Dirichlet boundary conditions.
  • 3.1 Basic Iterative Methods for PDE discretisations: High-frequency Poisson modes are damped rapidly by Jacobi, Gauss–Seidel, and SOR, whereas low-frequency modes decay slowly.This smoothing property is later used in multigrid methods.
  • 3.1.1 Convergence of the Jacobi iteration: Jacobi iteration becomes very slow as the grid grows because its spectral radius approaches one.This behavior is the direct consequence of the Jacobi eigenvalue analysis for the Poisson matrix.

3.1. BASIC ITERATIVE METHODS FOR PDE DISCRETISATIONS

The section relates Gauss–Seidel and Jacobi spectra, shows Gauss–Seidel’s non-diagonalisability for the one-dimensional Poisson problem, and compares the asymptotic convergence of Jacobi, Gauss–Seidel, and optimally relaxed SOR. Although SOR improves the rate, none of these methods is suitable as a standalone solver for large PDE systems.

  • Jacobi and Gauss–Seidel: If μ is an eigenvalue of Jacobi, then μ^2 is an eigenvalue of Gauss–Seidel, so Gauss–Seidel converges twice as fast as Jacobi.This spectral relation underpins the comparison of the two methods’ convergence factors.
  • Jacobi and Gauss–Seidel: For the one-dimensional Poisson problem, Gauss–Seidel has only ⌊n/2⌋ distinct nonzero simple eigenvalues, while 0 has algebraic multiplicity ⌈n/2⌉ and geometric multiplicity one.The unequal algebraic and geometric multiplicities imply Jordan blocks of size greater than one, making the iteration matrix non-diagonalisable and non-normal.
  • Asymptotic convergence: As n grows, Jacobi and Gauss–Seidel retain O(1/n^2) convergence rates, whereas optimally relaxed SOR improves to O(1/n) per iteration.Jacobi’s spectral radius approaches 1, Gauss–Seidel roughly doubles the effective rate, and optimal SOR remains poor relative to multigrid or Krylov methods.
  • Numerical comparison: 64× fewer iterations are required by optimally relaxed SOR than Jacobi for n = 49 to reduce the initial error by 10^-6.The corresponding iteration counts are kJ ≈7000, kGS ≈3500, and kωopt ≈110.
  • Limitations: None of Jacobi, Gauss–Seidel, or optimally relaxed SOR is suitable as a standalone solver for large PDE systems.For the one-dimensional tridiagonal example, direct Thomas factorisation is cheaper, while iterative methods become more valuable in two and three dimensions.

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

This section characterizes graph-Laplacian spectra and uses Jacobi iteration as a microscope for graph structure. It connects the iteration to random-walk diffusion and explains how slowly decaying modes support community detection.

  • 3.2 Basic Iterative Methods for Graph Laplacians: Every graph Laplacian has smallest eigenvalue 0 with the all-ones vector as a corresponding eigenvector.Each Laplacian row sums to 0.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: For path graphs, the Laplacian is a rank−2 perturbation of the one-dimensional Poisson matrix, enabling eigenvalue bounds.The bounds follow by writing L = ALap + D and applying the stated symmetric-matrix eigenvalue result.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: For path graphs, the extreme Laplacian eigenvalues receive explicit bounds from the perturbation-based corollary.The section separately records the resulting bounds for the extreme eigenvalues.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Cycle graphs have Fourier-mode eigenvectors with Laplacian eigenvalues λ_k = 2 −2 cos(2πk/n) and Jacobi eigenvalues cos(2πk/n).The star graph instead has spectrum {0, 1(×n −2), n}.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Classical iterations remain useful for interpreting graph structure even though practical graph algorithms typically use more advanced methods.Their role is primarily explanatory, connecting iterative behavior with classical graph-theoretic interpretations.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: On the four-vertex path, Jacobi gives endpoint transitions with probability 1 and interior transitions to either neighbor with probability 1/2.This follows from the transition-probability interpretation of the Jacobi iteration matrix.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Jacobi iteration on a graph Laplacian is equivalent to synchronous random-walk evolution, with eigenvalues describing diffusion rates.The updates spread mass along the path like a probability distribution or heat diffusion, also admitting a lazy-random-walk interpretation.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: The slowest-decaying Jacobi modes are Laplacian eigenvectors associated with the smallest nonzero eigenvalues and are smooth within dense regions but change across weak connections.Spectral graph theory uses this structure, including the Fiedler vector, to detect communities.

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

For graph Laplacians, Jacobi iteration is both a power method and a low-pass graph filter: it damps high-frequency oscillations while preserving smooth modes linked to connectivity and community structure. Deflating the trivial random-walk mode makes Jacobi converge toward the Fiedler vector, although convergence is generally too slow for efficient graph solving.

  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: The graph Laplacian has λ2 ≈0.438447 as its smallest nonzero eigenvalue, and the associated Fiedler vector changes sign across bridge (3, 4).Thresholding the vector at zero or using a sweep cut yields the two graph clusters.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Jacobi iteration on the random-walk matrix is a low-pass graph filter: high-frequency eigencomponents are attenuated rapidly, while smooth modes persist.The spectral multiplier is |µ_i|^k = |1 − λ_i^sym|^k.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: With a nonzero initial component along the dominant eigenvector, power iteration converges up to sign, with eigenvector error decaying like |λ2/λ1|^k.For symmetric matrices, the eigenvalue error decays twice as fast; convergence is slow when the two largest eigenvalue moduli are close.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: For the random-walk operator P, projecting out the trivial eigenvector 1 and normalizing in the D-inner product yields convergence to the Fiedler vector at rate |µ3/µ2|^k.Thus Jacobi is power iteration on P with deflation.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Although Jacobi iterations are conceptually valuable for understanding graph diffusion and spectral structure, their convergence is too slow for efficient graph solvers.The framework connects the random-walk operator, its symmetric diagonalization, power iteration, and the Fiedler vector’s role in spectral clustering.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Jacobi’s eigenmode analysis shows that high-frequency graph oscillations are quickly damped, whereas smooth global modes remain.These iterations act as diffusion or averaging processes on graphs.
  • 3.2. BASIC ITERATIVE METHODS FOR GRAPH LAPLACIANS: Laplacian eigenvectors with small eigenvalues represent slow modes that capture global connectivity patterns and reveal community structure.Jacobi therefore enhances the low-frequency information underlying spectral graph partitioning.

3.3 Basic Iterative Methods in Machine Learning

The section reframes classical iterative methods as optimization tools: Richardson iteration is gradient descent for quadratic objectives, while spectral properties govern convergence, stability, and learning dynamics. These ideas extend to regression, nonlinear losses, preconditioning, and neural-network spectral bias.

  • Richardson and Gradient Descent: Gradient descent on quadratic objectives is exactly Richardson iteration because the gradient of the quadratic loss is Ax − b.This identifies a classical linear-system iteration as a machine-learning optimization method.
  • Spectral Behaviour: For symmetric positive definite A, convergence requires 0 < τ < 2/λmax, while the condition number controls the rate and slows modes associated with eigenvalues near zero.The slow-mode bottleneck parallels smooth modes in PDEs and graph Laplacians.
  • Step-Size Selection: The minimax-optimal Richardson step balances the extreme spectral modes, giving q(τ⋆) = κ^-1 rather than optimizing against the largest eigenvalue alone.The two-eigenvalue example with λmin = 1 and λmax = 9 uses τ⋆ = 0.2 and equal mode gains of 0.8.
  • Preconditioning and Regression: Jacobi is Richardson with preconditioner D^-1, and preconditioning accelerates gradient descent by scaling the update direction.In linear regression, ill-conditioned X⊤X causes slow convergence, motivating preconditioning and stochastic gradient descent.
  • Nonlinear Losses and Spectral Bias: For nonlinear convex losses, the Hessian governs local convergence: small eigenvalues create slow learning, while large eigenvalues require small stable steps.In neural networks, gradient descent exhibits spectral bias, learning low-frequency global patterns before fine details.

3.3. BASIC ITERATIVE METHODS IN MACHINE LEARNING

Basic iterative methods in machine learning act as spectral filters, just as they do in PDE solvers, although low- and high-frequency roles are reversed. Their convergence depends on the eigenvalue spectrum, explaining slow directions, conditioning, preconditioning, normalization, and spectral bias.

  • Spectral filtering: Iterative methods filter spectral modes, with the spectrum of A or X⊤X determining which modes persist.In PDE solvers, large eigenvalues correspond to high-frequency error damped quickly; in learning, they correspond to low-frequency patterns learned quickly.
  • Convergence and conditioning: Convergence slows in directions associated with small eigenvalues, creating spectral bottlenecks.The eigenvalue spectrum therefore controls convergence speed.
  • Convergence and conditioning: In machine learning, spectral behavior explains conditioning, the need for preconditioning and normalization, and spectral bias.High-frequency patterns correspond to smaller eigenvalues, producing spectral bias.
  • Unified framework: Classical iterative methods provide a common spectral framework for PDE solvers, optimization, and machine learning.Richardson iteration is gradient descent on quadratic problems, while Jacobi is Richardson with a diagonal preconditioner.

3.4 Exercises

The exercises apply iterative-method theory to explicit matrices, convergence, spectra, and numerical experiments. They also extend these ideas to graph Laplacians, least-squares gradient descent, and the power method.

  • Problems 1–2: Problems 1–2 derive the 1D Gauss–Seidel iteration matrix and analyze Jacobi error components by frequency.The exercises ask for a closed form and matrix pattern for Gauss–Seidel, then show that Jacobi quickly damps high-frequency errors while low-frequency errors persist.
  • Problem 5: Problem 5 verifies predicted convergence rates for Jacobi, Gauss–Seidel, and optimally weighted SOR on the 1D Poisson matrix with n = 19.Students compute spectral radii and iteration counts, implement the methods, and compare observed asymptotic factors with predictions.
  • Problem 6: Problem 6 computes the spectrum of the five-vertex path graph and identifies the Fiedler vector and its induced bipartition.The computed eigenpairs are checked against λk(L) = 2 − 2 cos(πk/n).
  • Problem 7: Problem 7 uses singular values σ1 = 4, σ2 = 2, and σ3 = 1 to study least-squares conditioning, step sizes, convergence, and spectral bias.It asks for κ(X⊤X), admissible and minimax-optimal step sizes, iterations for a factor of 100, and the fastest and slowest learned modes.
  • Problem 8: Problem 8 derives power-method convergence, including eigenvalue rate |λ2/λ1|2k, a 1D Poisson example, and convergence to the Fiedler vector after projection.The exercises compare eigenvector and eigenvalue convergence, estimate iterations for n = 19, and explain the effect of removing the dominant all-ones eigenvector.

Krylov Methods for Symmetric Systems · Overview

The chapter replaces slow fixed-splitting iterations with adaptive Krylov methods built from successively richer matrix-generated subspaces. It develops CG and Lanczos for symmetric systems, then connects them to spectral estimation, graph clustering, and regularised least squares.

  • Overview: Krylov methods replace fixed-factor stationary iterations with adaptive approximation spaces generated by r(0), Ar(0), and higher matrix powers.For PDE matrices, stationary convergence factors behave like 1 − O(h2), making refinement especially slow.
  • 4.3 Krylov Methods and Machine Learning (optional): Ridge regression, truncated SVD, and CG act as spectral filters, while early CG stopping regularises noisy least-squares solutions without SVD or an additional parameter.CG resolves well-determined directions first, adapts its polynomial filter to the data, and requires a stopping rule because semi-convergence eventually increases admitted noise.
  • 4.1 The Conjugate Gradient method: For the 1D Poisson matrix, steepest descent needs on the order of n2 iterations, whereas CG requires only n iterations.The broader comparison gives stationary PDE methods O(n2) scaling and CG O(√κ(A)) iteration requirements.
  • 4.1 The Conjugate Gradient method: For SPD systems, CG selects A-conjugate search directions and minimizes the A-norm error over x(0) + Kk(A, r(0)).In exact arithmetic, A-conjugacy eliminates error components independently, while the iterate is the A-orthogonal projection onto the affine Krylov space.
  • 4.1 The Conjugate Gradient method: At most d ≤ n CG steps suffice in exact arithmetic, where d is the degree of A’s minimal polynomial and, for symmetric A, the number of distinct eigenvalues.Floating-point loss of orthogonality may delay or prevent exact termination.
  • 4.2 The Lanczos Method and the Symmetric Eigenproblem: Lanczos approximates eigenvalues and eigenvectors through Ritz values and vectors of a tridiagonal projection, with extremal values converging before interior ones.A few steps can bracket λmin and λmax, providing the spectral information needed for κ(A) estimates and CG bounds.
  • 4.2.1 Lanczos and CG are the same recurrence: CG and Lanczos generate the same orthonormal Krylov basis, so a CG solve can estimate λmin, λmax, and κ(A) without extra matrix–vector products.Lanczos exhausts the Krylov space after d steps, explaining CG’s finite termination from the spectral viewpoint.
  • 4.2.2 Application: computing the Fiedler vector: Deflating the constant vector makes the Fiedler eigenvalue extremal, enabling Lanczos to recover graph community structure through thresholding.In the example, three steps produced the exact nonzero eigenvalues and machine-precision Fiedler vector, separating {1, 2, 3} from {4, 5, 6}.

Krylov Methods for Non-Symmetric Systems · Overview

The chapter develops Krylov methods for nonsymmetric and often non-normal systems, replacing Lanczos with Arnoldi and using GMRES to minimize residuals directly. It extends the framework to PageRank, eigenvalue computation, and large least-squares problems through scalable matrix-vector methods, damping bounds, and LSQR.

  • Overview: Arnoldi builds an orthonormal Krylov basis with an upper Hessenberg projection, reducing to Lanczos when A is symmetric.The decomposition is AV_k = V_k+1H_k, and applying A to K_k produces at most one extra basis vector.
  • 5.1 GMRES and Arnoldi for non-normal matrices: GMRES chooses the Krylov-space iterate with the smallest 2-norm residual, equivalently solving min_y ∥βe1 − H_ky∥2.Its residual is orthogonal to AK_k, decreases monotonically, and terminates in at most d steps when d is the relative minimal-polynomial degree.
  • 5.1.2 Convergence: why eigenvalues are not enough: For non-normal A, eigenvalues alone can mispredict GMRES convergence because ill-conditioned eigenvectors make spectral bounds vacuous.The field of values controls ∥q(A)∥2 directly; if it lies in a disk excluding zero, the polynomial q_k(z) = (1 − z/c)^k gives a factor bounded by (s/|c|)^k.
  • 5.2 Networks: PageRank as an Eigenproblem: PageRank is the positive, unit-1-norm dominant eigenvector of a stochastic irreducible Google matrix, which can be applied without storing its dense rank-one perturbation.The Google matrix represents directed-web ranking and is handled through matrix-vector operations rather than full eigendecomposition.
  • 5.2.1 The role of the damping factor: |λ2(G)| ≤ α bounds power-iteration convergence independently of graph size, while larger α improves link fidelity but slows iteration.Arnoldi complements power iteration by producing several Ritz pairs and rigorous a posteriori residuals, including |h_m+1,m||eᵀ_m y|.
  • 5.3 Krylov Methods for Large Least Squares (optional): Large least-squares problems should avoid normal equations because forming X⊤X destroys sparsity and can lose rank information in floating-point arithmetic.For a 20 000 × 5000 matrix at 0.2% density, X⊤X becomes 7.7% dense; in double precision, rank information is already annihilated at ε = 10^-9.
  • 5.3.1 Golub–Kahan bidiagonalisation and LSQR: LSQR mirrors GMRES by minimizing ∥γe1 − B_ky∥2, using Golub–Kahan bidiagonalization and short recurrences with O(m + n) work and storage.The bidiagonal matrix is triangularized by one Givens rotation per step, avoiding the need to retain V_k.
  • 5.3.2 Early stopping and the truncated SVD: Early stopping obtains useful regularized least-squares solutions from a few dozen matrix-vector products instead of dense factorization.This connects Krylov computation with truncated-SVD-style filtering for large noisy data problems.

Preconditioning and Accelerating Solvers · Overview

Preconditioning accelerates Krylov methods by replacing a poorly conditioned operator with one whose spectrum is clustered, while Schwarz and multigrid extend this idea through local and coarse corrections. The chapter develops these constructions from basic preconditioners to mesh-independent multilevel solvers.

  • Overview: Preconditioning replaces A with M^-1A, where M approximates A and is cheap to invert, to cluster eigenvalues and bound conditioning independently of mesh size.For PDE matrices, unpreconditioned CG may require O(h^-1) iterations as the mesh is refined.
  • 6.1.2 A general framework for preconditioners: Jacobi preconditioning is cheap and parallel but cannot improve the 1D discrete Laplacian, whose condition number remains κ(A_h) ∼ h^-2 under mesh refinement.Jacobi rescales every eigenvalue by the same constant because the Laplacian diagonal is uniform.
  • 6.1.3 One framework for the whole chapter: The chapter unifies its preconditioners as sums of independent local solves on selected subspaces, with a global subspace added to capture errors missed by purely local corrections.This construction supports parallel application and organizes the progression from local methods toward coarse-space corrections.
  • 6.2 Schwarz Domain Decomposition: Schwarz methods replace pointwise updates with overlapping subdomain solves, and their 1D contraction improves as overlap widens while remaining independent of h when overlap is a fixed domain fraction.Additive Schwarz restricts residuals, solves local problems, prolongates corrections, and blends overlap contributions; its local solves are independent and suit Krylov preconditioning.
  • 6.3.1 Weighted Jacobi and the smoothing factor: Weighted Jacobi with ω = 2/3 reduces every high-frequency mode by at least a factor 1/3 per step, but leaves smooth error for coarse-grid correction.Thus relaxation is a smoother rather than a scalable solver, and effective multigrid combines local smoothing with global coarse correction.
  • 6.3.2 Coarse-grid correction: Smoothing and coarse correction are complementary: the smoother cannot remove smooth error, while coarse correction removes only error components represented by the coarse space.Each ingredient is ineffective alone, motivating their combination in two-grid and multigrid cycles.
  • 6.3.3 From coarse correction to the two-grid method: The multigrid V-cycle has O(n) work, error reduction independent of h, and iteration counts that remain essentially bounded under mesh refinement.Its reduction is governed by the smoothing factor rather than κ(A), exemplified by 1/3 per weighted-Jacobi step in 1D.
Loading 2608.21234v1…