Add support to `LinearMaps.jl`

[GiggleLiu]

LinearMaps.jl can map a matrix vector multiplication to a linear operator.

Is it possible to support expmv fallback for a general linear operator? I need it very much in my current project.

see this pr acroy/Expokit.jl#30

[matteoacrossi]

@GiggleLiu I added support for LinearMaps.jl in this pull request of the currently active ExpmV.jl repository.

[GiggleLiu]

Cool! This interface looks good and the code is qualified. I will definitely have a try!

I'd like to post some benchmark under this issue latter, by comparing its performance in my application with Expokit's realization.

Thanks for your efforts!

[GiggleLiu]

First, your implementation is correct!

But somehow less efficient than Expokit

vector dimension 4096

Expokit

BenchmarkTools.Trial: 
  memory estimate:  2.24 MiB
  allocs estimate:  134
  --------------
  minimum time:     190.535 ms (0.00% GC)
  median time:      194.090 ms (0.00% GC)
  mean time:        199.734 ms (1.57% GC)
  maximum time:     254.911 ms (24.27% GC)
  --------------
  samples:          26
  evals/sample:     1

ExpmV

BenchmarkTools.Trial: 
  memory estimate:  31.86 MiB
  allocs estimate:  258890
  --------------
  minimum time:     12.443 s (0.16% GC)
  median time:      12.443 s (0.16% GC)
  mean time:        12.443 s (0.16% GC)
  maximum time:     12.443 s (0.16% GC)
  --------------
  samples:          1
  evals/sample:     1

Vector dimension 1024, but more time consuming operator

Expokit

BenchmarkTools.Trial: 
  memory estimate:  176.94 MiB
  allocs estimate:  3677
  --------------
  minimum time:     2.393 s (4.12% GC)
  median time:      2.464 s (7.00% GC)
  mean time:        2.478 s (7.61% GC)
  maximum time:     2.577 s (11.43% GC)
  --------------
  samples:          3
  evals/sample:     1

ExpmV

BenchmarkTools.Trial: 
  memory estimate:  209.85 MiB
  allocs estimate:  93982
  --------------
  minimum time:     3.602 s (4.90% GC)
  median time:      3.651 s (6.74% GC)
  mean time:        3.651 s (6.74% GC)
  maximum time:     3.700 s (8.53% GC)
  --------------
  samples:          2
  evals/sample:     1

The allocation is very high, would you please check the type stability?

The following profile information may be helpful to you.

17171 ./task.jl:85; (::getfield(Atom, Symbol("##108#113")){Dict{String,Any}})()
 17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:89; macro expansion
  17171 /home/leo/.julia/packages/Atom/Pab0Z/src/repl.jl:76; hideprompt(::getfield(Atom, Symbol("##109#114")){String,Int64,String})
   17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:90; #109
    17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:46; withpath
     17171 /home/leo/.julia/packages/CodeTools/8CjY/src/utils.jl:30; withpath(::getfield(Atom, Symbol("##110#115")){String,Int64,String}, ::String)
      17171 /home/leo/.julia/packages/Atom/Pab0Z/src/eval.jl:91; (::getfield(Atom, Symbol("##110#115")){String,Int64,String})()
       17171 /home/leo/.julia/packages/CodeTools/8CjY/src/eval.jl:30; include_string(::Module, ::String, ::String, ::Int64)
        17171 ./loading.jl:1002; include_string(::Module, ::String, ::String)
         17171 ...x64/build/usr/share/julia/stdlib/v1.0/Profile/src/Profile.jl:25; top-level scope
          2535  /home/leo/.julia/dev/Yao/src/Blocks/TimeEvolution.jl:23; apply!(::DefaultRegister{1,Complex{Float64}}, ::TimeEvolution{10,Float64,ChainBlock{10,Comp...
           2535 ./reducedim.jl:305; mat
            2463 ./reduce.jl:314; _mapreduce(::getfield(Yao.Blocks, Symbol("##41#42")), ::typeof(Base.mul_prod), ::IndexLine...
             2334 ./reduce.jl:31; mul_prod(::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMatrixCSC{Complex{Float64},Int...
              2334 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:141; *
               2334 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:151; spmatmul
                823 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:178; #spmatmul#60(::Symbol, ::Function, ::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMa...
                 823 ./complex.jl:268; *
                  562 ./float.jl:399; *
                943 ...ild/usr/share/julia/stdlib/v1.0/SparseArrays/src/linalg.jl:186; #spmatmul#60(::Symbol, ::Function, ::SparseMatrixCSC{Complex{Float64},Int64}, ::SparseMa...
                 679 ./complex.jl:266; +
                  679 ./float.jl:395; +
          14433 /home/leo/jcode/lab/TE_benchmark.jl:21; apply2!(::DefaultRegister{1,Complex{Float64}}, ::TimeEvolution{12,Float64,PutBlock{12,1,XGa...
           14433 /home/leo/.julia/packages/ExpmV/GFsI0/src/expmv_fun.jl:26; expmv
            14123 /home/leo/.julia/packages/ExpmV/GFsI0/src/expmv_fun.jl:30; #expmv#1(::Nothing, ::String, ::Bool, ::Bool, ::Function, ::Complex{Float64}, ::LuxurySpar...
             14123 ...leo/.julia/packages/ExpmV/GFsI0/src/select_taylor_degree.jl:63; macro expansion
              13297 ...eo/.julia/packages/ExpmV/GFsI0/src/select_taylor_degree.jl:98; #select_taylor_degree#3(::Int64, ::Int64, ::String, ::Bool, ::Bool, ::Function, ::Luxury...
               13296 /home/leo/.julia/packages/ExpmV/GFsI0/src/normAm.jl:30; normAm
                8389 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:112; norm1est(::Int64, ::LuxurySparse.PermMatrix{Complex{Float64},Int64,Array{Complex{Float64...
                 1525 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:33; A_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int64...
                  1525 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
                   1525 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMat...
                    680 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:638; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMa...
                     679 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:480; copy_transpose!
                      620 .../share/julia/stdlib/v1.0/LinearAlgebra/src/transpose.jl:197; copy_transpose!(::Array{Complex{Float64},2}, ::UnitRange{Int64}, ::UnitRange{Int64}, ...
                    732 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:646; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMa...
                 6862 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:35; A_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int64...
                  6862 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
                   6862 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMat...
                    3005 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:638; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
                     2997 ...usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:480; copy_transpose!
                      2762 .../share/julia/stdlib/v1.0/LinearAlgebra/src/transpose.jl:197; copy_transpose!(::Array{Complex{Float64},2}, ::UnitRange{Int64}, ::UnitRange{Int64},...
                       780  ./array.jl:771; setindex!
                       880  ./promotion.jl:425; getindex
                       1102 /home/leo/.julia/dev/LuxurySparse/src/PermMatrix.jl:46; getindex
                        882 ./array.jl:731; getindex
                    3364 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:646; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
                     1574 ./complex.jl:268; *
                      1322 ./float.jl:399; *
                     1530 ./complex.jl:266; +
                      1530 ./float.jl:395; +
                4878 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:182; norm1est(::Int64, ::LuxurySparse.PermMatrix{Complex{Float64},Int64,Array{Complex{Float64...
                 955  /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:42; At_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int6...
                  954 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
                   954 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMatr...
                 3923 /home/leo/.julia/packages/ExpmV/GFsI0/src/norm1est.jl:44; At_pow_n_B!(::Array{Complex{Float64},2}, ::LuxurySparse.PermMatrix{Complex{Float64},Int6...
                  3922 ...d/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:172; mul!
                   3922 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:578; generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermMat...
                    1705 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:638; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
                     1698 ...usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:480; copy_transpose!
                      1556 .../share/julia/stdlib/v1.0/LinearAlgebra/src/transpose.jl:197; copy_transpose!(::Array{Complex{Float64},2}, ::UnitRange{Int64}, ::UnitRange{Int64},...
                       618 /home/leo/.julia/dev/LuxurySparse/src/PermMatrix.jl:46; getindex
                    1934 .../usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:646; _generic_matmatmul!(::Array{Complex{Float64},2}, ::Char, ::Char, ::LuxurySparse.PermM...
                     962 ./complex.jl:268; *
                      792 ./float.jl:399; *
                     786 ./complex.jl:266; +
                      786 ./float.jl:395; +

[matteoacrossi]

Thanks a lot for your thorough testing. There is definitely some issue with allocations and I will take a deeper look at it.

[matteoacrossi]

Anyways, from the profile info that you posted it is evident how almost all the time is spent in the mul! function, which is internal. I don't think that much can be done by reducing allocations.