Replace KPM with Lanczos

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,13 +1,40 @@
-# 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[]
+    vs = typeof(v)[]
+
+    c = sqrt(real(dot(v, S, v)))
+    @assert c ≈ 1 "Initial vector not normalized"
+    v /= c
 
-    v /= sqrt(real(dot(v, S, v)))
     w = A * S * v
     α = real(dot(w, S, v))
     w = w - α * v
     push!(αs, α)
+    push!(vs, v)
 
     for _ in 2:niters
         β = sqrt(real(dot(w, S, w)))
@@ -19,35 +46,54 @@ function lanczos_ref(A, S, v; niters)
         v = vp
         push!(αs, α)
         push!(βs, β)
+        push!(vs, v)
     end
 
-    return SymTridiagonal(αs, βs)
+    Q = reduce(hcat, vs)
+    T = SymTridiagonal(αs, βs)
+    return (; Q, T)
 end
 
 
-
-# Use Lanczos iterations to find a truncated tridiagonal representation of A S,
-# up to similarity transformation. Here, A is any Hermitian matrix, while S is
-# both Hermitian and positive definite. Traditional Lanczos uses the identity
-# matrix, S = I. The extension to non-identity matrices S is as follows: Each
-# matrix-vector product A v becomes A S v, and each vector inner product w† v
-# becomes w† S v. The implementation below follows the notation of Wikipedia,
-# and is the most stable of the four variants considered by Paige [1]. Each
-# Lanczos iteration requires only one matrix-vector multiplication for A and S,
-# respectively.
-
-# Similar generalizations of Lanczos have been considered in [2] and [3].
+# 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.
+#
+# In this implementation, each Lanczos iteration requires only one matrix-vector
+# multiplication involving A and S, respectively. 
+#
+# Details:
+#
+# * 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. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
-# https://doi.org/10.1093%2Fimamat%2F10.3.373.
-# 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(mulA!, mulS!, v; niters)
+#    https://doi.org/10.1093%2Fimamat%2F10.3.373.
+
+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}[]
 
+    v = copy(v)
     vp  = zero(v)
     Sv  = zero(v)
     Svp = zero(v)
@@ -56,6 +102,7 @@ function lanczos(mulA!, mulS!, v; niters)
 
     mulS!(Sv, v)
     c = sqrt(real(dot(v, Sv)))
+    @assert c ≈ 1 "Initial vector not normalized"
     @. v /= c
     @. Sv /= c
 
@@ -64,10 +111,36 @@ function lanczos(mulA!, mulS!, v; niters)
     @. w = w - α * v
     mulS!(Sw, w)
     push!(αs, α)
-
-    for _ in 2:niters
-        β = sqrt(real(dot(w, Sw)))
-        iszero(β) && break
+    push!(lhs_adj_Q, lhs' * v)
+
+    # 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)
@@ -77,9 +150,10 @@ function lanczos(mulA!, mulS!, v; niters)
         @. v = vp
         push!(αs, α)
         push!(βs, β)
+        push!(lhs_adj_Q, lhs' * v)
     end
 
-    return SymTridiagonal(αs, βs)
+    return SymTridiagonal(αs, βs), reduce(hcat, lhs_adj_Q)
 end
 
 
@@ -107,6 +181,9 @@ function eigbounds(swt, q_reshaped, niters)
     end
 
     v = randn(ComplexF64, 2L)
-    tridiag = lanczos(mulA!, mulS!, v; niters)
-    return eigmin(tridiag), eigmax(tridiag)
+    Sv = zero(v)
+    mulS!(Sv, v)
+    v /= sqrt(v' * Sv)
+    T, _ = lanczos(mulA!, mulS!, v; min_iters=niters)
+    return eigmin(T), eigmax(T)
 end