introduce PseudoInverse<:SpectrumSolver

Author: ytdHuang

benchmarks/correlations_and_spectrum.jl

@@ -11,6 +11,8 @@ function benchmark_correlations_and_spectrum!(SUITE)
     ω_l = range(0, 3, length = 1000)
     t_l = range(0, 333 * π, length = 1000)
 
+    PI_solver = PseudoInverse()
+
     function _calculate_fft_spectrum(H, tlist, c_ops, A, B)
         corr = correlation_2op_1t(H, nothing, tlist, c_ops, A, B; progress_bar = Val(false))
         ωlist, spec = spectrum_correlation_fft(tlist, corr)
@@ -23,5 +25,8 @@ function benchmark_correlations_and_spectrum!(SUITE)
     SUITE["Correlations and Spectrum"]["Spectrum"]["Exponential Series"] =
         @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a)
 
+    SUITE["Correlations and Spectrum"]["Spectrum"]["Pseudo Inverse"] =
+        @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a, solver = $PI_solver)
+
     return nothing
 end

docs/src/resources/api.md

@@ -235,6 +235,7 @@ correlation_2op_1t
 spectrum_correlation_fft
 spectrum
 ExponentialSeries
+PseudoInverse
 ```
 
 ## [Metrics](@id doc-API:Metrics)

src/spectrum.jl

@@ -1,5 +1,5 @@
 export spectrum, spectrum_correlation_fft
-export SpectrumSolver, ExponentialSeries#, PseudoInverse
+export SpectrumSolver, ExponentialSeries, PseudoInverse
 
 abstract type SpectrumSolver end
 
@@ -15,6 +15,17 @@ end
 
 ExponentialSeries(; tol = 1e-14, calc_steadystate = false) = ExponentialSeries(tol, calc_steadystate)
 
+@doc raw"""
+    PseudoInverse(; alg::SciMLLinearSolveAlgorithm = KrylovJL_GMRES())
+
+A solver which solves [`spectrum`](@ref) by finding the inverse of Liouvillian [`SuperOperator`](@ref) using the `alg`orithms given in [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/).
+"""
+struct PseudoInverse{MT<:SciMLLinearSolveAlgorithm} <: SpectrumSolver
+    alg::MT
+end
+
+PseudoInverse(; alg::SciMLLinearSolveAlgorithm = KrylovJL_GMRES()) = PseudoInverse(alg)
+
 @doc raw"""
     spectrum(H::QuantumObject,
         ωlist::AbstractVector,
@@ -32,6 +43,7 @@ S(\omega) = \int_{-\infty}^\infty \lim_{t \rightarrow \infty} \left\langle \hat{
 
 See also the following list for `SpectrumSolver` docstrings:
 - [`ExponentialSeries`](@ref)
+- [`PseudoInverse`](@ref)
 """
 function spectrum(
     H::QuantumObject{MT1,HOpType},
@@ -91,18 +103,47 @@ function _spectrum(
     return spec
 end
 
-#= function _spectrum(
+function _spectrum(
     L::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
     ωlist::AbstractVector,
     A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
     B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
     solver::PseudoInverse;
     kwargs...,
 ) where {T1,T2,T3}
-    allequal((L.dims, A.dims, B.dims)) || throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    allequal((L.dims, A.dims, B.dims)) ||
+        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+
+    ωList = convert(Vector{_FType(L)}, ωlist) # Convert it to support GPUs and avoid type instabilities
+    Length = length(ωList)
+    spec = Vector{_FType(L)}(undef, Length)
+
+    # calculate vectorized steadystate, multiply by operator B on the left (spre)
+    ρss = mat2vec(steadystate(L))
+    b = (spre(B) * ρss).data
+
+    # multiply by operator A on the left (spre) and then perform trace operation
+    D = prod(L.dims)
+    _tr = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(_CType(L), D)) # same as vec(system_identity_matrix)
+    _tr_A = transpose(_tr) * spre(A).data
+
+    cache = nothing
+    I_cache = I(D^2)
+    for (idx, ω) in enumerate(ωList)
+        if idx == 1
+            cache = init(LinearProblem(L.data - 1im * ω * I_cache, b), solver.alg, kwargs...)
+            sol = solve!(cache)
+        else
+            cache.A = L.data - 1im * ω * I_cache
+            sol = solve!(cache)
+        end
+
+        # trace over the Hilbert space of system (expectation value)
+        spec[idx] = -2 * real(dot(_tr_A, sol.u))
+    end
 
     return spec
-end =#
+end
 
 @doc raw"""
     spectrum_correlation_fft(tlist, corr; inverse=false)

test/core-test/correlations_and_spectrum.jl

@@ -9,21 +9,25 @@
 
     ω_l2 = range(0, 3, length = 1000)
     spec2 = spectrum(H, ω_l2, c_ops, a', a)
+    spec3 = spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
 
     spec1 = spec1 ./ maximum(spec1)
     spec2 = spec2 ./ maximum(spec2)
+    spec3 = spec3 ./ maximum(spec3)
 
     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)
     idxs1 = test_func1 .> 0.05
     idxs2 = test_func2 .> 0.05
     @test sum(abs2.(spec1[idxs1] .- test_func1[idxs1])) / sum(abs2.(test_func1[idxs1])) < 0.01
     @test sum(abs2.(spec2[idxs2] .- test_func2[idxs2])) / sum(abs2.(test_func2[idxs2])) < 0.01
+    @test all(spec2 .≈ spec3)
 
     @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, corr)
         @inferred spectrum(H, ω_l2, c_ops, a', a)
+        @inferred spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
     end
 
     @test_throws ErrorException FFTCorrelation()