Adapting the Lanczos algorithm to matrices with almost continuous spectra

TL;DR

This paper adapts the Lanczos algorithm to better approximate matrices with dense, continuous spectra, improving efficiency and accuracy in large-scale PDE discretizations, especially for wave propagation problems.

Contribution

It introduces a quadratic terminator using Kren--Nudelman semi-infinite strings to model the continuous spectral measure more effectively than traditional Krylov methods.

Findings

Significant error reductions in PDE discretizations

Enhanced computation of MIMO transfer functions

Improved wave propagation simulations

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\u{i}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\u{i}n--Nudelman spectrum, we obtain a qualitative improvement in finite-difference approximations, effectively transforming standing-wave artifacts into outgoing propagating waves.