A Krylov projection algorithm for large symmetric matrices with dense spectra

Abstract

We consider the approximation of $B^T(A+sI{)}^{-1}B$ for large s.p.d. $A\in {\mathbb{R}}^{n\times n}$ with dense spectrum and $B\in {\mathbb{R}}^{n\times p}$, $p\ll n$. We target the computations of Multiple-Input Multiple-Output (MIMO) transfer functions for large-scale discretizations of problems with continuous spectral measures, such as linear time-invariant (LTI) PDEs on unbounded domains. Traditional Krylov methods, such as the Lanczos or CG algorithm, are known to be optimal for the computation of $(A+sI{)}^{-1}B$ with real positive $s$, resulting in an adaptation to the distinctively discrete and nonuniform spectra. However, the adaptation is damped for matrices with dense spectra. It was demonstrated in [Zimmerling, Druskin, Simoncini, Journal of Scientific Computing 103(1), 5 (2025)] that averaging Gauß and Gauß-Radau quadratures computed using the block-Lanczos method significantly reduces approximation errors for such problems. Here, we introduce an adaptive Kreĭn-Nudelman extension to the (block) Lanczos recursions, allowing further acceleration at negligible $o(n)$ cost. Similar to the Gauß-Radau quadrature, a low-rank modification is applied to the (block) Lanczos matrix. However, unlike the Gauß-Radau quadrature, this modification depends on $\sqrt{s}$ and can be considered in the framework of the Hermite-Padé approximants, which are known to be efficient for problems with branch-cuts, that can be good approximations to dense spectral intervals. Numerical results for large-scale discretizations of heat-diffusion and quasi-magnetostatic Maxwell's operators in unbounded domains confirm the efficiency of the proposed approach.

Block Lanczos, Quadrature, Transfer function, Kreĭn-Nudelman, Hermite-Padé