Working!

examples/09_Disorder_KPM.jl

@@ -47,51 +47,46 @@ plot_intensities(res)
 
 sys_inhom = to_inhomogeneous(repeat_periodically(sys, (10, 10, 1)))
 
-# Use [`symmetry_equivalent_bonds`](@ref) to iterate over all nearest neighbor
-# bonds of the inhomogeneous system. Modify each AFM exchange with a noise term
-# that has variance of 1/3. The newly minimized energy configuration allows for
-# long wavelength modulations on top of the original 120° order.
+# Iterate over all [`symmetry_equivalent_bonds`](@ref) for nearest neighbor
+# sites in the inhomogeneous system. Set the AFM exchange to include a noise
+# term that has a standard deviation of 1/3. The newly minimized energy
+# configuration allows for long wavelength modulations on top of the original
+# 120° order.
 
-for (site1, site2, offset) in symmetry_equivalent_bonds(sys_inhom, Bond(1,1,[1,0,0]))
+for (site1, site2, offset) in symmetry_equivalent_bonds(sys_inhom, Bond(1, 1, [1, 0, 0]))
     noise = randn()/3
     set_exchange_at!(sys_inhom, 1.0 + noise, site1, site2; offset)
 end
 
 minimize_energy!(sys_inhom, maxiters=5_000)
 plot_spins(sys_inhom; color=[S[3] for S in sys_inhom.dipoles], ndims=2)
 
-# Traditional spin wave theory calculations become impractical for large system
-# sizes. Significant acceleration is possible with the [kernel polynomial
-# method](https://arxiv.org/abs/2312.08349). Enable it by selecting
-# [`SpinWaveTheoryKPM`](@ref) in place of the traditional
-# [`SpinWaveTheory`](@ref). Using KPM, the cost of an [`intensities`](@ref)
-# calculation becomes linear in system size and scales inversely with the width
-# of the line broadening `kernel`. Error tolerance is controlled through the
-# dimensionless `tol` parameter. A relatively small value, `tol = 0.01`, helps
-# to resolve the large intensities near the ordering wavevector. The alternative
-# choice `tol = 0.1` would be twice faster, but would introduce significant
-# numerical artifacts.
+# Spin wave theory with exact diagonalization becomes impractical for large
+# system sizes. Significant acceleration is possible with an iterative Krylov
+# space solver. With [`SpinWaveTheoryKPM`](@ref), the cost of an
+# [`intensities`](@ref) calculation scales linearly in the system size and
+# inversely with the width of the line broadening `kernel`. Good choices for the
+# dimensionless error tolerance are `tol=0.05` (more speed) or `tol=0.01` (more
+# accuracy).
 #
 # Observe from the KPM calculation that disorder in the nearest-neighbor
 # exchange serves to broaden the discrete excitation bands into a continuum.
 
-swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.01)
+swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.05)
 res = intensities(swt, path; energies, kernel)
 plot_intensities(res)
 
 # Now apply a magnetic field of magnitude 7.5 (energy units) along the global
-# ``ẑ`` axis. This field fully polarizes the spins. Because gap opens, a larger
-# tolerance of `tol = 0.1` can be used to accelerate the KPM calculation without
-# sacrificing much accuracy. The resulting spin wave spectrum shows a sharp mode
-# at the Γ-point (zone center) that broadens into a continuum along the K and M
-# points (zone boundary).
+# ``ẑ`` axis. This field fully polarizes the spins and opens a gap. The spin
+# wave spectrum shows a sharp mode at the Γ-point (zone center) that broadens
+# into a continuum along the K and M points (zone boundary).
 
 set_field!(sys_inhom, [0, 0, 7.5])
 randomize_spins!(sys_inhom)
 minimize_energy!(sys_inhom)
 
 energies = range(0.0, 9.0, 150)
-swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.1)
+swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.05)
 res = intensities(swt, path; energies, kernel)
 plot_intensities(res)
 
@@ -105,8 +100,10 @@ for site in eachsite(sys_inhom)
     noise = randn()/6
     sys_inhom.gs[site] = [1 0 0; 0 1 0; 0 0 1+noise]
 end
+randomize_spins!(sys_inhom)
+minimize_energy!(sys_inhom)
 
-swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.1)
+swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.05)
 res = intensities(swt, path; energies, kernel)
 plot_intensities(res)
 

src/KPM/Lanczos.jl

@@ -1,4 +1,26 @@
-# Reference implementation, without consideration of speed.
+# Reference implementation of the generalized Lanczos method, without
+# consideration of speed. The input A can be any Hermitian matrix. The optional
+# argument S is an inner-product that must be both Hermitian and positive
+# definite. Intuitively, Lanczos finds a low-rank approximant A S ≈ Q T Q† S,
+# where T is tridiagonal and Q is orthogonal in the sense that Q† S Q = I. The
+# first column of Q is the initial vector v, which must be normalized. The
+# Lanczos method is useful in two ways. First, eigenvalues of T are
+# representative of extremal eigenvalues of A (in an appropriate S-measure
+# sense). Second, one obtains a very good approximation to the matrix-vector
+# product, f(A S) v ≈ Q f(T) e₁, valid for any function f [1].
+# 
+# The generalization of Lanczos to an arbitrary measure S is implemented as
+# follows: Each matrix-vector product A v is replaced by A S v, and each vector
+# inner product w† v is replaced by w† S v. Similar generalizations of Lanczos
+# have been considered in [2] and [3].
+#
+# 1. N. Amsel, T. Chen, A. Greenbaum, C. Musco, C. Musco, Near-optimal
+#    approximation of matrix functions by the Lanczos method, (2013)
+#    https://arxiv.org/abs/2303.03358v2.
+# 2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),
+#    https://doi.org/10.1090/s0025-5718-1982-0669648-0
+# 3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),
+#    https://doi.org/10.1016/j.commatsci.2011.02.021.
 function lanczos_ref(A, S, v; niters)
     αs = Float64[]
     βs = Float64[]
@@ -33,35 +55,40 @@ function lanczos_ref(A, S, v; niters)
 end
 
 
-# Apply Lanczos iterations to a Hermitian matrix A, with the option to specify
-# an inner-product measure S that must be both Hermitian and positive definite.
-# Loosely speaking, Lanczos finds a low-rank approximant A S ≈ Q T Q† S, where T
-# is tridiagonal and Q is orthogonal in the measure S, i.e. Q† S Q = I. The
-# first column of Q is the initial vector v. The approximant is useful purposes
-# of multiplying by v. For example, f(A S) v ≈ Q f(T) e₁ is a near-optimal
-# approximation within the Krylov space, for any function f [1]
+# An optimized implementation of the generalized Lanczos method. See
+# `lanczos_ref` for explanation, and a reference implementation. This optimized
+# implementation follows "the most stable" of the four variants considered by
+# Paige [1], as listed on Wikipedia. Here, however, we allow for an arbitary
+# measure S. In practice, this means that each matrix-vector product A v is
+# replaced by A S v, and each inner product w† v is replaced by w† S v.
 #
-# Returns T, as well as Q† lhs, if left-hand side vectors are specified.
+# In this implementation, each Lanczos iteration requires only one matrix-vector
+# multiplication involving A and S, respectively. 
 #
-# The implementation below follows the notation of Wikipedia, and is the most
-# stable of the four variants considered by Paige [2].
+# Details:
 #
-# The generalization to non-identity matrices S is implemented as follows: Each
-# matrix-vector product A v becomes A S v, and each vector inner product w† v
-# becomes w† S v. Each Lanczos iteration requires only one matrix-vector
-# multiplication for A and S, respectively. Similar generalizations of Lanczos
-# have been considered in [3] and [4].
+# * The return value is a pair, (T, lhs† Q). The symbol `lhs` denotes "left-hand
+#   side" column vectors, packed horizontally into a matrix. If the `lhs`
+#   argument is ommitted, its default value will be a matrix of width 0.
+# * If the initial vector `v` is not normalized, then this normalization will be
+#   performed automatically, and the scale `√v S v` will be absorbed into the
+#   return value `lhs† Q`.
+# * The number of iterations will never exceed length(v). If this limit is hit
+#   then, mathematically, the Krylov subspace should be complete, and the matrix
+#   T should be exactly similar to the matrix A S. In practice, however,
+#   numerical error at finite precision may be severe. Further testing is
+#   required to determine whether the Lanczos method is appropriate in this
+#   case, where the matrices have modest dimension. (Direct diagonalization may
+#   be preferable.)
+# * After `min_iters` iterations, this procedure will estimate the spectral
+#   bandwidth Δϵ between extreme eigenvalues of T. Then `Δϵ / resolution` will
+#   be an upper bound on the total number of iterations (which includes the
+#   initial `min_iters` iterations as well as subsequent ones).
 #
-# 1. N. Amsel, T. Chen, A. Greenbaum, C. Musco, C. Musco, Near-optimal
-#    approximation of matrix functions by the Lanczos method, (2013)
-#    https://arxiv.org/abs/2303.03358v2.
-# 2. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
+# 1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
 #    https://doi.org/10.1093%2Fimamat%2F10.3.373.
-# 3. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),
-#    https://doi.org/10.1090/s0025-5718-1982-0669648-0
-# 4. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),
-#    https://doi.org/10.1016/j.commatsci.2011.02.021.
-function lanczos(mulA!, mulS!, v; niters, lhs=zeros(length(v), 0))
+
+function lanczos(mulA!, mulS!, v; min_iters, resolution=Inf, lhs=zeros(length(v), 0), verbose=false)
     αs = Float64[]
     βs = Float64[]
     lhs_adj_Q = Vector{ComplexF64}[]
@@ -85,9 +112,34 @@ function lanczos(mulA!, mulS!, v; niters, lhs=zeros(length(v), 0))
     push!(αs, α)
     push!(lhs_adj_Q, lhs' * v)
 
-    for _ in 2:niters
-        β = sqrt(real(dot(w, Sw)))
-        iszero(β) && break
+    # The maximum number of iterations is length(v)-1, because that is the
+    # dimension of the full vector space. Break out early if niters has been
+    # reached, which will be set according to the spectral bounds Δϵ after
+    # iteration i == min_iters.
+    niters = typemax(Int)
+    for i in 1:length(v)-1
+        if i == min_iters
+            T = SymTridiagonal(αs, βs)
+            Δϵ = eigmax(T) - eigmin(T)
+            niters = max(min_iters, fld(Δϵ, resolution))
+            niters += mod(niters, 2) # Round up to nearest even integer
+            if verbose
+                println("Δϵ=$Δϵ, niters=$niters")
+            end
+        end
+
+        i >= niters && break
+
+        # Let β = w† S w. If β < 0 then S is not a valid positive definite
+        # measure. If β = 0, this would formally indicate that v is a linear
+        # combination of only a subset of the eigenvectors of A S. In this case,
+        # it is valid to break early for purposes of approximating the
+        # matrix-vector product f(A S) v. 
+        β² = real(dot(w, Sw))
+        iszero(β²) && break
+        β² < 0 && error("S is not a positive definite measure")
+
+        β = sqrt(β²)
         @. vp = w / β
         @. Svp = Sw / β
         mulA!(w, Svp)
@@ -131,6 +183,6 @@ function eigbounds(swt, q_reshaped, niters)
     Sv = zero(v)
     mulS!(Sv, v)
     v /= sqrt(v' * Sv)
-    T, _ = lanczos(mulA!, mulS!, v; niters)
+    T, _ = lanczos(mulA!, mulS!, v; min_iters=niters)
     return eigmin(T), eigmax(T)
 end

src/KPM/SpinWaveTheoryKPM.jl

@@ -2,24 +2,16 @@
     SpinWaveTheoryKPM(sys::System; measure, regularization=1e-8, tol)
 
 A variant of [`SpinWaveTheory`](@ref) that uses the kernel polynomial method
-(KPM) to perform [`intensities`](@ref) calculations [1]. Instead of explicitly
-diagonalizing the dynamical matrix, KPM approximates intensities using
-polynomial expansion truncated at order ``M``. The reduces the computational
-cost from ``𝒪(N^3)`` to ``𝒪(N M)``, which is favorable for large system sizes
-``N``.
+(KPM) to calculate [`intensities`](@ref) [1]. This method avoids direct matrix
+diagonalization, which scales cubically in the system size ``N``. Instead, using
+``M`` iterative matrix-vector multiplications, the cost becomes ``𝒪(N M +
+M^2)``. This provides a very significant acceleration when ``N`` is large.
 
-The polynomial order ``M`` will be determined from the line broadening kernel
-and the specified error tolerance `tol`. Specifically, for each wavevector,
-``M`` scales like the spectral bandwidth of excitations, divided by the energy
-resolution of the broadening kernel, times the negative logarithm of `tol`.
-
-The error tolerance `tol` should be tuned empirically for each calculation.
-Reasonable starting points are `1e-1` (more speed) or `1e-2` (more accuracy).
-
-!!! warning "Missing intensity at small quasi-particle energy"  
-    The KPM-calculated intensities are unreliable at small energies ``ω``. In
-    particular, KPM may mask intensities that arise near the Goldstone modes of
-    an ordered state with continuous symmetry.
+The number of iterations is typically not too large: `M ≈ -2 log10(tol) Δϵ /
+fwhm` where `Δϵ` is the estimated spectral bandwidth of excitations, `fwhm` is
+the full width at half maximum of the broadening kernel, and `tol` is a
+dimensionless tolerance parameter. Good choices are `0.05` (more speed) or
+`0.01` (more accuracy).
 
 ## References
 
@@ -185,15 +177,15 @@ end
 function intensities(swt_kpm::SpinWaveTheoryKPM, qpts; energies, kernel::AbstractBroadening, kT=0.0, verbose=false)
     qpts = convert(AbstractQPoints, qpts)
     data = zeros(eltype(swt_kpm.swt.measure), length(energies), length(qpts.qs))
+    # return intensities!(data, swt_kpm, qpts; energies, kernel, kT, verbose)
     return intensities2!(data, swt_kpm, qpts; energies, kernel, kT, verbose)
 end
 
 
 function intensities2!(data, swt_kpm::SpinWaveTheoryKPM, qpts; energies, kernel::AbstractBroadening, kT=0.0, verbose=false)
-    iszero(kT) || error("KPM does not yet support finite kT")
     qpts = convert(AbstractQPoints, qpts)
 
-    (; swt, lanczos_iters) = swt_kpm
+    (; swt, tol, lanczos_iters) = swt_kpm
     (; sys, measure) = swt
     cryst = orig_crystal(sys)
 
@@ -266,12 +258,22 @@ function intensities2!(data, swt_kpm::SpinWaveTheoryKPM, qpts; energies, kernel:
             return w
         end
 
+        resolution = (kernel.fwhm/2) / (-log10(tol))
+
         for ξ in 1:Nobs
             mulA!(v, view(u, :, ξ))
             mulS!(Sv, v)
             c = sqrt(Sv' * v)
             v ./= c
-            tridiag, lhs_adj_Q = lanczos(mulA!, mulS!, v; lhs=u, niters=lanczos_iters)
+            tridiag, lhs_adj_Q = try
+                lanczos(mulA!, mulS!, v; lhs=u, min_iters=lanczos_iters, resolution, verbose)
+            catch e
+                if e.msg == "S is not a positive definite measure"
+                    rethrow(ErrorException("Not an energy-minimum; wavevector q = $q unstable."))
+                else
+                    rethrow()
+                end
+            end
 
             (; values, vectors) = eigen(tridiag)
 
@@ -306,4 +308,3 @@ function intensities2!(data, swt_kpm::SpinWaveTheoryKPM, qpts; energies, kernel:
 
     return Intensities(cryst, qpts, collect(energies), data)
 end
-

src/SpinWaveTheory/DispersionAndIntensities.jl

@@ -97,7 +97,7 @@ function excitations!(T, tmp, swt::SpinWaveTheory, q)
     try
         return bogoliubov!(T, tmp)
     catch _
-        error("Instability at wavevector q = $q")
+        rethrow(ErrorException("Not an energy-minimum; wavevector q = $q unstable."))
     end
 end
 

test/test_kpm.jl

@@ -19,15 +19,15 @@
 
     mulA!(w, v) = mul!(w, A, v)
     mulS!(w, v) = mul!(w, S, v)
-    T, lhs_adj_Q = Sunny.lanczos(mulA!, mulS!, copy(v); lhs, niters=10)
+    T, lhs_adj_Q = Sunny.lanczos(mulA!, mulS!, copy(v); lhs, min_iters=10)
     ev2 = eigvals(T)
 
     @test ev1 ≈ ev2
     @test lhs' * Q ≈ lhs_adj_Q
 
     # Check that extremal eigenvalues match to reasonable accuracy
 
-    T, _ = Sunny.lanczos(mulA!, mulS!, copy(v); niters=20)
+    T, _ = Sunny.lanczos(mulA!, mulS!, copy(v); min_iters=20)
     @test isapprox(eigmin(T), eigmin(A * S); atol=1e-3)
     @test isapprox(eigmax(T), eigmax(A * S); atol=1e-3)
 end