Adapting the Lanczos algorithm to matrices with almost continuous spectra

TL;DR

This work addresses efficient approximation of $F(s)=B^T(A+sI)^{-1}B$ for large SPD matrices with dense spectra, common in discretizations of PDEs on unbounded domains. The authors replace classical block Lanczos quadrature with a Kreĭn--Nudelman square-root terminator, implemented as an adaptive tail parameter $\phi$, to model the continuous spectral measure and suppress artificial reflections. The method yields a low-rank modification to the Lanczos matrix and, via a port-Hamiltonian interpretation, promotes spectral adaptation to scattering poles, improving both accuracy and physical fidelity (e.g., converting standing to outgoing waves in wave problems). Numerical experiments on 2D diffusion and 3D electromagnetic diffusion show substantial error reductions over Gauss and Gauss–Radau quadratures, with the KN approach also delivering qualitative improvements for wave propagation. The framework suggests a scalable path to accurate transfer-function and spectral-density approximations for operators with continuous spectra and motivates several extensions to other matrix functions and randomized techniques.

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.

Figures (9)

• Figure 1: This figure shows a contour plot of respectively $\log |1+\widehat{ {\cal F}}^{\phi}_m(s)|$ (Kreĭn--Nudelman) and $\log |1+ {\cal F}_m(s)|$ (Gauß) with $p=1$ (SISO) for $m=90$ (left) and $m=290$ (right) as a function of $\sqrt{s}$ in the complex plane for the 3D Maxwell problem from subsection \ref{['sec:EMdiff']}. To show both Riemann sheets in a single figure, the contours are shown as a function of the real and imaginary part of $\sqrt{s}$. The Gauß quadrature ${\cal F}_m(s)$ is a single-valued function, thus it is even with respect to the imaginary axis, and as a function of $\sqrt{s}$ has poles on the imaginary axis. Due to introduction of $\sqrt{s}$ in the definition of the Kreĭn--Nudelman quadrature $\widehat{ {\cal F}}^\phi_m(s)$ it becomes a two-valued function (similar to the transfer function for the continuous problem), with the branch-cut at the imaginary axis and (scattering) poles at the second Riemann sheet $\Re(\sqrt {s})<0$. The dissipation optimization moves poles away from the branch cut adaptively; i.e. poorly converged Ritz values in the dense spectral area are moved further than converged well-separated Ritz values. This distance diminishes for the larger $m$ together with the quadrature approximation error, thus the Kreĭn--Nudelman poles approach the ones of the Gauß quadrature as $m$ grows. The spectral measure at the branch cut is smooth in the area of poor convergence, collectively approximating the Stieltjes moments at dense intervals of $A$'s spectrum.
• Figure 2: Grid, heat conductivity $\sigma({{\itbf x}})$ and transducer locations of the heat diffusion testcase. In the SISO results for $p=1$, only transducer 1 is used. Only the first geometrically increasing grid steps are shown.
• Figure 3: Converge result for the 2D diffusion case in the SISO setting ($p=1$).
• Figure 4: 3D rendering of the simulated configuration with two inclusions for the tri-axial electromagnetic borehole tool, p=6.
• Figure 5: Converge results for the 3D EM case in the MIMO setting ($p=6$). As before "GR average" is arithmetic average of Gauß and Gauß--Radau quadratures.

Theorems & Definitions (8)

• Remark 1
• Remark 2
• Proposition 1
• Remark 3
• Proposition 2
• proof
• proof
• Remark 4