Recycling Krylov Subspaces

• Recycling Krylov subspaces is a technique that reuses previously computed Krylov bases to accelerate the solution of sequences of related linear systems.
• It employs augmentation, deflation, and block strategies to incorporate spectral and solution information from past iterations, reducing computational overhead.
• Applications in time-dependent PDEs, bilevel learning, model reduction, and inverse imaging demonstrate significant gains in convergence and efficiency.

Recycling Krylov Subspaces

Recycling Krylov subspaces is a suite of methodologies enabling the reuse or adaptive augmentation of previously constructed Krylov subspaces to accelerate the iterative solution of related linear systems, matrix function evaluations, or nested optimization problems. This class of techniques is motivated by sequences of problems with slowly varying system matrices or right-hand sides, as encountered in time-dependent PDEs, shape optimization, parameter sweeps, inverse problems, model reduction, and bilevel learning. Recycling exploits approximate spectral or solution information from previous solves, strategically augmenting the search space for subsequent problems to reduce iteration counts and computational cost.

1. Foundational Principles of Krylov Subspace Recycling

Recycling Krylov subspaces is rooted in the classical theory of Krylov subspace methods, such as GMRES, MINRES, BiCG, and CG, which iteratively construct approximations in spaces of the form $$\mathcal{K}_m(A, r_0) = \mathrm{span}\{r_0, Ar_0, \ldots, A^{m-1}r_0\}$$. In practical applications, one often encounters sequences of linear systems,

$A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$

where $A^{(i)}$ and $b^{(i)}$ change slowly. Rather than restarting iterative solvers from scratch for each such system, recycling subspace methods extract useful spectral directions such as approximate invariant subspaces (e.g., Ritz or harmonic Ritz vectors) and transfer them, typically via augmentation or deflation, into the correction space for subsequent systems (Soodhalter et al., 2020, Parks et al., 2016).

Canonical recycling methods—such as GCRO-DR, RMINRES, and their block and flexible variants—explicitly maintain an augmented subspace $U \oplus \mathcal{K}_m((I-Q)A, (I-Q)r_0)$, where $U$ contains recycled basis vectors and $Q$ projects onto $A U$ (Parks et al., 2016, Soodhalter et al., 2020, Burke et al., 2023). The correction $x_m$ is computed by minimizing the residual over this space, requiring the solution of small projected least-squares subproblems.

A major conceptual advance is the recognition that recycling can be performed in both block and flexible Krylov settings, and that the selection of the recycled subspace is crucial for robust acceleration and must be regularly updated via spectral, subspace, or optimization-based criteria (Carlberg et al., 2015, Gaul et al., 2012, Pasha et al., 2023).

2. Methodologies for Recycling: Algorithms and Subspace Selection

The operational mechanism of recycling varies according to the specific method, system properties (e.g., symmetry, definiteness), and the application context:

• Augmentation: An extra subspace $U$ is added to the correction space. In each solve, $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$0 where $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$1 spans the new Krylov subspace, and $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$2 minimizes the chosen residual or error norm. This is the core of GCRO-DR and similar methods (Soodhalter et al., 2020, Parks et al., 2016).
• Deflation: A projector $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$3, with $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$4 orthonormal, is applied to "remove" slow-converging modes; the problem $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$5 is solved, and the solution is corrected to account for the contribution in the deflated space (Blondeau et al., 2023).
• Block and Flexible Schemes: Block Krylov methods build bases for multiple right-hand sides or shifted systems, and block recycling is achieved by augmenting or deflating in the block setup, often pairing Arnoldi or block-Lanczos recurrences with block-constrained subspace update (Parks et al., 2016, Burke, 2022, Soodhalter, 2014).
• Short-Recurrence and Storage-Efficient Approaches: When storage is prohibitive, methods such as short-recurrence IDR-based recycling (Neuenhofen, 2015) or memory-efficient block-partitioned schemes (Neuenhofen et al., 2016) are employed, leveraging polynomial recurrences or short-representation techniques.
• Biorthogonalization and Recycling for Nonsymmetric Systems: For nonsymmetric and dual system sequences, biorthogonal recycling subspaces (left and right) are maintained, as in RBiCG and RBiCGSTAB (Ahuja et al., 2010, Ahuja et al., 2014).
• Subspace Selection and Truncation: Strategies include extracting harmonic Ritz vectors, goal-oriented POD vectors (targeting error metrics or quantities of interest), SVD- or Schur-based truncations, or spectral components linked to specific function action approximations (Carlberg et al., 2015, Burke et al., 2023, Burke et al., 2023). The selection mechanism is often tailored to the application (e.g., minimizing hypergradient error in bilevel learning (Ehrhardt et al., 2024)).

A recurring theme is the emphasis on small, well-chosen augmentation or deflation spaces (few tens of vectors) to keep overhead manageable and maximize convergence benefit (Neuenhofen, 2015, Gaul et al., 2012).

3. Advanced Strategies: Hypergradient-Error and Model-Driven Recycling

Specialized recycling strategies have emerged to address nested or inverse problems, shifted system sequences, and model reduction:

• Bilevel Optimization and Hypergradient Error: In bilevel learning, the dominant cost is repeatedly solving similar linear systems for hypergradient computation. An innovative approach is recycling via Ritz generalized singular vectors obtained from a partial GSVD of $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$6, where $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$7 encodes upper-level derivative structure. This targets the error propagated into the hypergradient explicitly rather than the residual norm, facilitating a more direct and computationally effective stopping criterion (Ehrhardt et al., 2024).
• Model Reduction and Dual System Recycling: In interpolatory model reduction and some PDE contexts, the need to solve many slowly varying dual linear systems with both left and right actions motivates biorthogonal recycling, updating both left and right subspaces using block Lanczos recurrences and extracting harmonic Ritz pairs (Ahuja et al., 2010, Ahuja et al., 2014).
• Streaming and Dynamic Data: In large-scale dynamic or streaming inverse problems, the recycled subspace is adaptively compressed—via SVD truncation, greedy, or sparsity-based selection—each time the subspace grows beyond preset storage bounds. This maintains computational feasibility while retaining approximation quality over long sequences (Pasha et al., 2023).
• Block/Shifted System and Sylvester Equivalence: For sequences of shifted systems, block-Krylov recycling leverages the Sylvester-operator formulation to enable recycling even with unrelated right-hand sides or widely varying shifts. Shift-dependent harmonic Ritz vectors are used to update the recycle space in each cycle (Burke, 2022, Soodhalter, 2014, Soodhalter et al., 2013).

4. Implementation, Complexity, and Performance Considerations

Key implementation issues and trade-offs are dictated by the recycling approach, system structure, and hardware context:

• Orthogonalization and Storage Costs: Full Gram-Schmidt or Householder orthogonalization with possibly large recycled spaces is expensive; block or randomized sketching methods (e.g., in block-GCRO-DR, GMRES-SDR) dramatically reduce inner product costs and enhance cache/memory utilization (Parks et al., 2016, Burke et al., 2023, Burke et al., 2023).
• Augmentation Size vs. Overhead: Increasing the recycled space size often yields diminishing returns due to increased orthogonalization overhead (typically $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$8 per iteration), and must be balanced for each application (Gaul et al., 2012, Neuenhofen, 2015).
• Convergence Acceleration: Across many tested regimes (imaging, CFD, model reduction, tomography), iteration count savings of $A^{(i)}x^{(i)} = b^{(i)},\quad i = 1,2,\ldots ,$9–$A^{(i)}$0 over non-recycling methods are frequently reported, with wall-clock time or matvec reductions typically ranging from 30–80% (Ehrhardt et al., 2024, Parks et al., 2016, Gaul et al., 2012, Carlberg et al., 2015).
• Flexibility with Non-Hermitian, Indefinite, or Nonsymmetric Systems: Recycling extends to IDR, BiCGSTAB, FOM, and non-symmetric GMRES, with algorithmic modifications to preserve short recurrences and low storage (Ahuja et al., 2014, Neuenhofen, 2015, Ahuja et al., 2010).
• Stopping Criteria: In advanced contexts, surrogate error estimators (e.g., via low-rank GSVD for hypergradients) are used instead of classical norm-based stopping, enabling reliable inexactness control with respect to specific objective functions (Ehrhardt et al., 2024, Blondeau et al., 2023).

5. Applications and Numerical Validation

Recycling Krylov subspaces has been validated across diverse fields:

• Bilevel Learning and Inverse Imaging: Accelerated inner solves for bilevel hypergradient computation, reducing the iteration count and computational time by an order of magnitude in imaging applications (Ehrhardt et al., 2024).
• Shape and Topology Optimization: Application to evolving mesh PDE-constrained optimization, using FE-function based subspace transfer and warm-start Krylov-Schur refinement, achieving 20–30% savings in Krylov steps across mesh evolution cycles (Bolten et al., 2020).
• CFD and Aerostructural Solvers: Block and flexible recycling methods embedded in large-scale CFD and partitioned aerostructural adjoint solvers yield reductions of 19–39% in matrix-vector product count (Blondeau et al., 2023, Amritkar et al., 2015).
• Model Order Reduction: In IRKA and moment-matching reduction, biorthogonal recycling achieves 50–70% iteration savings and 30–50% runtime reductions (Ahuja et al., 2010, Ahuja et al., 2014).
• Matrix Functions in Lattice QCD: Augmented and sketched block FOM-type recycling yields 1–2 order-of-magnitude gains for sign, inverse, or fractional matrix exponentiation (Burke et al., 2022, Burke et al., 2023, Burke et al., 2023).
• Streaming/Dynamic Tomography: Recycled MM-GKS and related approaches enable memory-constrained computation for huge, streamed inverse problems while retaining high-quality reconstructions (Pasha et al., 2023).

6. Advanced Topics and Future Directions

• Randomized Sketching and Restarted/Block Variants: Recent methods deploy random projection (sketching) to permit recycling without full basis orthogonalization, enabling favorable complexity scaling on modern hardware and for very large $A^{(i)}$1 (Burke et al., 2023, Burke et al., 2023).
• Hierarchical and Two-Level Recycling: For sequences of nonlinear or shifted systems, hierarchical recycle spaces (e.g., principal/global and local/shift-specific) are constructed, leveraging both structural similarity and local spectral adaptation (Kilmer et al., 2023).
• Goal-Oriented and Model-Driven Subspace Selection: Optimization-directed truncation, e.g., via goal-oriented POD, aligns the recycle space with practical error metrics or quantities of interest, outperforming generic spectral approaches in reduced-order modeling, optimization, and control (Carlberg et al., 2015).
• Function Action Approximations: Closed-form Bubnov–Galerkin projections paired with Schur/SVD-based subspace recycling avoid the need for shifted linear solves in evaluating matrix functions, offering stable and fast alternatives (Burke et al., 2022, Burke et al., 2023).

The continued development of recycling Krylov methods is driven by application demands for scalability, adaptability, and robust control over inexactness, with ongoing research in randomized/sketched implementations, recycling in highly nonlinear or structural evolution regimes, and integration with learning-based or data-driven subspace selection strategies.