Replace KPM with Lanczos

[kbarros]

The Lanczos algorithm promises these benefits relative to KPM:

• Higher accuracy at a fixed number of iterations (i.e., at a fixed Krylov subspace dimension).
• Does not require an initialization phase that bounds the spectrum prior to polynomial expansion.
• Ability to detect instabilities (negative eigenvalues) in the dynamical matrix.
• Allows for a sharp cutoff between positive and negative eigenvalues near zero. Whereas KPM is susceptible to a loss of intensity near Goldstone modes, Lanczos seems to avoid this problem. (See below.)

As currently implemented, Lanczos is currently slower than KPM per iteration. But many fewer iterations are required, so Lanczos still seems favorable in overall speed. And future optimizations may be possible, e.g., operating on multiple right-hand side vectors in parallel.

A possible disadvantage of Lanczos is the accumulation of numerical roundoff error (loss of orthogonality in basis vectors for the Krylov space), and the effects of this are harder to quantify. See https://arxiv.org/abs/2410.11090 for a modern mathematical review of the Lanczos method. A backward stability analysis, in Sec. 4, suggests that the instabilities of Lanczos are not catastrophic for the desired observables. Empirically, however, Lanczos can still exhibit severe artifacts in estimating the intensities of bands near zero energy (quantified in examples below). Note, however, that these are cases where KPM also fails.

Another consideration is that Lanczos requires direct diagonalization of a tridiagonal matrix. After $k$ Lanczos iterations, this cost should scale like $\mathcal{O}(k^2)$, which may dominate if $k$ sufficiently large.

[kbarros]

Here is an example where floating point roundoff leads to a large error in the intensity. This happens very small energy modes lead to ill-conditioning.

cryst = Sunny.bcc_crystal()
sys = System(cryst, [1 => Moment(; s=1, g=2)], :dipole; dims=(1, 1, 1))
sys = to_inhomogeneous(sys)
set_exchange_at!(sys, -1.0, (1, 1, 1, 1), (1, 1, 1, 2); offset=[0, 0, 0]) 

polarize_spins!(sys, [0, 0, 1])
energy_per_site(sys)

kernel = lorentzian(fwhm=0.2)
energies = range(0, 6, 1000)
qs = [[0.35, 0.35, 0]]

measure = ssf_custom((q, ssf) -> real(ssf[1, 1]), sys) # Sˣˣ(q,ω) only

swt = SpinWaveTheoryKPM(sys; measure, lanczos_iters=4, tol=1, regularization=1e-8)
res = intensities(swt, qs; energies, kernel)

swt0 = SpinWaveTheory(sys; measure, regularization=1e-8)
res0 = intensities(swt0, qs; energies, kernel)

lines(vec(res.data))
lines!(vec(res0.data))

Note that the blue curve is missing some intensity near zero energy.

In this particular case, the magnetic cell consists of two sites and a single bond. The dynamical matrix has size 4×4. The two modes should energies of 2 and 0, slightly modified by the default regularization of $ϵ = 10^{-8}$,

D = Sunny.dynamical_matrix(swt.swt, Sunny.Vec3(only(qs)))
It = Diagonal([1, 1, -1, -1])
eigvals(It * D) # ±(2 + ϵ), ±ϵ

However, there is an inversion ("particle-hole") symmetry in this problem that Lanczos doesn't exactly respect. Specifically, because $D_q = D_{-q}$, and similarly for the right-hand side vector representing the observable, all the diagonal elements $(α_1, ..., α_k)$ of the Lanczos tridiagonal matrix should be exactly zero. In practice, however, the Lanczos calculated tridiagonal matrix has significant deviation from 0 in its bottom-right matrix element, $α_4 = T_{4,4}$. Having broken this symmetry, the diagonalization of $T$ allows significant "mixing" between eigenvectors in the negative and positive parts of the spectrum near zero. In practice, this mixing causes huge error in the estimated intensities near zero energy.

A workaround for this specific problem could be to modify the Lanczos iterations to exactly respect inversion symmetry, e.g., following the algorithm in Fig. 4 of Grüning et al, Implementation and testing..., arXiv:1102.3909. However, other applications of LSWT will not have this exact inversion symmetry, and a general solution is lacking.

One idea is to double the dimension of the matrix size by working with both wavevectors $\pm q$ simultaneously. Then, in principle, the symmetrized Lanczos version would always apply, and the positive and negative parts of the spectrum become redundant. However, this still leaves the possibility of mixing between the eigenvectors of the positive part of the spectrum, and it's not clear there is a win.

[kbarros]

Here is an example where Lanczos misses intensity even though there is not significant floating point roundoff. This example builds a 3x3 system, nearest neighbor ferromagnetic interactions, where 40% sites are replaced by vacancies.

import Random

function create_dilute_sys(n, x; s=1, mode=:dipole)
    a = b = 1.0
    c = 10.0
    latvecs = lattice_vectors(a, b, c, 90, 90, 90)
    positions = [[0, 0, 0]]
    cryst = Crystal(latvecs, positions)
    J₁ = 1.0
    sys = System(cryst, [1 => Moment(; s, g=2)], mode; dims=(n, n, 1))
    set_exchange!(sys, J₁, Bond(1, 1, [1, 0, 0]))

    sys_inhom = to_inhomogeneous(sys)
    nvacancies = Int(round(n*n*x))

    vacancy_sites = sample(eachsite(sys_inhom), nvacancies; replace=false)
    for site in vacancy_sites
        Sunny.set_vacancy_at!(sys_inhom, site)
    end

    randomize_spins!(sys_inhom) 
    minimize_energy!(sys_inhom; maxiters=10_000)
    return sys_inhom, cryst
end

kernel = lorentzian(fwhm=0.2)
energies = range(0, 6, 1000)
qs = [[0.35, 0.35, 0]]

Random.seed!(2)
n = 3
x = 0.4

sys, cryst = create_dilute_sys(n, x)

measure = ssf_custom((q, ssf) -> real(ssf[1, 1]), sys)

swt = SpinWaveTheoryKPM(sys; measure, lanczos_iters=2, tol=1, regularization=1e-8)
res_Lanczos = intensities(swt, qs; energies, kernel)

swt_lswt = SpinWaveTheory(sys; measure, regularization=1e-8)
res_LSWT = intensities(swt_lswt, qs; energies, kernel)

lines(energies, vec(res_Lanczos.data), label="Lanczos")
lines!(energies, vec(res_LSWT.data), label="Reference")
axislegend()

To emphasize the failure, just two Lanczos iterations are being executed. The is, the Krylov subspace has dimension 2. This yields a single estimated excitation band with energy 0.9. It is not clear whether there is any flexibility to modify the placement of this band. Because the estimated energy significantly deviates from 0, a significant amount of intensity is missing.

To verify that this is not an issue with FP roundoff, note that increasing to lanczos_iters=4 yields quite good agreement between Lanczos and direct diagonalization.

Note, also, strong sensitivity to the choice of regularization. Reducing from the default of 1e-8 down to regularization=1e-2 shows instead this behavior:

Here it is clear Lanczos tries to "split the difference" between the two bands with intensity, although it is apparently using a weighted average that deemphasizes the divergence near zero energy.