Lanczos quadrature in $A = \tilde{I} D$ can significantly underestimate the scattering intensity of low-energy modes. This is illustrated by some examples at #339.
To fix this, we may need an approximation scheme that effectively works in the extended Krylov space,
$$\mathcal{K}_{l,m} = \mathrm{span} \{A^{-(l-1)} v, \dots, A^{-1} v, v, A v, \dots, A^{m-1} v \}$$.
The mat-vec multiplies involving $A^{-1} = D^{-1} \tilde{I}$ can become efficient after paying a one-time cost to build an LDL decomposition of $D$.
There exist recurrences to build the extended Krylov space, e.g. Jagels, Jbilou, and Reichel, J. Comput. Appl. Math. 381, 113027 (2021). But to avoid the need for an expensive reorthogonalization step, I would be interested to try a simplified scheme that decomposes,
$$f(A) = f_0(A) + f_\infty(A).$$
The parts $f_0$ and $f_\infty$ would be amenable to Lanczos quadrature using $A^{-1}$ and $A$, respectively.
Define,
$$
\begin{aligned}
f_0(x) &= \omega(x / \sigma) f(x) \\
f_\infty(x) &= [1 - \omega(x / \sigma)] f(x),
\end{aligned}
$$
involving the localized mask function
$$
\omega(x) = \frac{1}{1 + x^{p}},
$$
at some energy scale $\sigma$ and exponent $p$ (an even integer).
The idea is that $f_\infty(A)$ is amenable to approximation in the original Krylov space, since the $1 - \omega(x/\sigma)$ prefactor masks the problematic behavior at $x \sim 0$.
The remaining term is naturally expressed in the new variable $y = 1 / x$. Define also $g(y) = f(1/y)$. Then,
$$f_0(x) = \omega(1 / \sigma y) g(y) = [1 - \omega(\sigma y)] g(y).$$
In this form, we see that the $1 - \omega(\sigma y)$ prefactor masks the problematic behavior at $y \sim 0$ (equivalently, $x \sim \pm \infty$) for polynomial approximation in $y$. Therefore $f_0(x)$ is amenable to approximation in the Krylov space built from the inverse matrix, i.e., $1 / A$.
Together, these two approximations can be understood as a single approximation to $f(A)$ in the extended Krylov space.
The exponent $p$ should be sufficiently large to damp singular behavior at the origin, but otherwise should be small for smoothness. The $f(x)$ of interest behaves like $\Theta(x) / (\beta x)$ at $x \sim 0$, so a good exponent might be $p = 4$. A symmetric choice of cross-over scale would be $\sigma = \sqrt{\lambda_\min \lambda_\max}$, i.e., the geometric mean of the extreme eigenvalues of $A$. But $\sigma$ could be reduced if one wants to use fewer inverse-Krylov steps due to the expense of applying $A^{-1}$.
Perhaps a practical strategy is to allocate a fixed wall clock budget to running block-Lanczos on $A$ and separately $A^{-1}$. Then one could dynamically optimize $\sigma$ to get the best overall approximation to $f(A)$. Conjecture: Is there some sense in which this strategy is locally optimal within the restricted approximation family? Another thought is that either Lanczos recurrence could be dynamically extended if it is identified that the Ritz eigenpairs are not well converged.
Lanczos quadrature in$A = \tilde{I} D$ can significantly underestimate the scattering intensity of low-energy modes. This is illustrated by some examples at #339.
To fix this, we may need an approximation scheme that effectively works in the extended Krylov space,
The mat-vec multiplies involving$A^{-1} = D^{-1} \tilde{I}$ can become efficient after paying a one-time cost to build an LDL decomposition of $D$ .
There exist recurrences to build the extended Krylov space, e.g. Jagels, Jbilou, and Reichel, J. Comput. Appl. Math. 381, 113027 (2021). But to avoid the need for an expensive reorthogonalization step, I would be interested to try a simplified scheme that decomposes,
The parts$f_0$ and $f_\infty$ would be amenable to Lanczos quadrature using $A^{-1}$ and $A$ , respectively.
Define,
involving the localized mask function
at some energy scale$\sigma$ and exponent $p$ (an even integer).
The idea is that$f_\infty(A)$ is amenable to approximation in the original Krylov space, since the $1 - \omega(x/\sigma)$ prefactor masks the problematic behavior at $x \sim 0$ .
The remaining term is naturally expressed in the new variable$y = 1 / x$ . Define also $g(y) = f(1/y)$ . Then,
In this form, we see that the$1 - \omega(\sigma y)$ prefactor masks the problematic behavior at $y \sim 0$ (equivalently, $x \sim \pm \infty$ ) for polynomial approximation in $y$ . Therefore $f_0(x)$ is amenable to approximation in the Krylov space built from the inverse matrix, i.e., $1 / A$ .
Together, these two approximations can be understood as a single approximation to$f(A)$ in the extended Krylov space.
The exponent$p$ should be sufficiently large to damp singular behavior at the origin, but otherwise should be small for smoothness. The $f(x)$ of interest behaves like $\Theta(x) / (\beta x)$ at $x \sim 0$ , so a good exponent might be $p = 4$ . A symmetric choice of cross-over scale would be $\sigma = \sqrt{\lambda_\min \lambda_\max}$ , i.e., the geometric mean of the extreme eigenvalues of $A$ . But $\sigma$ could be reduced if one wants to use fewer inverse-Krylov steps due to the expense of applying $A^{-1}$ .
Perhaps a practical strategy is to allocate a fixed wall clock budget to running block-Lanczos on$A$ and separately $A^{-1}$ . Then one could dynamically optimize $\sigma$ to get the best overall approximation to $f(A)$ . Conjecture: Is there some sense in which this strategy is locally optimal within the restricted approximation family? Another thought is that either Lanczos recurrence could be dynamically extended if it is identified that the Ritz eigenpairs are not well converged.