Lanczos-based spectrum method

benchmarks/correlations_and_spectrum.jl

@@ -19,6 +19,8 @@ function benchmark_correlations_and_spectrum!(SUITE)
 
     PI_solver = PseudoInverse()
 
+    L_solver = Lanczos()
+
     SUITE["Correlations and Spectrum"]["FFT Correlation"] =
         @benchmarkable _calculate_fft_spectrum($H, $t_l, $c_ops, $(a'), $a)
 
@@ -28,5 +30,8 @@ function benchmark_correlations_and_spectrum!(SUITE)
     SUITE["Correlations and Spectrum"]["Spectrum"]["Pseudo Inverse"] =
         @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a, solver = $PI_solver)
 
+    SUITE["Correlations and Spectrum"]["Spectrum"]["Lanczos"] =
+        @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a, solver = $L_solver)
+
     return nothing
 end

docs/src/resources/api.md

@@ -248,6 +248,7 @@ spectrum_correlation_fft
 spectrum
 ExponentialSeries
 PseudoInverse
+Lanczos
 ```
 
 ## [Entropy and Metrics](@id doc-API:Entropy-and-Metrics)

src/QuantumToolbox.jl

@@ -115,14 +115,14 @@ include("time_evolution/time_evolution_dynamical.jl")
 
 # Others
 include("correlations.jl")
-include("spectrum.jl")
 include("wigner.jl")
 include("spin_lattice.jl")
 include("arnoldi.jl")
 include("entropy.jl")
 include("metrics.jl")
 include("negativity.jl")
 include("steadystate.jl")
+include("spectrum.jl")
 include("visualization.jl")
 
 # deprecated functions

src/spectrum.jl

@@ -1,5 +1,5 @@
 export spectrum, spectrum_correlation_fft
-export SpectrumSolver, ExponentialSeries, PseudoInverse
+export SpectrumSolver, ExponentialSeries, PseudoInverse, Lanczos
 
 abstract type SpectrumSolver end
 
@@ -26,6 +26,20 @@
 
 PseudoInverse(; alg::SciMLLinearSolveAlgorithm = KrylovJL_GMRES()) = PseudoInverse(alg)
 
+@doc raw"""
+    Lanczos(; tol = 1e-8, maxiter = 5000, verbose = 0)
+
+A solver which solves [`spectrum`](@ref) by using a non-symmetric Lanczos variant of the algorithm in [Koch2011](https://www.cond-mat.de/events/correl11/manuscripts/koch.pdf).
+The nonsymmetric Lanczos algorithm is adapted from Algorithm 6.6 in [Saad2011](https://www-users.cse.umn.edu/~saad/eig_book_2ndEd.pdf).
+The running estimate is updated via a [Wallis-Euler recursion](https://en.wikipedia.org/wiki/Continued_fraction).
+"""
+Base.@kwdef mutable struct Lanczos{T<:Real,IT<:Int} <: SpectrumSolver
+    tol::T = 1e-8
+    maxiter::IT = 5000
+    verbose::IT = 0
+    steadystate_solver::Union{Nothing,SteadyStateSolver} = nothing
+end
+
 @doc raw"""
     spectrum(H::QuantumObject,
         ωlist::AbstractVector,
@@ -44,6 +58,7 @@
 See also the following list for `SpectrumSolver` docstrings:
 - [`ExponentialSeries`](@ref)
 - [`PseudoInverse`](@ref)
+- [`Lanczos`](@ref)
 """
 function spectrum(
     H::QuantumObject{HOpType},
@@ -144,6 +159,142 @@
     return spec
 end
 
+function _spectrum(
+    L::QuantumObject{SuperOperator},
+    ωlist::AbstractVector,
+    A::QuantumObject{Operator},
+    B::QuantumObject{Operator},
+    solver::Lanczos;
+)
+    check_dimensions(L, A, B)
+
+    # Define type shortcuts
+    fT = _FType(L)
+    cT = _CType(L)
+
+    # Calculate |v₁> = B|ρss>
+    ρss =
+        solver.steadystate_solver === nothing ? mat2vec(steadystate(L)) :
+        mat2vec(steadystate(L; solver = solver.steadystate_solver))
+    vₖ = (spre(B) * ρss).data
+
+    # Define (possibly GPU) vector type
+    vT = typeof(vₖ)
+
+    # Calculate <w₁| = <I|A
+    D = prod(L.dimensions)
+    Ivec = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(cT, D)) # same as vec(system_identity_matrix)
+    wₖ = transpose(vT(Ivec)) * spre(A).data
+
+    # Store the norm of the Green's function before renormalizing |v₁> and <w₁|
+    gfNorm = abs(wₖ * vₖ)
+    vₖ ./= sqrt(gfNorm)
+    wₖ ./= sqrt(gfNorm)
+
+    # Handle input frequency range
+    ωList = vT(convert(Vector{fT}, ωlist))  # Make sure they're real frequencies and potentially on GPU
+    Length = length(ωList)
+
+    # Current and previous estimates of the spectrum
+    lanczosFactor = vT(zeros(cT, Length))
+    lanczosFactor₋₁ = vT(zeros(cT, Length))
+
+    # Tridiagonal matrix elements
+    αₖ = cT(0)
+    βₖ = cT(-1)
+    δₖ = cT(+1)
+
+    # Current and up to second-to-last A and B Euler sequences
+    A₋₂ = vT(ones(cT, Length))
+    A₋₁ = vT(zeros(cT, Length))
+    Aₖ = vT(zeros(cT, Length))
+    B₋₂ = vT(zeros(cT, Length))
+    B₋₁ = vT(ones(cT, Length))
+    Bₖ = vT(zeros(cT, Length))
+
+    # Maximum norm and residue
+    maxNorm = vT(zeros(cT, length(ωList)))
+    maxResidue = fT(0.0)
+
+    # Previous and next left/right Krylov vectors
+    v₋₁ = vT(zeros(cT, D^2))
+    v₊₁ = vT(zeros(cT, D^2))
+    w₋₁ = vT(zeros(cT, D^2))'
+    w₊₁ = vT(zeros(cT, D^2))'
+
+    # Frequency of renormalization
+    renormFrequency::typeof(solver.maxiter) = 1
+
+    # Loop over the Krylov subspace(s)
+    for k in 1:solver.maxiter
+        # k-th diagonal element
+        w₊₁ = wₖ * L.data
+        αₖ = w₊₁ * vₖ
+
+        # Update A(k), B(k) and continuous fraction; normalization avoids overflow
+        Aₖ .= (-1im .* ωList .+ αₖ) .* A₋₁ .- (βₖ * δₖ) .* A₋₂
+        Bₖ .= (-1im .* ωList .+ αₖ) .* B₋₁ .- (βₖ * δₖ) .* B₋₂
+        lanczosFactor₋₁ .= lanczosFactor
+        lanczosFactor .= Aₖ ./ Bₖ
+
+        # Renormalize Euler sequences to avoid overflow
+        if k % renormFrequency == 0
+            maxNorm .= max.(abs.(Aₖ), abs.(Bₖ))
+            Aₖ ./= maxNorm
+            Bₖ ./= maxNorm
+            A₋₁ ./= maxNorm
+            B₋₁ ./= maxNorm
+        end
+
+        # Check for convergence
+        maxResidue =
+            maximum(abs.(lanczosFactor .- lanczosFactor₋₁)) /
+            max(maximum(abs.(lanczosFactor)), maximum(abs.(lanczosFactor₋₁)))
+        if maxResidue <= solver.tol
+            if solver.verbose > 1
+                println("spectrum(): solver::Lanczos converged after $(k) iterations")
+            end
+            break
+        end
+
+        # (k+1)-th left/right vectors, orthogonal to previous ones
+        # Consider using explicit BLAS calls
+        v₊₁ = L.data * vₖ
+        v₊₁ .= v₊₁ .- αₖ .* vₖ .- βₖ .* v₋₁
+        w₊₁ .= w₊₁ .- αₖ .* wₖ .- δₖ .* w₋₁
+        v₋₁ .= vₖ
+        w₋₁ .= wₖ
+        vₖ .= v₊₁
+        wₖ .= w₊₁
+
+        # k-th off-diagonal elements
+        buf = wₖ * vₖ
+        δₖ = sqrt(abs(buf))
+        βₖ = buf / δₖ
+
+        # Normalize (k+1)-th left/right vectors
+        vₖ ./= δₖ
+        wₖ ./= βₖ
+
+        # Update everything for the next cycle
+        A₋₂ .= A₋₁
+        A₋₁ .= Aₖ
+        B₋₂ .= B₋₁
+        B₋₁ .= Bₖ
+    end
+
+    if solver.verbose > 0 && maxResidue > solver.tol
+        println("spectrum(): maxiter = $(solver.maxiter) reached before convergence!")
+        println("spectrum(): Max residue = $maxResidue")
+        println("spectrum(): Consider increasing maxiter and/or tol")
+    end
+
+    # Restore the norm
+    lanczosFactor .= gfNorm .* lanczosFactor
+
+    return -2 .* real(lanczosFactor)
+end
+
 @doc raw"""
     spectrum_correlation_fft(tlist, corr; inverse=false)
 

test/core-test/correlations_and_spectrum.jl

@@ -13,10 +13,12 @@
     ω_l2 = range(0, 3, length = 1000)
     spec2 = spectrum(H, ω_l2, c_ops, a', a)
     spec3 = spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
+    spec4 = spectrum(H, ω_l2, c_ops, a', a; solver = Lanczos())
 
     spec1 = spec1 ./ maximum(spec1)
     spec2 = spec2 ./ maximum(spec2)
     spec3 = spec3 ./ maximum(spec3)
+    spec4 = spec4 ./ maximum(spec4)
 
     test_func1 = maximum(real.(spec1)) * (0.1 / 2)^2 ./ ((ω_l1 .- 1) .^ 2 .+ (0.1 / 2)^2)
     test_func2 = maximum(real.(spec2)) * (0.1 / 2)^2 ./ ((ω_l2 .- 1) .^ 2 .+ (0.1 / 2)^2)
@@ -26,12 +28,14 @@
     @test sum(abs2.(spec2[idxs2] .- test_func2[idxs2])) / sum(abs2.(test_func2[idxs2])) < 0.01
     @test all(corr1 .≈ corr2)
     @test all(spec2 .≈ spec3)
+    @test all(spec2 .≈ spec4)
 
     @testset "Type Inference spectrum" begin
         @inferred correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
         @inferred spectrum_correlation_fft(t_l, corr1)
         @inferred spectrum(H, ω_l2, c_ops, a', a)
         @inferred spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
+        @inferred spectrum(H, ω_l2, c_ops, a', a; solver = Lanczos())
     end
 
     # tlist and τlist checks