A Krylov projection algorithm for large symmetric matrices with dense spectra
Vladimir Druskin, Jörn Zimmerling
TL;DR
This work tackles the efficient computation of the MIMO transfer function $F(s)=B^T(A+sI)^{-1}B$ for large, dense-spectrum SPD matrices arising from unbounded-domain PDE discretizations. It introduces a adaptive Kreĭn-Nudelman extension to block Lanczos recursions, adding absorbing-boundary-like modifications parameterized by $(oldsymbol heta,oldsymbol ealfn{ heta})$, and connects these to Hermite-Padé approximants to better handle branch cuts in the spectral measure. The authors derive two-sided bounds, develop optimization strategies (reflection matching and targeted smoothing of the spectral measure), and provide a practical framework to choose KN parameters that accelerate convergence beyond averaged Gauß/Gauß-Radau quadratures. Numerical experiments on 2D diffusion and 3D quasi-magnetostatic Maxwell problems demonstrate substantial acceleration and robustness of the KN approach across real and imaginary shifts, validating its potential for large-scale MIMO transfer computations in unbounded-domain PDE contexts. The work also establishes a theoretical bridge between Stieltjes-string representations, continued fractions, and Hermite-Padé approximants, suggesting promising extensions to more general matrix functions and parametric model reduction.
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.