Adapting the Lanczos algorithm to matrices with almost continuous spectra

Abstract: We consider the approximation of $B^T (A+sI)^{-1} B$ where $A\in\mathbb{R}^{n\times n}$ is large, symmetric positive definite, and has a dense spectrum, and $B\in\mathbb{R}^{n\times p}$, $p\ll n$. Our target application is the computation of Multiple-Input Multiple-Output transfer functions arising from large-scale discretizations of problems with continuous spectral measures, such as linear time-invariant PDEs on unbounded domains. Traditional Krylov methods, such as Lanczos or conjugate gradients, focus on resolving individual eigenvalues of a dense discretization, while ignoring the underlying continuous spectral measure that these points approximate. We argue that it is more efficient to model the inherent branch cut of the original operator than to exhaustively resolve the artificial point spectrum induced by discretization. We place this problem in a framework, known in the physics literature as the square-root terminator. To overcome its limitations, we formulate a quadratic terminator using Kreĭn--Nudelman semi-infinite strings, with parameters chosen adaptively to maximize relative energy outflow. This approach results in a low-rank modification to the (block) Lanczos matrix, dependent on $\sqrt{s}$, with an additional $O(n)$ cost. We demonstrate significant error reductions for large-scale self-adjoint PDE discretizations in unbounded domains, including two- and three-dimensional Maxwell's equations in diffusive regimes. The method proves particularly advantageous in computing state-space solutions for wave propagation, specifically for 2D wave and 3D Maxwell's operators. Implicitly replacing the conventional Lanczos spectral decomposition with a representation in terms of the continuous Kreĭn--Nudelman spectrum, we obtain a qualitative improvement in finite-difference approximations, effectively transforming standing-wave artifacts into outgoing propagating waves.