Rational Krylov Space Methods

• Rational Krylov space methods are projection-based algorithms that use adaptive shift-and-invert strategies to efficiently approximate functions of large operators.
• They extend polynomial Krylov techniques by incorporating user-selected poles to capture spectral features and accelerate convergence in both linear and nonlinear problems.
• Applications span solving Lyapunov equations, computing matrix exponentials, and model order reduction, highlighting their practical impact in scientific computing.

Rational Krylov space methods constitute a broad class of projection-based algorithms designed to efficiently approximate vector- and matrix-valued functions of large-scale operators, with applications spanning numerical linear algebra, scientific computing, and control theory. Extending classical polynomial Krylov subspace techniques, rational Krylov methods replace powers of an operator with sequences of shift-and-invert steps at user-selected "poles", yielding subspaces that adaptively resolve spectral features and accelerate convergence for both linear and nonlinear problems. The pole (shift) selection mechanisms and adaptable recurrence structures are central to their superior approximation properties.

1. Definitions and Theoretical Structure

A rational Krylov subspace associated with a matrix $A\in\mathbb{C}^{n\times n}$, seed vector $b\in\mathbb{C}^n$, and poles $\{\xi_1,\xi_2,\dots,\xi_m\}$ is defined as

$$\mathcal{K}_m(A, b, \{\xi_j\}) = \text{span}\left\{b,\, (A-\xi_1 I)^{-1} b,\, (A-\xi_2 I)^{-1}(A-\xi_1 I)^{-1}b,\,\ldots,\,\left(\prod_{j=1}^m (A-\xi_j I)^{-1} \right) b \right\}$$

(Kolesnikov et al., 2014, Simoncini, 2016, Kindermann et al., 2023, Daas et al., 30 Jun 2025). The rational subspace can also be written as the span of functions $r(A)b$ where $r(z)$ is a rational function constructed from polynomials in the numerator and products of $(z-\xi_j)$ or $(z+\xi_j)$ in the denominator, capturing spectral behavior relative to chosen shifts.

Block and tensor generalizations, such as block rational Krylov subspaces and tensor rational block Arnoldi/Lanczos procedures, use shift-and-invert steps on block vectors or tensors to build subspaces suitable for multidimensional or multi-parameter problems (Barkouki et al., 2024, Casulli, 2023).

The choice and distribution of poles, whether adaptive or a priori, is crucial: well-chosen poles lead to superlinear convergence in rational approximation of matrix functions, solution of shifted systems, and model order reduction (MOR) (Beckermann et al., 2020, Lin, 2021, Barkouki et al., 2023).

2. Algorithmic Frameworks and Recurrences

Rational Krylov methods generalize Arnoldi and Lanczos: each iteration computes a new vector through a shifted linear system solve, orthogonalizes it against previous basis vectors, and updates projected matrices. The rational Arnoldi process builds an orthonormal basis $V_m$, satisfying rational recurrence relations

$A V_m K_m = V_{m+1} H_m$

where $K_m$ encodes the effect of shifts and $H_m$ is upper Hessenberg (Benzi et al., 2022, Barkouki et al., 2023). For non-Hermitian or biorthogonal settings, rational Lanczos variants construct biorthonormal bases $V, W$ satisfying oblique projections and compact tridiagonal (or more general pencil) representations (Buggenhout et al., 2018). Short recurrence schemes have been developed for biorthogonal rational projections, resulting in pentadiagonal or tridiagonal pencils suitable for large non-Hermitian eigenproblems and matrix function evaluation.

Minimal residual (MR)-type rational Krylov methods enforce a Petrov-Galerkin residual orthogonality, leading to least-squares solves in the projected space, and can handle large sequences of shifted systems efficiently (Daas et al., 30 Jun 2025). Mixed polynomial-rational Krylov methods incorporate both powers and shifted-invert vectors, yielding pentadiagonal projected operators and CG-type short recurrences (Kindermann et al., 2023).

Implicit restarting, compact representations, and two-level decompositions (e.g., R-CORK) enable memory-efficient rational Krylov expansions for large-scale rational eigenvalue problems, reducing both orthogonalization and solve costs through hierarchical block compression (Dopico et al., 2017).

3. Pole Selection, Adaptivity, and Residual Monitoring

Strategic pole selection—fixed, adaptive, or via greedy optimization—is central to rational Krylov convergence. Techniques include:

• Adaptive shift selection via residual monitoring: selecting the next pole as the shift corresponding to the largest current residual, thereby annihilating the worst-approximated direction in subsequent steps (Daas et al., 30 Jun 2025, Bergermann et al., 2023).
• Moments-matching and interpolation-based greedy selection: maximizing factors controlling the interpolation error (e.g., $$\Bigl|\prod_{j=1}^m(s-\sigma_j)/\prod_{j=1}^m(s-\lambda_j)\Bigr|$$) (Lin, 2021).
• Rational minimax or Zolotarev pole sequences: employing extremal distributions for best uniform approximation of target functions over spectral intervals (Benzi et al., 2022, Beckermann et al., 2020).
• Adaptive quadrature, field-of-values or pole selection based on the spectral properties of the projected or closed-loop matrices in Riccati/Lyapunov equations (Simoncini, 2016, Bertram et al., 2023, Kolesnikov et al., 2014).

Pseudocode in key works demonstrates flexible shift-picking, residual-based stopping, and performance guarantees under various pole-selection strategies.

4. Applications and Advanced Extensions

Rational Krylov methods have transformed algorithms for:

• Large-scale Lyapunov and algebraic Riccati equations: providing low-rank factorizations via rational subspace projection, adaptive shifts, and small dense solves (Kolesnikov et al., 2014, Simoncini, 2016, Bertram et al., 2023, Kürschner et al., 2018).
• Matrix exponential and Cauchy–Stieltjes function computations: block rational Lanczos/Arnoldi algorithms approximate $e^{tA} B$ and related spectral functions with scalable error bounds (Benzi et al., 2022, Barkouki et al., 2023, Bergermann et al., 2023).
• Linear inverse ill-posed problems: aggregation and rational-CG methods interpret regularized solvers as rational Krylov projections, yielding order-optimal regularization under discrepancy principles (Kindermann, 15 Jan 2026, Kindermann et al., 2023).
• Model order reduction and transfer function approximation: moment-matching and explicit remainder estimation via rational Krylov subspaces, supporting $H_\infty$-optimal reduction and adaptive greedy pole selection (Lin, 2021, Benzi et al., 2022).
• Multidimensional/tensor dynamical systems: tensor rational Krylov procedures enable low-rank solution of multidimensional Lyapunov and Sylvester equations, as well as balanced truncation for high-dimensional system reduction (Barkouki et al., 2024, Casulli, 2023).
• Rational eigenvalue problems: structure-preserving rational linearization (e.g., Even-IRA, R-CORK) facilitates efficient eigen-computations for polynomial/rational matrix functions (Dopico et al., 2017, Benner et al., 2020).

Sophisticated residual analysis, error bounds, and convergence diagnostics have been developed for all these scenarios, guaranteeing small errors with low-dimensional subspaces for well-chosen poles.

5. Error Analysis, Optimality, and Regularization Theory

Error bounds for rational Krylov approximations are typically expressed as best uniform rational approximation errors over the spectrum of the operator, depending on the chosen poles (Benzi et al., 2022, Lin, 2021, Beckermann et al., 2020). Explicit remainder formulas involve both the pole sequence and the spectrum (or Ritz values) of the projected operator, revealing systematic interpolation and quadrature connections. For matrix functions, optimality in various operator-weighted norms is precisely characterized, and rational Krylov methods can achieve near-minimax error rates in practice.

For ill-posed problems in Hilbert spaces, rational Krylov sequences—when combined with discrepancy-principle stopping rules and sufficiently large shift parameters—form order-optimal regularization schemes, matching the rates of Tikhonov and Landweber methods, but with accelerated practical convergence (Kindermann, 15 Jan 2026, Kindermann et al., 2023).

6. Computational Complexity, Storage, and Scalability

The per-iteration costs of rational Krylov methods mirror their polynomial analogs except that (a) each iteration requires the solution to a shifted linear system, and (b) projected matrices often have more complex (block Hessenberg, tridiagonal, or pentadiagonal) structure (Kolesnikov et al., 2014, Benzi et al., 2022, Kindermann et al., 2023). When adaptive shift reuse and preconditioning (e.g., AMG, incomplete factorization) are employed, rational Krylov methods can scale linearly (or near-linearly) in $n$—a major advantage for very large problems (Bergermann et al., 2023).

Memory requirements are reduced by block and compact representations, adaptive truncation, and hierarchically structured orthogonalization, with best practices established for tensor and block variants (Barkouki et al., 2024, Dopico et al., 2017, Buggenhout et al., 2018).

Numerical experiments consistently demonstrate rapid convergence—often achievable with subspace dimensions an order of magnitude smaller than polynomial methods—for properly placed poles and adaptive error control.

7. Limitations, Open Problems, and Future Directions

Key theoretical and practical challenges remain in rational Krylov methods:

• Optimal shift selection for very general spectral distributions, closed-loop operators (Riccati), and high-dimensional tensor products requires further research.
• Stability and biorthogonality loss in non-Hermitian and block algorithms necessitate robust reorthogonalization or SVD-based correction (Barkouki et al., 2024, Buggenhout et al., 2018).
• Systematic regularization theory for mixed rational–polynomial spaces, and adaptive regularization parameter selection, especially in data-driven and nonlinear problems, is an active area (Kindermann, 15 Jan 2026).
• Extending explicit error bounds, remainder formulas, and greedy pole selection strategies to MIMO systems, high-order Markov functions, and fully tensorized polynomial spaces is ongoing (Lin, 2021, Casulli, 2023).
• Efficient preconditioning and solution strategies for large-scale shift-and-invert problems—especially for tensor, block, and non-Hermitian settings—require further algorithmic refinement and analysis.

These limitations suggest fertile directions for both theoretical advancement and practical algorithm development in the rational Krylov framework, as evidenced by the ongoing research documented in recent literature.