A Krylov projection algorithm for large symmetric matrices with dense spectra

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{\ss} and Gau\ss -Radau quadratures computed using the block-Lanczos method significantly reduces approximation errors for such problems. Here, we introduce an adaptive Kre\u{i}n-Nudelman extension to the (block) Lanczos recursions, allowing further acceleration at negligible $o(n)$ cost. Similar to the Gau\ss -Radau quadrature, a low-rank modification is applied to the (block) Lanczos matrix. However, unlike the Gau\ss -Radau quadrature, this modification depends on $\sqrt{s}$ and can be considered in the framework of the Hermite-Pad\'e 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.