Implement Lanczos-based `spectrum` method

[albertomercurio]

For a given steady-state correlation function $\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle$, spectrum calculates the quantity

S(\omega) = \int_{-\infty}^\infty \left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle e^{-i \omega \tau} d \tau

— see also the documentation on Emission spectrum.

Compared to spectrum_correlation_fft, spectrum works directly in the frequency domain by using a SpectrumSolver. In fact, in the vectorized representation, $S(\omega)$ can be expressed as

S(\omega) = 2 \mathrm{Im}
\bigg\langle I \bigg\vert \hat{A}
\frac{1}{\omega + i\mathcal{L}}
\hat{B} \bigg\vert \rho_{\mathrm{ss}} \bigg\rangle

Currently, there two such solvers:

• ExponentialSeries, which employes a full spectral decomposition;
• PseudoInverse, which computes the propagator $(\omega + i\mathcal{L})^{-1}$ as the solution of a sparse linear system at each frequency.

PseudoInverse is currently the fastest option for large systems, but it requires (pseudo-)inverting $\omega + i\mathcal{L}$ at each frequency value.

Another possible avenue for the computation of the spectrum would be to employ the well-known Lanczos method, which is currently a gold standard in condensed matter calculations. The method is extensively discussed in [1] for systems governed by a Hamiltonian $H$, in which response functions are calculated on top of the ground state $\vert \psi_{0} \rangle$:

\begin{equation}
    S(\omega)
    = 2  \mathrm{Im} \Braket{\psi_0 | \hat{A}\frac{1}{\omega-\hat{H}}\hat{B} | \psi_0}
\end{equation}

In a nutshell, for $\hat{A} = \hat{B}^{\dagger}$, $H$ is iteratively expressed in a Krylov basis built on top of $\vert v_{1} \rangle = \hat{B} \vert \psi_{0} \rangle$ by successive applications of $H$. In this basis the Hamiltonian becomes tridiagonal:

\begin{equation}
    \hat{H}_{\mathcal{K}(\ket{v_1})} =
    \begin{pmatrix}
        a_1    & b_2 & 0   & 0   & \cdots \\
        b_2    & a_2 & b_3 & 0   &        \\
        0      & b_3 & a_3 & b_4 &        \\
        0      & 0   & b_4 & a_4 &        \\
        \vdots &     &     &     & \ddots
    \end{pmatrix}\,.
\end{equation}

Therefore, the desidered expectation value of the propagator $(\omega - H)^{-1}$ can be read off its top-left matrix element $\left( \omega - \hat{H} \right)^{-1}_{11}$.

Thanks to the tridiagonal structure of the propagator, the top-left matrix element can be recursively expressed as a continuous fraction for each Krylov iteration by using a suitable expression for the inverse of a $2 \times 2$ block matrix:

\begin{equation}
    \left( \omega - \hat{H} \right)^{-1}_{11} 
    = \cfrac{1}{\omega-a_1-\cfrac{b_2^2}{\omega-a_2-\cfrac{b_3^2}{\omega-a_3-\cdots}}} \,.
\end{equation}

The algorithm stops when the continuous fraction reaches the desired precision or the maximum number of iterations. Early stops are possible in case of loss of orthogonality, a known drawback of the Lanczos algorithm; however, they are not as catastrophic as for the case of pure diagonalization, they can be easily detected as $b_n$ becomes $0$, and they just result in a termination of the continued fraction.

Cases in which $\hat{A} \neq \hat{B}^{\dagger}$ can be handled by starting the Krylov basis from $(\hat{B} \pm \hat{A}^{\dagger}) \vert \psi_{0} \rangle$ and expanding the sums in the expectation value.

In our case, the goal would be to implement a SpectrumSolver for spectrum called Lanczos, which implements a non-symmetric variant of the algorithm due to the fact that $\mathcal{L}$, unlike $\mathcal{H}$, is non-Hermitian. It should support a usage like:

spectrum(H, ωlist, c_ops, A, B; solver = Lanczos())

and provide, at a minimum, a tol parameter and a maxiter parameter, respectively controlling the desired precision and the maximum number of Krylov iterations.

Moreover, this new implementation should be added to the benchmarks, in order to track its performances together with the other solvers.

---

1. Koch, E. (2011). The Lanczos Method. In E. Pavarini, E. Koch, D. Vollhardt, & A. Lichtenstein (Eds.), The LDA+DMFT approach to strongly correlated materials (Vol. 1, pp. 8.1-8.30). Forschungszentrum Jülich. ↩