layout: true class: typo, typo-selection --- count: false class: nord-dark, middle, center # 🚀🌈 Preconditioned Krylov Subspace Methods for Large-Scale Nonsymmetric Linear Systems @luk036 👨‍💻 · 2026 📅 --- ### 📋 Agenda .pull-left[ **Part 1: Introduction** 🎯 - Background - Motivation **Part 2: Krylov Methods** 🔄 - Subspace Methods - Preconditioning ] .pull-right[ **Part 3: Applications** 🔬 - Numerical Results - Conclusions ] --- ### 📈🔥 The Core Problem: Large-Scale Linear Systems * We aim to solve a system of linear equations: 🎯 * ${\color{blue} A}x = {\color{coral} b}$ * Where ${\color{blue} A} \in \mathbb{R}^{n \times n}$ is a large, sparse, and nonsymmetric matrix, $x \in \mathbb{R}^{n}$ is the unknown solution vector, and ${\color{coral} b} \in \mathbb{R}^{n}$ is the known right-hand side vector. * **Characteristics of these systems:** * **Large:** Sizes range from 10⁵ to 10⁷ unknowns, typical of modern industrial applications. 🏗️ * **Sparse:** The majority of entries in ${\color{blue} A}$ are zero. 🕸️ * **Often Nonsymmetric:** This arises from various reasons such as the multiphysics nature of PDEs, sophisticated boundary or jump conditions, or the discretization methods themselves (e.g., convection-diffusion equations, Petrov-Galerkin methods). 🔄 --- * **Why Iterative Methods?** * Direct methods (e.g., Cholesky decomposition, Gaussian elimination) are often too computationally expensive and memory-intensive for large, sparse systems. 💸 * Iterative methods, like Krylov subspace methods, find approximate solutions incrementally, reducing computational and storage costs. ⏳ --- ### 🧱✨ KSP: Building Blocks for Solutions * **Definition:** The $k$-th Krylov subspace ${\color{darkmagenta} K_k(A,v)}$ generated by matrix ${\color{blue} A}$ and vector ${\color{green} v}$ is: * ${\color{darkmagenta} K_k(A,v)} = \text{span}\{{\color{green} v}, {\color{blue} A}{\color{green} v}, {\color{blue} A}^2{\color{green} v}, \ldots, {\color{blue} A}^{k-1}{\color{green} v}\}$ * To solve ${\color{blue} A}x={\color{coral} b}$, ${\color{green} v}$ is typically the initial residual ${\color{red} r_0} = {\color{coral} b} - {\color{blue} A}x_0$. * **How they work:** * KSP methods incrementally find approximate solutions within ${\color{darkmagenta} K_k(A,v)}$. 🔄 * They avoid explicit matrix formation by primarily relying on matrix-vector products, enabling "matrix-free" methods. 🎯 * **Popular KSP Methods for Nonsymmetric Systems:** * GMRES 🏆 * TFQMR 🚀 * BiCGSTAB ⚡ * QMRCGSTAB 🌟 --- * **Core Iteration Procedures:** * **Arnoldi Iteration:** Constructs an *orthogonal* basis for $K_k(A,v)$. * The KSP method GMRES is based on Arnoldi iteration. * Has a $k$-term recurrence, meaning computational cost increases with $k$. 📉 * Often requires restarts (e.g., GMRES(30)) to manage cost, which can unfortunately undermine convergence. ⚠️ * **Bi-Lanczos Iteration:** Offers an alternative for constructing basis. * Enjoys a *three-term recurrence*, meaning fixed computational cost per iteration. 💰 * Requires two matrix-vector multiplications per iteration, but the basis is *not orthogonal*. 🔄 * Used by methods like BiCG and QMR, and their transpose-free variants (TFQMR, BiCGSTAB, QMRCGSTAB). * Can suffer from breakdown, though this is rare in practice. ⚠️ --- ### 🚀💡 The Power of Preconditioning * **Why are Preconditioners Needed?** 🤔 * The convergence of KSP methods can be improved significantly by using preconditioners. * The effectiveness of Krylov methods heavily depends on the *conditioning* of the matrix ${\color{blue} A}$. Ill-conditioned matrices lead to slow convergence. 🐌 * **What is a Preconditioner (${\color{orange} M}$)?** * A matrix or transformation ${\color{orange} M}$ whose inverse ${\color{orange} M}^{-1}$ *approximates* the inverse of ${\color{blue} A}$ (${\color{blue} A}^{-1}$). * The application of ${\color{orange} M}^{-1}$ to a vector must be computationally efficient. * **Goal of Preconditioning:** 🎯 * To transform the original system ${\color{blue} A}x={\color{coral} b}$ into an equivalent system (e.g., ${\color{orange} M}^{-1}{\color{blue} A}x = {\color{orange} M}^{-1}{\color{coral} b}$ or ${\color{blue} A}{\color{orange} M}^{-1}y={\color{coral} b}$) such that the preconditioned matrix (${\color{orange} M}^{-1}{\color{blue} A}$ or ${\color{blue} A}{\color{orange} M}^{-1}$) has a **smaller condition number** than ${\color{blue} A}$. * A smaller condition number directly translates to a faster rate of convergence for iterative solvers. 🚀 --- * **Types of Preconditioning (for nonsymmetric systems):** 🔀 * **Left Preconditioning:** Solve ${\color{orange} M}^{-1}{\color{blue} A}x = {\color{orange} M}^{-1}{\color{coral} b}$. * The convergence criterion is typically based on the preconditioned residual $‖{\color{orange} M}^{-1}{\color{red} r_k}‖$. * This can differ significantly from the *true residual* $‖{\color{red} r_k}‖$ if $‖{\color{orange} M}^{-1}‖$ is far from 1, leading to **false stagnation** or **premature termination**. ⚠️ * **Right Preconditioning (Recommended) 👍:** Solve ${\color{blue} A}{\color{orange} M}^{-1}y = {\color{coral} b}$, then $x = {\color{orange} M}^{-1}y$. * The crucial advantage is that it **does not alter the true residual** ${\color{red} r_k}$. * The stopping criteria can use $‖{\color{red} r_k}‖$ with little or no extra computational cost. * This makes right preconditioning generally **more reliable** for large-scale systems. ✅ * **Two-Sided/Symmetric Preconditioning:** Involves applying preconditioners to both sides (e.g., $QAP^{-1}(Px) = Qb$). While it can preserve symmetry for symmetric matrices, it also alters the norm of the residual for nonsymmetric matrices. 🔄 --- ### 🧰🌈 Our Toolkit: A Look at Preconditioner Families We systematically compared KSP methods with three classes of state-of-the-art preconditioners: .mermaid[ <pre> graph LR P["Preconditioner Families 🧰"] --> GS["Gauss-Seidel 🔄\nsimple, no setup"] P --> ILU["ILU 🏗️\nILU0 / ILUT / ILUTP / MILU"] P --> AMG["AMG 🌐\nBoomerAMG / ML"] style P fill:#f3e5f5,stroke:#7b1fa2,stroke-width:3px,color:#2e3440 style GS fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px,color:#2e3440 style ILU fill:#fff3e0,stroke:#e65100,stroke-width:3px,color:#2e3440 style AMG fill:#e3f2fd,stroke:#1565c0,stroke-width:3px,color:#2e3440 </pre> ] * **Gauss-Seidel (GS)** 🔄: * One of the simplest stationary iterative methods. * Easy to implement and requires virtually no setup time (in serial). * Parameter-free (when the relaxation parameter $\omega = 1$). * **Limitation:** Requires nonzero diagonal entries in $A$. * Generally outperformed by ILU0 for the tested large linear systems. --- * **Incomplete LU (ILU) Factorization** 🏗️: * Performs an approximate factorization $A \approx \tilde{L}\tilde{U}$, where $\tilde{L}$ and $\tilde{U}$ are much sparser than full LU factors. * **Variants studied:** * **ILU0 / ILU(0):** Simplest form, introduces no fill-in, preserving $A$'s sparsity pattern. * **ILUT:** (ILU with dual thresholding) Extends ILU(k) by preserving fills based on levels and thresholding numerical values. * **ILUTP:** (ILU with partial pivoting) Mitigates breakdown issues when $A$ has zeros on the diagonal (common in saddle-point problems). * **Drawback:** Number of nonzeros in $\tilde{L}$ and $\tilde{U}$ can grow super-linearly with matrix size, leading to super-linear growth in setup and solve times for very large systems. * **Multilevel ILU (MILU):** Better controls the growth of nonzeros, leading to *linear* growth of setup and solve times for large sparse systems. * **Key Advantage:** Succeeds for **saddle-point problems** where AMG may fail. * Empirically, MILU (ILUPACK) often outperformed SuperLU's ILUTP and ILU0 for large systems. --- * **Algebraic Multigrid (AMG)** 🌐: * The most sophisticated preconditioners, accelerating convergence by constructing a series of coarser representations. * More expensive than GS and ILU in terms of setup time and per-iteration cost, but offer significantly faster convergence. * **Types studied:** * **Classical AMG:** Implemented by BoomerAMG (part of HYPRE). Uses "classical coarsening" strategies. * **Smoothed Aggregation:** Implemented by ML. Coarsening is done by accumulating aggregates. * **BoomerAMG Tuning (Crucial for 3D problems):** * BoomerAMG offers various **coarsening strategies** (Falgout, PMIS, HMIS) and **interpolation techniques** (classical, extended+classical, FF1). * Falgout coarsening is good for 2D, but less efficient for 3D. PMIS and HMIS significantly improve memory and runtime for 3D problems. * FF1 interpolation performed better than extended+classical interpolation with HMIS/PMIS coarsening. * **Limitation:** Both AMG methods may **fail** for some cases, especially for saddle-point-like problems (matrices with many zero diagonal entries). --- ### ⚙️🔧 KSP Methods: Trade-offs & Mechanics | Method | Basis Construction | Mat-Vec Products (per iter.) | Recurrence Type | Storage (vectors) | :--------- | :----------------------- | :------------------------------- | :-------------- | :---------------- | **GMRES** | Arnoldi (orthogonal) | 1 ($AM^{-1}$) | $k$-term | $k + 5$ | **BiCGSTAB** | Transpose-free Bi-Lanczos (non-orthogonal) | 2 ($AM^{-1}$) | 3-term | $10$ | **TFQMR** | Same as BiCGSTAB | 2 ($AM^{-1}$) | 3-term | $8$ | **QMRCGSTAB** | Same as BiCGSTAB | 2 ($AM^{-1}$) | 3-term | $13$ --- * **Computational Cost & Convergence Rate:** * While GMRES requires one matrix-vector multiplication per iteration and the other methods require two, the reduction of error in one iteration of a bi-Lanczos-based KSP method is approximately equal to that of two iterations in an Arnoldi-based KSP method. Therefore, the costs of these methods are **comparable** in terms of the number of matrix-vector multiplications required to achieve similar error reduction. * GMRES, minimizing the 2-norm of the residual, *can* converge faster in terms of matrix-vector multiplications, especially with an effective preconditioner. 🚀 * **Robustness Considerations:** * Theoretically, Arnoldi iteration is more robust due to its use of an orthogonal basis, but restarts can undermine convergence. * Bi-Lanczos iterations may break down if $v_{k+1}^T w_{k+1} = 0$, but this is rare in practice, especially with a good preconditioner. * **Crucially, effective preconditioners can overcome many of the inherent disadvantages of each KSP method**, leading to smooth and rapid convergence. --- ### ✅ Choosing the Best Method: Recommendations from Empirical Study Our systematic comparison across 16 large-scale nonsymmetric linear systems (up to 10⁷ unknowns) derived from various PDE discretizations (FEM, AES-FEM, GFD, FDM, DG, Mixed FEM, Stokes) provides the following practical guidelines: 📊 .mermaid[ <pre> graph TD A["Large sparse system Ax = b\n(nonsymmetric)"] --> B{"Well-conditioned?\nAMG available?"} B -->|"Yes ✅"| C["GMRES + Classical AMG\n(right preconditioning)"] B -->|"No / saddle-point ⚠️"| D["BiCGSTAB + MILU\n(right preconditioning)"] C --> E["BoomerAMG: HMIS + FF1"] D --> F["Or ILUTP for\nmoderate sizes"] style A fill:#f3e5f5,stroke:#7b1fa2,stroke-width:3px,color:#2e3440 style B fill:#ffcdd2,stroke:#c62828,stroke-width:3px,color:#2e3440 style C fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px,color:#2e3440 style D fill:#fff3e0,stroke:#e65100,stroke-width:3px,color:#2e3440 style E fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px,color:#2e3440 style F fill:#fff3e0,stroke:#e65100,stroke-width:3px,color:#2e3440 </pre> ] --- * **Primary Recommendation: For Very Large, Reasonably Well-Conditioned Linear Systems** * **Use GMRES with Classical AMG as a Right Preconditioner**. * **Why?** GMRES (without restart) minimizes the $L2$ norm of the residual, and with an effective AMG, it converges quickly, keeping the cost of orthogonalization low. * **AMG Tuning is CRITICAL (especially for 3D problems):** * For BoomerAMG (classical AMG, e.g., in HYPRE), the most robust option found was **HMIS coarsening with FF1 interpolation**. This combination outperforms PMIS+FF1 for ill-conditioned systems and significantly outperforms smoothed-aggregation AMG (ML). * AMG-preconditioned KSP methods exhibit **nearly linear asymptotic growth** ($slope ≈ 1$) in runtime with respect to problem size, making them highly scalable for very large problems. --- * **Secondary Recommendation: If AMG is Unavailable or Fails (e.g., Saddle-Point Problems)** * **Use BiCGSTAB with Multilevel ILU (MILU) as a Right Preconditioner**. * **Why MILU?** It effectively controls fill-in, leading to nearly linear growth of setup and solve times, making it suitable for large sparse systems. Crucially, MILU succeeds for **saddle-point-like problems** (matrices with many zero diagonal entries), where AMG methods may fail. * **Why BiCGSTAB?** Without an effective AMG, BiCGSTAB typically outperforms GMRES and is often more efficient and robust than restarted GMRES. * **Alternative:** If MILU is unavailable (e.g., in PETSc), ILUTP (e.g., in SuperLU or MATLAB's built-in $ilu$ function) can be used for moderate-sized systems, as it also handles saddle-point problems. However, ILUTP's setup times may grow super-linearly for very large problems. --- * **General Principle: Prioritize Right Preconditioning** * **Right preconditioning is generally more reliable** than left preconditioning for large-scale systems. * It ensures that the stopping criteria based on the residual norm $‖r_k‖$ accurately reflect the true error, preventing misleading convergence behavior like false stagnation or premature termination. * **Recommendation:** When using software like PETSc, *explicitly set the option to use right preconditioning*, as left preconditioning is often the default. --- ### 🤖 AI for Solvers: Graph Neural Preconditioners * **Concept:** A novel approach proposing to use **Graph Neural Networks (GNNs)** to learn and approximate the matrix inverse ${\color{blue} A}^{-1}$ (which serves as the preconditioner ${\color{orange} M}$). 🧠 * **Motivation:** Inspired by the universal approximation property of neural networks and their widespread success in artificial intelligence. * **Why GNNs?** A sparse matrix $A$ can be naturally interpreted as a graph adjacency matrix. GNNs leverage this graph structure, performing convolutions in node neighborhoods, which is well-suited for sparse matrices. * **Integration with Krylov Solvers:** Since ${\color{orange} M}$ (the GNN) is a *nonlinear* operator 🤔, standard Krylov subspace theory (which relies on linearity) breaks down. Thus, GNP is used with **Flexible GMRES (FGMRES)**. 💡 * FGMRES is designed to accommodate preconditioners that change (or are nonlinear) in each Arnoldi step. --- * **Key Technical Contributions & Innovations:** 1. **Convergence Analysis for FGMRES with Nonlinear $M$:** Provides a novel *a posteriori* analysis. The residual $r_m$ is bounded by: * $‖r_m‖_2 \le \kappa_2(\mathbf{X})\epsilon^{(m)}(\mathbf{\Lambda})‖r_0‖_2 + ‖\mathbf{Q}_m\mathbf{Q}_m^{\top} - \widetilde{\mathbf{Q} }_m\widetilde{\mathbf{Q} }_m^{\top}‖_2‖r_0‖_2$ * The second term highlights the impact of the nonlinear preconditioner, showing how closely the subspace generated by the GNN ($\mathbf{Q}_m$) matches that of an ideal linear preconditioner ($\widetilde{\mathbf{Q} }_m$). 2. **Effective Training Data Generation:** * Addresses the challenge of sampling $x$ such that $b=Ax$ covers important spectral information (especially the bottom eigen-subspace). * Samples $x$ from a mixture of $\mathcal{N}(\mathbf{0}, \mathbf{\Sigma}_m^{\mathbf{x} })$ (derived from Arnoldi process for extreme eigenvectors) and $\mathcal{N}(\mathbf{0}, \mathbf{I}_n)$ (uniform sphere coverage). This mix balances ease of training with preconditioner. 3. **Scale-Equivariant GNN Architecture:** * Introduces a parameter-free scaling $s(\cdot) = \frac{\sqrt{n} }{\tau}\cdot$ and corresponding back-scaling $s^{-1}(\cdot) = \frac{\tau}{\sqrt{n} }\cdot$, where $\tau = ‖b‖_2$. * This design imposes an inductive bias that ensures **scale-equivariance** ($\mathbf{M}(\alpha\mathbf{b})=\alpha\mathbf{M}(\mathbf{b})$ for $\alpha \neq 0$), facilitating training when input $b$ vectors have vastly different scales. --- * **Performance Highlights (Evaluated on >800 SuiteSparse matrices):** * **Superior for Ill-Conditioned Problems:** Particularly effective for matrices with condition numbers $\geq 10^{16}$. * **High Robustness:** GNP exhibited very few failures (0.12% solution failures, 0% construction failures), significantly outperforming ILU (40% construction failures) and AMG (7% construction failures). * **Predictable Construction Time:** The GNN training time is nearly proportional to the matrix size, offering a much more predictable cost than ILU and AMG setup times, which can be highly erratic. * **Faster Execution:** GNP can be significantly faster than using GMRES as an inner preconditioner. --- ### 🔭🚀 Horizon: What's Next for Solvers? * **Current Study Limitations:** * The study primarily focused on serial performance and did not include comparisons of parallel performance or asymptotic complexity with respect to the number of cores. * Did not include comparisons with geometric multigrid (GMG) methods, which can sometimes outperform AMG. * Did not address the solution of nonsymmetric, rank-deficient linear systems, a challenging problem in its own right. --- * **Future Avenues for Graph Neural Preconditioners (GNP):** 1. **SPD Matrices:** Extend GNP to symmetric positive-definite (SPD) matrices. While flexible CG exists, flexible preconditioners can be fragile. Exploring split preconditioners or autoencoder networks might offer more robust solutions. 2. **Sequential Problems & Continual Learning:** Adapt GNP for sequences of evolving matrices, common in time-stepping or sequential problems. Continual learning can amortize the initial preconditioner construction cost by continuously adapting the GNN for slowly evolving data distributions. 3. **Scaling to Larger Matrices:** Current GPU memory limits restrict matrix size and network complexity. Future work can explore multi-GPU and distributed training of GNNs, potentially using neighborhood sampling to mitigate the "neighborhood explosion" problem in GNNs for larger graphs. 4. **Hyperparameter Optimization & Architecture Improvement:** While the reported results used a fixed set of hyperparameters, significant performance improvements are expected with problem-specific tuning and advances in GNN architectures. --- ### 🎯 Key Takeaways for Iterative Solver Success * **Preconditioners are Essential:** They are indispensable for achieving efficient and robust convergence of Krylov subspace methods when solving large, sparse linear systems. ✅ * **Right Preconditioning is Key:** Always prefer right preconditioning to ensure that convergence criteria accurately reflect the true residual, thereby avoiding false stagnation or premature termination. 👍 * **Primary Recommendation:** For large, reasonably well-conditioned linear systems, the most efficient method is often **GMRES coupled with a Classical Algebraic Multigrid (AMG) preconditioner**, particularly BoomerAMG with HMIS coarsening and FF1 interpolation. AMG's near-linear scalability offers significant advantages for increasing problem sizes. 🥇 * **Secondary Recommendation:** If AMG is unavailable or if the problem involves **saddle-point-like matrices** (with many zero diagonal entries), **BiCGSTAB with Multilevel ILU (MILU) as a right preconditioner** is a robust and effective alternative. 🥈 * **Emerging AI-Driven Solutions:** **Graph Neural Preconditioners (GNP)** represent a promising new direction 🤔. 🔮 --- count: false class: nord-dark, middle, center ## Q&A 🎤