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

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

• The paper presents numerical linear algebra as a common framework for solving large, structured linear algebra problems arising from PDE discretization.
• Key applications include graph analysis, machine learning, and data assimilation, each utilizing ideas like sparsity, spectral analysis, conditioning, Krylov subspaces, and preconditioning.

Scope and organizing thesis

“Advanced Linear Algebra with Applications – Part I” is a set of master’s-level lecture notes by Victorita Dolean and Jemima M. Tabeart that presents numerical linear algebra as a common computational framework for PDE discretization, graph analysis, machine learning, and data assimilation (2608.21234). Its central thesis is that these applications repeatedly generate large, structured linear-algebra problems whose practical treatment depends on a small collection of ideas: sparsity, spectral analysis, conditioning, Krylov subspaces, and preconditioning.

The notes are organized progressively. They begin with norms, matrix decompositions, floating-point arithmetic, conditioning, direct factorization, and stationary iterations. They then derive sparse matrices from finite-difference PDE discretizations and connect these structures to graph Laplacians, machine-learning design matrices, Hessians, and covariance models. Subsequent chapters analyze Jacobi, Gauss–Seidel, SOR, and Richardson spectrally; introduce CG and Lanczos for symmetric systems; develop Arnoldi and GMRES for nonsymmetric systems; discuss LSQR for rectangular least-squares problems; and conclude with Schwarz methods, coarse-grid correction, and multigrid.

The pedagogical strategy is explicitly application-oriented. Classical algorithms are not treated merely as historical techniques: Jacobi is interpreted as a diffusion operator, Richardson as gradient descent, Lanczos as a tool for Fiedler-vector computation, and early-stopped CG or LSQR as a regularization mechanism. Reproducible Python experiments are intended to complement the mathematical exposition.

Numerical foundations: norms, spectra, conditioning, and stability

The first chapter establishes the standard hierarchy of numerical linear algebra concepts. Vector norms $\ell^1$, $\ell^2$, and $\ell^\infty$ are introduced together with norm equivalence in finite-dimensional spaces. The notes correctly distinguish induced matrix norms from the Frobenius norm and emphasize submultiplicativity as the property that enables perturbation and convergence estimates.

The discussion of decompositions covers eigendecomposition, Schur decomposition, and the SVD. The notes correctly stress that eigendecomposition is not universally available in a numerically useful form, whereas Schur decomposition and the SVD exist for every square or rectangular matrix over the appropriate field. This distinction is important for later chapters: spectral reasoning may be analytically convenient, but robust numerical algorithms must generally rely on orthogonal or unitary transformations.

The treatment of floating-point arithmetic introduces the standard relative-error model and machine precision, with double precision quoted as approximately $1.11\times 10^{-16}$. Rounding, truncation, cancellation, overflow, and non-associativity are presented as distinct error mechanisms. The discussion also points toward reduced-precision and mixed-precision computation, which is particularly relevant to machine learning and data assimilation. Mixed-precision numerical linear algebra has been surveyed in detail elsewhere (Bonavita et al., 2022).

The conceptual distinction between conditioning and stability is one of the strongest parts of the chapter. The condition number $\kappa(A)=\|A\|\|A^{-1}\|$ measures sensitivity of the problem, while backward stability measures whether an algorithm solves a nearby problem. The notes correctly make the key point that a backward-stable algorithm can still produce a large forward error when the underlying problem is ill-conditioned. This is illustrated with nearly singular matrices and with the Hilbert matrix.

Direct methods are then introduced through LU and Cholesky factorization. Gaussian elimination has $\mathcal{O}(n^3)$ factorization cost and $\mathcal{O}(n^2)$ triangular-solve cost in the dense case. Partial pivoting is presented as essential for practical backward stability, producing $PA=LU$. For SPD matrices, Cholesky reduces the cost to approximately $n^3/3$ flops and halves storage relative to general LU. These observations motivate the transition to iterative methods for large sparse systems.

The notes’ cost argument is straightforward and effective: a sparse matrix–vector product costs $\mathcal{O}(\operatorname{nnz}(A))$, whereas factorization can introduce substantial fill-in. This distinction becomes decisive for PDE systems with millions of unknowns, graph Laplacians, and sparse machine-learning operators.

Stationary iterations and spectral convergence

Jacobi, Gauss–Seidel, SOR, and Richardson are unified through the splitting $\ell^2$0 and the fixed-point iteration

$\ell^2$1

The fundamental convergence result is correctly stated: convergence for every initial vector is equivalent to $\ell^2$2, where $\ell^2$3. The error satisfies $\ell^2$4, so the spectral radius controls asymptotic convergence. A subordinate norm satisfying $\ell^2$5 provides a practical sufficient condition and an explicit error estimate.

The applications to strictly diagonally dominant and SPD matrices are useful. Jacobi and Gauss–Seidel converge under strict row diagonal dominance, while Gauss–Seidel converges for SPD matrices. For SOR, the necessary condition $\ell^2$6 is correctly identified, and the notes state the classical SPD convergence result over this interval. Richardson iteration is analyzed for SPD matrices, including the optimal step size and contraction factor

$\ell^2$7

The analysis of the one-dimensional Poisson problem gives the notes their first substantial quantitative result. For

$\ell^2$8

the eigenvalues are

$\ell^2$9

and $\ell^\infty$0. The Jacobi iteration eigenvalues are

$\ell^\infty$1

so

$\ell^\infty$2

Thus Jacobi requires $\ell^\infty$3 iterations to reduce error by a fixed asymptotic factor. The notes identify the mechanism precisely: high-frequency modes are damped rapidly, whereas smooth low-frequency modes persist.

For the stated one-dimensional Poisson setting, the Gauss–Seidel factor is the square of the Jacobi factor,

$\ell^\infty$4

which improves the constant but not the asymptotic order. Optimally relaxed SOR achieves

$\ell^\infty$5

and

$\ell^\infty$6

The numerical example with $\ell^\infty$7 is particularly informative:

Method	Spectral radius	Estimated iterations for a $\ell^\infty$8 reduction
Jacobi	$\ell^\infty$9	approximately $1.11\times 10^{-16}$0
Gauss–Seidel	$1.11\times 10^{-16}$1	approximately $1.11\times 10^{-16}$2
Optimal SOR	$1.11\times 10^{-16}$3	approximately $1.11\times 10^{-16}$4

The implication is not simply that SOR is “better.” Even 110 iterations are inferior to the linear-time Thomas algorithm for this tridiagonal problem. The real motivation for iterative methods appears when direct factorization loses sparsity or becomes infeasible, particularly in two and three dimensions.

Sparse matrices from PDEs, graphs, machine learning, and data assimilation

The second chapter derives sparse matrices from finite differences. Taylor expansion produces the standard first- and second-derivative formulas, and the one-dimensional Dirichlet Poisson problem yields a tridiagonal matrix with stencil $1.11\times 10^{-16}$5. The stencil determines both sparsity and bandwidth, linking discretization choices directly to solver complexity.

For the two-dimensional Poisson equation, the five-point stencil gives

$1.11\times 10^{-16}$6

The resulting matrix has at most five nonzeros per row, is symmetric, and is SPD under homogeneous Dirichlet conditions. However, under lexicographic ordering its bandwidth grows like the one-dimensional grid width. With $1.11\times 10^{-16}$7 unknowns, naive banded LU has cost $1.11\times 10^{-16}$8, whereas nested dissection can reduce the direct-solver complexity to approximately $1.11\times 10^{-16}$9.

The graph-Laplacian discussion makes the PDE analogy explicit. For an undirected graph,

$\kappa(A)=\|A\|\|A^{-1}\|$0

is symmetric positive semidefinite, with

$\kappa(A)=\|A\|\|A^{-1}\|$1

The multiplicity of its zero eigenvalue equals the number of connected components, while the Fiedler value measures algebraic connectivity. The path graph produces a near-tridiagonal operator, whereas the star graph produces the spectrum

$\kappa(A)=\|A\|\|A^{-1}\|$2

This comparison demonstrates that graph topology affects both sparsity patterns and spectral behavior. The notes then connect graph Laplacians to spectral clustering, graph neural networks, electrical networks, and diffusion-based covariance models.

Machine-learning examples focus on least squares and logistic regression. The normal equations

$\kappa(A)=\|A\|\|A^{-1}\|$3

show that sparse design matrices can produce denser Gram matrices. The notes correctly identify two hazards: forming $\kappa(A)=\|A\|\|A^{-1}\|$4 can destroy sparsity, and it squares the condition number,

$\kappa(A)=\|A\|\|A^{-1}\|$5

The logistic-regression Hessian $\kappa(A)=\|A\|\|A^{-1}\|$6 provides a corresponding second-order example. The broader data-assimilation discussion emphasizes that empirical covariance estimation can create spurious long-range correlations, motivating localization and sparse approximations. Covariance localization and diffusion-based covariance models are active topics in data assimilation (Skrunes et al., 2023).

The sparse-matrix toolbox covers CSR and CSC storage, bandwidth, profile, fill-in, and reordering. A particularly important practical claim is that the inverse of a sparse matrix can be dense even when its LU or Cholesky factors remain relatively sparse. This is why numerical algorithms apply factorizations or operators rather than explicitly forming inverses.

Figure 1: The sparsity pattern of a matrix, its inverse, and its LU factors illustrates why explicit inversion is usually more expensive than factorization.

Iterative methods as filters

The third chapter interprets stationary methods through their action on spectral modes. This is a major conceptual contribution of the notes because it provides a common language for PDE relaxation, graph diffusion, and optimization.

For the one-dimensional Laplacian, Jacobi acts on the sine eigenvectors with factors

$\kappa(A)=\|A\|\|A^{-1}\|$7

High-frequency modes have small $\kappa(A)=\|A\|\|A^{-1}\|$8 and are damped quickly; low-frequency modes have $\kappa(A)=\|A\|\|A^{-1}\|$9 and decay slowly. Weighted Jacobi with $\mathcal{O}(n^3)$0 minimizes the worst-case high-frequency damping factor and gives

$\mathcal{O}(n^3)$1

This is the smoothing property used by multigrid. The notes correctly distinguish smoothing from convergence: a method can be an effective smoother while remaining an inefficient standalone solver.

For graph Laplacians, Jacobi becomes a random-walk or diffusion operator. With $\mathcal{O}(n^3)$2, repeated application of $\mathcal{O}(n^3)$3 damps high graph frequencies and preserves smooth modes. The Fiedler vector is the slowest nontrivial mode and encodes a natural graph partition. The two-triangle example is numerically concrete: the Fiedler value is approximately $\mathcal{O}(n^3)$4, and the associated vector changes sign across the bridge connecting the two communities.

Richardson iteration is identified exactly with gradient descent on a quadratic objective. For an SPD Hessian or normal matrix, the error in eigenmode $\mathcal{O}(n^3)$5 is multiplied by

$\mathcal{O}(n^3)$6

The optimal fixed step is

$\mathcal{O}(n^3)$7

with worst-case factor $\mathcal{O}(n^3)$8. This establishes conditioning as the direct analogue of spectral stiffness in PDEs.

The notes extend this analysis to spectral bias in machine learning. They are careful to state the necessary qualification: the association between frequency and eigenvalue depends on the data representation, kernel, or feature distribution. Under common settings, low-frequency or global patterns correspond to larger covariance or kernel eigenvalues and are learned earlier, while fine-scale components associated with smaller eigenvalues are learned later. The shared mechanism is spectral filtering, although the interpretation of “frequency” differs between a spatial PDE basis and a learned data basis.

CG, Lanczos, and regularization

The fourth chapter replaces fixed splittings with adaptive Krylov subspaces. For SPD $\mathcal{O}(n^3)$9, CG minimizes the quadratic functional over

$\mathcal{O}(n^2)$0

where

$\mathcal{O}(n^2)$1

The notes give the standard CG recurrence, residual orthogonality, $\mathcal{O}(n^2)$2-conjugacy of search directions, and the best-approximation property in the $\mathcal{O}(n^2)$3-norm. In exact arithmetic, termination occurs in at most the degree of the minimal polynomial, and therefore in at most $\mathcal{O}(n^2)$4 steps.

The principal quantitative advantage is the Chebyshev estimate

$\mathcal{O}(n^2)$5

For Poisson discretizations with $\mathcal{O}(n^2)$6, this changes the iteration scale from $\mathcal{O}(n^2)$7 for steepest descent to $\mathcal{O}(n^2)$8 for CG. The bound is conservative when eigenvalues cluster, but it captures the decisive square-root improvement.

Lanczos is presented as the eigenvalue counterpart of CG. It constructs an orthonormal basis of a Krylov subspace and projects a large symmetric matrix onto a small tridiagonal matrix $\mathcal{O}(n^2)$9. Ritz values approximate extremal eigenvalues first, allowing estimates of spectral endpoints and condition numbers using only matrix–vector products. The notes further explain the exact-arithmetic equivalence between CG and Lanczos: the normalized CG residuals form the Lanczos basis up to alternating signs, and the Lanczos recurrence coefficients can be recovered from CG scalars.

For graph Laplacians, deflation of the known nullspace vector $PA=LU$0 makes the Fiedler eigenvalue extremal on the deflated subspace. In the two-triangle example, three Lanczos steps recover the nonzero spectral information to the displayed precision, including the Fiedler value $PA=LU$1 and the largest eigenvalue $PA=LU$2. This is a strong illustrative result, although it relies on a very small graph with only a few distinct eigenvalues.

The machine-learning subsection interprets least-squares solvers as spectral filters. The SVD solution amplifies noise through factors $PA=LU$3. Ridge regression replaces these with the bounded filter

$PA=LU$4

Early-stopped CG produces a polynomial filter $PA=LU$5 that is data-adaptive and resolves well-determined directions first. The notes correctly emphasize semi-convergence: the approximation error decreases initially, but noise amplification eventually increases. Early stopping is therefore regularization, not merely an incomplete solve, and a stopping rule remains necessary.

GMRES, nonnormality, PageRank, and LSQR

The fifth chapter addresses nonsymmetric systems, where CG and Lanczos no longer apply. Arnoldi constructs

$PA=LU$6

with $PA=LU$7 upper Hessenberg. GMRES selects the iterate minimizing the residual over the current Krylov affine space. By orthonormality, the large residual minimization reduces exactly to

$PA=LU$8

The notes correctly identify full-GMRES residual monotonicity, the polynomial representation $PA=LU$9, and finite termination in the degree of the relative minimal polynomial. Incremental Givens rotations are used to update the reduced least-squares problem efficiently.

The discussion of nonnormality is especially important. For diagonalizable $n^3/3$0, an eigenvalue-based bound contains the factor $n^3/3$1, which can be enormous. Thus clustered eigenvalues do not guarantee rapid GMRES convergence. The companion-matrix example with all eigenvalues equal to $n^3/3$2 demonstrates stagnation for two iterations followed by exact termination at the third. The field of values provides a more robust alternative, with the Crouzeix–Palencia inequality yielding the disk bound

$n^3/3$3

when $n^3/3$4 lies in a disk centered at $n^3/3$5 of radius $n^3/3$6 that excludes the origin. The implication is precise: for nonnormal systems, pseudospectral or numerical-range information can be more predictive than eigenvalues alone.

PageRank is formulated as the dominant eigenvector of the Google matrix

$n^3/3$7

Although $n^3/3$8 is dense, it can be applied without formation:

$n^3/3$9

for probability vectors $\mathcal{O}(\operatorname{nnz}(A))$0. The cost is one sparse product with $\mathcal{O}(\operatorname{nnz}(A))$1 plus vector operations. All non-dominant eigenvalues satisfy $\mathcal{O}(\operatorname{nnz}(A))$2, so the power iteration converges with asymptotic factor at most $\mathcal{O}(\operatorname{nnz}(A))$3, independently of graph size. This gives a clear modeling–computational tradeoff: increasing $\mathcal{O}(\operatorname{nnz}(A))$4 makes the ranking more faithful to the link graph but slows convergence approximately like $\mathcal{O}(\operatorname{nnz}(A))$5.

The final subsection introduces Golub–Kahan bidiagonalization and LSQR for rectangular least squares. LSQR reduces the problem to

$\mathcal{O}(\operatorname{nnz}(A))$6

with a small bidiagonal matrix $\mathcal{O}(\operatorname{nnz}(A))$7. This avoids explicitly forming $\mathcal{O}(\operatorname{nnz}(A))$8, preserving the conditioning of $\mathcal{O}(\operatorname{nnz}(A))$9 and its sparse operator structure. The notes’ warning example is numerically important: if

$\ell^2$00

then forming $\ell^2$01 in double precision can destroy rank information when $\ell^2$02, already occurring around $\ell^2$03. The practical implication is unambiguous: CG on the normal equations may be useful analytically, but LSQR or related bidiagonalization methods are preferable computationally.

Preconditioning, Schwarz methods, and multigrid

The final chapter treats preconditioning as spectral transformation. Left, right, and symmetric preconditioning are distinguished, and for SPD problems the symmetrically preconditioned operator

$\ell^2$04

is used to retain SPD structure. The preconditioned CG recurrence requires only solves with $\ell^2$05, not explicit formation of $\ell^2$06.

The notes give a useful negative result for Jacobi preconditioning. For the one-dimensional Laplacian, the diagonal is a scalar multiple of the identity, so Jacobi preconditioning merely rescales all eigenvalues and leaves the condition number unchanged. This illustrates that diagonal scaling cannot address global low-frequency error.

The general subspace-correction framework

$\ell^2$07

unifies pointwise Jacobi, block methods, Schwarz preconditioners, and coarse-space methods. Additive Schwarz uses overlapping subdomain solves and combines the corrections using a partition of unity. In the one-dimensional continuous model, alternating Schwarz has contraction factor

$\ell^2$08

for overlapping subdomains $\ell^2$09 and $\ell^2$10. The formula exposes the role of overlap: as the overlap shrinks, $\ell^2$11; increasing overlap accelerates information transfer.

The multigrid analysis returns to the smoothing property. Weighted Jacobi with $\ell^2$12 reduces high-frequency modes by at least a factor of $\ell^2$13 per sweep, but leaves smooth modes nearly unchanged. Coarse-grid correction represents and removes those smooth modes. With prolongation $\ell^2$14, restriction $\ell^2$15, and Galerkin coarse operator $\ell^2$16, the correction is

$\ell^2$17

The notes correctly emphasize that coarse-grid correction alone cannot converge. Its error-propagation operator is a projection, with eigenvalues $\ell^2$18 and $\ell^2$19; therefore unresolved components remain untouched. Smoothing and coarse correction are complementary rather than interchangeable.

Recursively applying this process yields a V-cycle. Under the standard assumptions for elliptic problems, the work per cycle is $\ell^2$20 in one dimension, and the convergence factor is approximately independent of mesh size. The resulting multigrid preconditioner can reduce the $\ell^2$21 CG iteration growth to an essentially mesh-independent iteration count. The notes state this conclusion correctly at the conceptual level, although a rigorous mesh-independent theorem requires assumptions on the operator, transfer operators, smoothing scheme, coefficient regularity, and boundary treatment.

Limitations and open questions

The manuscript’s principal limitation is that it is lecture-note material rather than a validated research study. Most numerical claims are theoretical examples or hand-selected small systems; there is no systematic benchmarking against contemporary sparse direct solvers, AMG packages, restarted GMRES variants, flexible preconditioning, or modern LSQR implementations.

Several statements require correction or qualification.

• The displayed forward-difference formula labels the truncation remainder as $\ell^2$22 after retaining an $\ell^2$23 term. The approximation is first-order accurate, with residual $\ell^2$24 only in the unscaled Taylor remainder.
• The claim that Gauss–Seidel eigenvalues are generally the squares of Jacobi eigenvalues is not valid for arbitrary splittings. It requires additional structural assumptions, such as those satisfied by the standard consistently ordered Poisson model.
• The statement that reordering is a similarity transformation is inaccurate for the usual symmetric row-and-column permutation $\ell^2$25. Such a permutation preserves eigenvalues, but a row-only permutation generally does not.
• Sparse $\ell^2$26 does not imply that $\ell^2$27 remains significantly sparse. Its sparsity depends on feature co-occurrence; in many applications the Gram matrix is effectively dense.
• The PageRank discussion must distinguish row- and column-stochastic conventions carefully. The eigenvector orientation, stochasticity condition, and teleportation term must be consistent throughout.
• The notes contain numerous malformed LaTeX expressions, missing braces, inconsistent matrix dimensions, and notation changes. These are not merely cosmetic defects in a numerical-linear-algebra text because they obscure whether a formula is mathematically intended or computationally implementable.
• The multigrid claims are stated at a level appropriate for an introductory course but omit the assumptions needed for mesh-independent convergence, especially for variable-coefficient, anisotropic, nonsymmetric, or unstructured problems.
• The machine-learning discussion of spectral bias is suggestive but not universal. The identification of low-frequency functions with large kernel or covariance eigenvalues depends on the representation and data distribution; it should not be presented as an unconditional equivalence.

The main open technical question is how the proposed unified exposition should be extended from model operators and textbook examples to nonsmooth coefficients, anisotropic PDEs, indefinite systems, nonsymmetric preconditioners, stochastic optimization, and covariance operators arising from realistic data-assimilation workflows.

Conclusion

The notes provide a coherent numerical-linear-algebra curriculum centered on the interaction between sparsity and spectrum. Their strongest contribution is conceptual unification: finite-difference Laplacians, graph diffusion, gradient descent, Krylov solvers, PageRank, LSQR, Schwarz methods, and multigrid are presented as variations on structured operator application, spectral filtering, and subspace correction.

The quantitative narrative is consistent across the chapters. Direct factorization is robust but vulnerable to cost and fill-in; stationary methods are cheap but spectrally limited; CG and GMRES exploit expanding Krylov spaces; LSQR avoids the numerical liabilities of normal equations; and preconditioning changes the spectrum so that Krylov convergence becomes scalable. Despite technical and editorial issues requiring correction, the manuscript offers a useful foundational framework for advanced study of numerical linear algebra in PDEs, machine learning, networks, and data assimilation (2608.21234).