Use FastExpm.jl for matrix exponential of SparseSuperOperator (#112)

---------

Co-authored-by: Stefan Krastanov <[email protected]>

Project.toml

@@ -1,10 +1,11 @@
 name = "QuantumOpticsBase"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "0.4.11"
+version = "0.4.12"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -17,6 +18,7 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 [compat]
 Adapt = "1, 2, 3.3"
 FFTW = "1.2"
+FastExpm = "1.1.0"
 FillArrays = "0.13, 1"
 LRUCache = "1"
 QuantumInterface = "0.3.0"

src/operators_dense.jl

@@ -226,6 +226,15 @@ function expect(op::DataOperator{B1,B2}, state::DataOperator{B2,B2}) where {B1,B
     result
 end
 
+"""
+    exp(op::DenseOpType)
+
+Operator exponential used, for example, to calculate displacement operators.
+Uses LinearAlgebra's `Base.exp`.
+
+If you only need the result of the exponential acting on a vector,
+consider using much faster implicit methods that do not calculate the entire exponential.
+"""
 function exp(op::T) where {B,T<:DenseOpType{B,B}}
     return DenseOperator(op.basis_l, op.basis_r, exp(op.data))
 end

src/operators_sparse.jl

@@ -1,5 +1,6 @@
 import Base: ==, *, /, +, -, Broadcast
 import SparseArrays: sparse
+import FastExpm: fastExpm
 
 const SparseOpPureType{BL,BR} = Operator{BL,BR,<:SparseMatrixCSC}
 const SparseOpAdjType{BL,BR} = Operator{BL,BR,<:Adjoint{<:Number,<:SparseMatrixCSC}}
@@ -45,6 +46,21 @@ function expect(op::SparseOpPureType{B1,B2}, state::Operator{B2,B2}) where {B1,B
     result
 end
 
+"""
+    exp(op::SparseOpType; opts...)
+
+Operator exponential used, for example, to calculate displacement operators.
+Uses [`FastExpm.jl.jl`](https://github.com/fmentink/FastExpm.jl) which will return a sparse
+or dense operator depending on which is more efficient.
+All optional arguments are passed to `fastExpm` and can be used to specify tolerances.
+
+If you only need the result of the exponential acting on a vector,
+consider using much faster implicit methods that do not calculate the entire exponential.
+"""
+function exp(op::T; opts...) where {B,T<:SparseOpType{B,B}}
+    return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+end
+
 function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis}
     @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
     @assert isperm(perm)
@@ -99,4 +115,4 @@ mul!(result::Bra{B2},b::Bra{B1},M::SparseOpPureType{B1,B2},alpha,beta) where {B1
 /(op::EyeOpType,a::T) where {T<:Number} = sparse(op)/a
 tensor(a::EyeOpType, b::SparseOpType) = tensor(sparse(a),b)
 tensor(a::SparseOpType, b::EyeOpType) = tensor(a,sparse(b))
-tensor(a::EyeOpType, b::EyeOpType) = tensor(sparse(a),sparse(b))
+tensor(a::EyeOpType, b::EyeOpType) = tensor(sparse(a),sparse(b))

src/superoperators.jl

@@ -1,4 +1,5 @@
 import QuantumInterface: AbstractSuperOperator
+import FastExpm: fastExpm
 
 """
     SuperOperator <: AbstractSuperOperator
@@ -221,10 +222,27 @@ end
 """
     exp(op::DenseSuperOperator)
 
-Operator exponential which can for example used to calculate time evolutions.
+Superoperator exponential which can, for example, be used to calculate time evolutions.
+Uses LinearAlgebra's `Base.exp`.
+
+If you only need the result of the exponential acting on an operator,
+consider using much faster implicit methods that do not calculate the entire exponential.
 """
 Base.exp(op::DenseSuperOpType) = DenseSuperOperator(op.basis_l, op.basis_r, exp(op.data))
 
+"""
+    exp(op::SparseSuperOperator; opts...)
+
+Superoperator exponential which can, for example, be used to calculate time evolutions.
+Uses [`FastExpm.jl.jl`](https://github.com/fmentink/FastExpm.jl) which will return a sparse
+or dense operator depending on which is more efficient.
+All optional arguments are passed to `fastExpm` and can be used to specify tolerances.
+
+If you only need the result of the exponential acting on an operator,
+consider using much faster implicit methods that do not calculate the entire exponential.
+"""
+Base.exp(op::SparseSuperOpType; opts...) = SuperOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+
 # Array-like functions
 Base.size(A::SuperOperator) = size(A.data)
 @inline Base.axes(A::SuperOperator) = axes(A.data)

test/test_operators.jl

@@ -114,7 +114,10 @@ OP_cnot = embed(b8, [1,3], op_cnot)
 @test reduced(OP_cnot, [1,3])/2. == op_cnot
 @test_throws AssertionError embed(b2⊗b2, [1,1], op_cnot)
 
-@test_throws ArgumentError exp(sparse(op1))
+b = FockBasis(40)
+alpha = 1+5im
+H = alpha * create(b) - conj(alpha) * destroy(b)
+@test exp(sparse(H); threshold=1e-10) ≈ displace(b, alpha)
 
 @test one(b1).data == Diagonal(ones(b1.shape[1]))
 @test one(op1).data == Diagonal(ones(b1.shape[1]))

test/test_superoperators.jl

@@ -188,6 +188,7 @@ end
 
 # tout, ρt = timeevolution.master([0.,1.], ρ₀, H, J; reltol=1e-7)
 # @test tracedistance(exp(dense(L))*ρ₀, ρt[end]) < 1e-6
+# @test tracedistance(exp(sparse(L))*ρ₀, ρt[end]) < 1e-6
 
 @test dense(spre(op1)) == spre(op1)
 
@@ -217,4 +218,15 @@ Ldense .+= L
 @test_throws ErrorException cos.(Ldense)
 @test_throws ErrorException cos.(L)
 
+b = FockBasis(20)
+L = liouvillian(identityoperator(b), [destroy(b)])
+@test exp(sparse(L)).data ≈ exp(dense(L)).data
+N = exp(log(2) * sparse(L)) # 50% loss channel
+ρ = N * dm(coherentstate(b, 1))
+@test (1 - real(tr(ρ^2))) < 1e-10 # coherent state remains pure under loss
+@test tracedistance(ρ, dm(coherentstate(b, 1/sqrt(2)))) < 1e-10
+ρ = N * dm(fockstate(b, 1))
+@test (0.5 - real(tr(ρ^2))) < 1e-10 # one photon state becomes maximally mixed
+@test tracedistance(ρ, normalize(dm(fockstate(b, 0)) + dm(fockstate(b, 1)))) < 1e-10
+
 end # testset