`exp` does not work on many sparse `LinearAlgebra` matrices

[Krastanov]

julia> exp(Bidiagonal([1.,2.,3.],[4.,5.],:U))
ERROR: MethodError: no method matching exp(::Bidiagonal{Float64,Array{Float64,1}})

It fails similarly for some other sparse matrices.

It works if I explicitly call collect on the sparse matrix, but I was expecting Julia's promotion mechanism / multiple dispatch to make that unnecessary?

Is this some form of protection against accidentally making gigantic dense matrices? Even if that is the case, I would have expected there to be a matrix exponentiation function that "just works", even if in some cases that involves an expensive operation.

Observed both in julia 1.2 and 1.3.

[mbauman]

The exp method on dense matrices actually requires dense matrices as it depends upon BLAS and LAPACK to do some of its work. We simply don't have a generic exp that can work on all matrix types.

Julia won't do any promotion/conversion behind your back when calling functions. Either we have a method for it or we don't. In this case we don't.

[Krastanov]

Julia won't do any promotion/conversion behind your back when calling functions. Either we have a method for it or we don't. In this case we don't.

I guess the fact that exp(Int(1)) returns a float, made me expect exp(sparsematrix) to return a dense matrix if necessary. In other words, does it not make sense to just default to the expensive slow dense algorithm instead of expecting the user to call collect?

Simply a method for some common abstract supertype of (Bi/Tri)diagonal that can be overloaded in the future if someone devises a smart algorithm. Something like exp(m::AbstractSparse) = exp(collect(m))?

This seems especially important given that without it, code that tries to use exp without runtime checks of types can not be written. For instance, I have a function f that uses exp and I want my function f to work on both dense and sparse matrices. Now I have to put an if statement that checks the type, or write two versions of f, or write my own internal exp, all of which are easy, but feel silly to do when Julia already can do multiple dispatch for exp.

[mbauman]

Yes, you're absolutely right that we could achieve this like:

exp(A::AbstractMatrix{T}) where {T<:BlasFloat} = exp!(Matrix(A))

That actually wouldn't be such a horrible implementation, but of course certain types could do much much better. My point is that this isn't automatic — we need to decide to do this and someone needs to implement it.

[Krastanov]

What would be the process to petition for this to be included? Should I just try to make a pull request (obviously, the first pull request I do will probably get it completely wrong)?

[Krastanov]

Related question: I guess the implementation that you sketched actually needs to go in something like Base, not LinearAlgebra? The source code tree looks a bit intimidating, so I would appreciate a pointer to where roughly to place it.

[antoine-levitt]

I think it makes sense to leave these operations (exp, eigen, etc.) undefined, since they create dense matrices anyway, and can't be implemented more efficiently than by converting the input to dense. It's good to make the user be explicit here, so that it's clear that there aren't any clever tricks going on. AFAICT the only defined operation that creates fully dense matrices from sparse ones is broadcasting, ie a .+ 1, and even the utility of that is debatable...

[Krastanov]

But then the following very common situation is made unnecessarily complicated:

function some_heavy_computation(dense_or_sparse_matrix)
    intermediate_result = computation_that_can_be_very_fast_if_sparse(dense_or_sparse_matrix)
    return exp(intermediate_result)
end

In this simple case a single convert would be a clean solution, but the computation sequence can be more complicated and cumbersome to write down, especially when one defines a ton of small one-liner functions to make the computation more semantically readable.

[antoine-levitt]

I'm curious, what's your use case? I would imagine the computational cost for exp to completely dwarf any intermediate computation anyway, so I would just do everything as dense.

[Krastanov]

@antoine-levitt, my case involves some quantum mechanics simulation where I have to construct a matrix (the Hamiltonian) by first dealing with a lot of sparse matrices (products of Diagonal, Bidiagonal, etc). When that expression is finally constructed (and it can be a cumbersome one), I call exp.

The computational/memory cost, as you guessed, is not that significant compared to what one needs to expend anyway for the exponentiation. It is measurable, but not worth worrying about at this point. However, the legibility of the code definitely suffers.

[ryanengle]

FWIW, @Krastanov, @antoine-levitt, my research group has encountered a similar problem related to quantum mechanics involving sparse matrices.

[fab4ap]

Looks like I have the same problem trying to calculate Planck's black body radiation: ERROR: LoadError: MethodError: no method matching exp(::Array{Transpose{Float64,Array{Float64,1}},1})

Had to resort to Octave/Matlab.

[ViralBShah]

Should be a simple PR for now to collect it and call exp at least for the dense case. We should definitely add that.

[ViralBShah]

If exp(A') is the same as exp(A)' that's an easier rule to add that may work in more cases.

[fredrikekre]

ERROR: LoadError: MethodError: no method matching exp(::Array{Transpose{Float64,Array{Float64,1}},1})

That is a weird type to take exp of -- what do you want to achieve by taking exp of a Vector with Transposes of Vectors? You probably meant to do it elementwise somehow? Note that Julia's exp is a matrix exponential defined only for square matrices whereas matlab (and presumably octave) applies it element wise.

[fab4ap]

Probably. Let me tell you the context, in case it helps to make Julia more popular. I'm a Research Ecologist and discovered that Wein's/Planck's black body radiation responses are very similar to tree population growth response. So, I wanted to implement that in Julia to take advantage of some graphing abilities for a publication. I'm not really savvy about the myriad Types in Julia - and perhaps not willing to invest the time to figure it all out. It worked in matlab and R (but I wasn't happy with the graphs) - so I tried to make it work in Julia (recently I've tried to switch to Julia from Python for data manipulation and analysis). When an end-user like me encounters errors like that, it is rather mystifying - and reminiscent of class-based compiled languages that I've been trying to avoid. I perhaps should have known better, but I plead ignorance about the intricate Types. String, integer and floating point incorporates > 90% of what I do! As I start to encounter more of such errors in Julia, I wish it was not so complex.

On Thu, Mar 26, 2020 at 3:40 PM Fredrik Ekre notifications@github.com wrote:

ERROR: LoadError: MethodError: no method matching exp(::Array{Transpose{Float64,Array{Float64,1}},1})

That is a weird type to take exp of -- what do you want to achieve by taking exp of a Vector with Transposes of Vectors? You probably meant to do it elementwise somehow? Note that Julia's exp is a matrix exponential https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.exp-Tuple%7BUnion%7BDenseArray%7B%23s4,2%7D,%20Base.ReinterpretArray%7B%23s4,2,S,A%7D%20where%20S%20where%20A%3C:Union%7BSubArray%7BT,N,A,I,true%7D%20where%20I%3C:Union%7BTuple%7BVararg%7BReal,N%7D%20where%20N%7D,%20Tuple%7BAbstractUnitRange,Vararg%7BAny,N%7D%20where%20N%7D%7D%20where%20A%3C:DenseArray%20where%20N%20where%20T,%20DenseArray%7D,%20Base.ReshapedArray%7B%23s4,2,A,MI%7D%20where%20MI%3C:Tuple%7BVararg%7BBase.MultiplicativeInverses.SignedMultiplicativeInverse%7BInt64%7D,N%7D%20where%20N%7D%20where%20A%3C:Union%7BBase.ReinterpretArray%7BT,N,S,A%7D%20where%20S%20where%20A%3C:Union%7BSubArray%7BT,N,A,I,true%7D%20where%20I%3C:Union%7BTuple%7BVararg%7BReal,N%7D%20where%20N%7D,%20Tuple%7BAbstractUnitRange,Vararg%7BAny,N%7D%20where%20N%7D%7D%20where%20A%3C:DenseArray%20where%20N%20where%20T,%20DenseArray%7D%20where%20N%20where%20T,%20SubArray%7BT,N,A,I,true%7D%20where%20I%3C:Union%7BTuple%7BVararg%7BReal,N%7D%20where%20N%7D,%20Tuple%7BAbstractUnitRange,Vararg%7BAny,N%7D%20where%20N%7D%7D%20where%20A%3C:DenseArray%20where%20N%20where%20T,%20DenseArray%7D,%20SubArray%7B%23s4,2,A,I,L%7D%20where%20L%20where%20I%3C:Tuple%7BVararg%7BUnion%7BInt64,%20AbstractRange%7BInt64%7D,%20Base.AbstractCartesianIndex%7D,N%7D%20where%20N%7D%20where%20A%3C:Union%7BBase.ReinterpretArray%7BT,N,S,A%7D%20where%20S%20where%20A%3C:Union%7BSubArray%7BT,N,A,I,true%7D%20where%20I%3C:Union%7BTuple%7BVararg%7BReal,N%7D%20where%20N%7D,%20Tuple%7BAbstractUnitRange,Vararg%7BAny,N%7D%20where%20N%7D%7D%20where%20A%3C:DenseArray%20where%20N%20where%20T,%20DenseArray%7D%20where%20N%20where%20T,%20Base.ReshapedArray%7BT,N,A,MI%7D%20where%20MI%3C:Tuple%7BVararg%7BBase.MultiplicativeInverses.SignedMultiplicativeInverse%7BInt64%7D,N%7D%20where%20N%7D%20where%20A%3C:Union%7BBase.ReinterpretArray%7BT,N,S,A%7D%20where%20S%20where%20A%3C:Union%7BSubArray%7BT,N,A,I,true%7D%20where%20I%3C:Union%7BTuple%7BVararg%7BReal,N%7D%20where%20N%7D,%20Tuple%7BAbstractUnitRange,Vararg%7BAny,N%7D%20where%20N%7D%7D%20where%20A%3C:DenseArray%20where%20N%20where%20T,%20DenseArray%7D%20where%20N%20where%20T,%20SubArray%7BT,N,A,I,true%7D%20where%20I%3C:Union%7BTuple%7BVararg%7BReal,N%7D%20where%20N%7D,%20Tuple%7BAbstractUnitRange,Vararg%7BAny,N%7D%20where%20N%7D%7D%20where%20A%3C:DenseArray%20where%20N%20where%20T,%20DenseArray%7D%20where%20N%20where%20T,%20DenseArray%7D%7D%20where%20%23s4%3C:Union%7BComplex%7BFloat32%7D,%20Complex%7BFloat64%7D,%20Float32,%20Float64%7D%7D defined only for square matrices whereas matlab (and presumably octave) applies it element wise https://se.mathworks.com/help/matlab/ref/exp.html.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/julia/issues/32899#issuecomment-604643957, or unsubscribe https://github.com/notifications/unsubscribe-auth/AI2ARBYQ7ZPG2ZEE6WFX5FTRJOVSBANCNFSM4ILZEQXQ .

[ViralBShah]

It might help to post the code on discourse.julialang.org to see what you're trying to do. Julia is not the same as matlab/octave and a few things have to be done differently as it has a different underlying philosophy. It may appear complex if you took your matlab/octave code and just tried to run it in Julia. For example, as @fredrikekre pointed out, exp here is probably not doing what your intention is, and perhaps the data structures are coming out somewhat differently than you expect them to.

While in general one can always expect things to have better error messages, in this particular case, the error message is accurate.

[fab4ap]

OK - thanks. Probably worth my time to understand how to do it in Julia - now I know it is a valid error and not a bug/limitation.

On Fri, Mar 27, 2020 at 11:32 AM Viral B. Shah notifications@github.com wrote:

It might help to post the code on discourse.julialang.org to see what you're trying to do. Julia is not the same as matlab/octave and a few things have to be done differently as it has a different underlying philosophy. It may appear complex if you took your matlab/octave code and just tried to run it in Julia. For example, as @fredrikekre https://github.com/fredrikekre pointed out, exp here is probably not doing what your intention is, and perhaps the data structures are coming out somewhat differently than you expect them to.

While in general one can always expect things to have better error messages, in this particular case, the error message is accurate.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/julia/issues/32899#issuecomment-605064407, or unsubscribe https://github.com/notifications/unsubscribe-auth/AI2ARBYT525333L3OVS4XADRJTBJTANCNFSM4ILZEQXQ .

[StefanKarpinski]

If you want to do elementwise exp then you should write exp.(A) instead of exp(A).

[fab4ap]

I tried that already. But still I get: ERROR: LoadError: MethodError: no method matching exp(::Transpose{Float64,Array{Float64,1}}) Closest candidates are: exp(!Matched::Float16) at math.jl:1114 exp(!Matched::Complex{Float16}) at math.jl:1115 exp(!Matched::Missing) at math.jl:1167

Oh well. I'll probe into it when I have time. Thanks

On Fri, Mar 27, 2020 at 1:49 PM Stefan Karpinski notifications@github.com wrote:

If you want to do elementwise exp then you should write exp.(A) instead of exp(A).

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/julia/issues/32899#issuecomment-605151865, or unsubscribe https://github.com/notifications/unsubscribe-auth/AI2ARB23WKNW7PKHQ7E7OGTRJTRJFANCNFSM4ILZEQXQ .

[StefanKarpinski]

Are you sure?

julia> using LinearAlgebra

julia> A = rand(3, 3)
3×3 Array{Float64,2}:
 0.359941  0.298302  0.682141
 0.410766  0.461424  0.771319
 0.324901  0.890749  0.184061

julia> exp(transpose(A))
ERROR: MethodError: no method matching exp(::Transpose{Float64,Array{Float64,2}})
Closest candidates are:
  exp(::Complex{Float16}) at math.jl:1131
  exp(::Missing) at math.jl:1183
  exp(::BigFloat) at mpfr.jl:604
  ...
Stacktrace:
 [1] top-level scope at REPL[12]:1

julia> exp.(transpose(A))
3×3 Array{Float64,2}:
 1.43324  1.50797  1.38389
 1.34757  1.58633  2.43695
 1.97811  2.16262  1.20209

The error message you're reporting doesn't make sense for broadcasted exp unless you have a vector of transposed matrices, which seems unlikely.

[mbauman]

@fab4ap's original data structure is a vector of transposed vectors — so performing exp element wise on the outer container is going to apply exp to each transposed vector. You need to go down two levels; broadcast only goes inside one container.

We should probably take this to discourse as we're fairly far afield from the original problem.

[fab4ap]

Super apologies. I revisited it at leisure and discovered that I had mistakenly used Wlen = [0.0:0.01:20].1e-6; instead of Wlen = (0.0:0.01:20).1e-6; When I corrected that and also did element-wise subtraction with the .- operator, all was well. Again my apologies - will be more careful before posting in the future. Here is the correct snippet:

c = 3.74210^8; # Speed of light in vacuum in m/sec h = 6.62510^-34; # Planck's constant bk = 1.3810^-23; # Boltzmann constant T = [475, 625, 700]; # Temperature in Kelvin Wlen = (0.0:0.01:20).1e-6; # Wavelengths

= Planck's Blackbody Radiation Law

    2*h*c^2              1

PL = ------- ----------------- Wlen^5 exp(hc/WlenbkT) - 1 =# for i=1:3 PL=(2hcc)./( (Wlen.^5).( exp.((h.c)./(bk.T[i].*Wlen)) .- 1) ); # Planck's law

...

end

On Fri, Mar 27, 2020 at 4:13 PM Matt Bauman notifications@github.com wrote:

@fab4ap https://github.com/fab4ap's original data structure is a vector of transposed vectors — so performing exp element wise on the outer container is going to apply exp to each transposed vector. You need to go down two levels; broadcast only goes inside one container.

We should probably take this to discourse as we're fairly far afield from the original problem.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/julia/issues/32899#issuecomment-605295139, or unsubscribe https://github.com/notifications/unsubscribe-auth/AI2ARBZWUNIFLB4P73OT22LRJUCHTANCNFSM4ILZEQXQ .

[jonniedie]

What ever happened to this? exp(A::AbstractSparseMatrix)--and quite a few other matrix types--still seem not to be implemented. Here are two I noticed:

julia> using SparseArrays

julia> exp(sprand(3, 3, 0.5))
ERROR: MethodError: no method matching exp(::SparseMatrixCSC{Float64, Int64})
.
.
.

julia> exp(ones(Int, 3, 3) .// 1)
ERROR: MethodError: no method matching exp(::Matrix{Rational{Int64}})
.
.
.

It seems like it would be useful to have a fallback to Matrix{Float64} for these cases.

[jonniedie]

Any reason we can't use something like this from StaticArrays as a fallback for general matrices?

[birkmichael]

What if we port this code for fast sparse exponentiation from Matlab to Julia?

[stevengj]

I'm skeptical of the idea of automatically converting to a dense matrix as a fallback. (We don't even do this for inv, much less exp.) Converting a sparse matrix to a dense one is really something that the user should decide explicitly, as @antoine-levitt commented above.

@birkmichael, that code has already been ported to Julia, by the original authors no less: FastExpm.jl. But it's not really a general-purpose solution from what I can tell — whether it preserves sparsity, or whether it converges at a reasonable rate, seems to be highly problem-dependent. A standard-library default needs to be bulletproof.

Basically, if you have a sparse matrix and you want to exponentiate it, you really need to make an informed decision for your problem — there isn't a great default choice for Julia to make.

1. Is the matrix small and you don't care too much about performance? Then convert it to dense (or use a dense matrix to start with).
2. Do you just need to multiply exp(A) by a vector? Then consider an iterative Krylov method, e.g. from ExponentialUtilities.jl or ExponentialAction.jl or KrylovKit.jl
3. If you need the whole exp(A) matrix, try out FastExpm.jl … or maybe other algorithms in future packages … and see if they work well for you.

[Krastanov]

This discourse discussion on LinearSolve and A\b syntax echos some of the discussion here: https://discourse.julialang.org/t/a-b-vs-linearsolve-jl-at-juliacon-2022/85030

In summary, both exp(nondense_matrix_type) and \ might need to involve an informed decision from the user on what is the appropriate algorithm for the matrix type. On the side of \ it seems the community has settled on "do something fast and maybe naive in Base and have better, more general solutions in LinearSolve". On the side of exp the current sentiment seems to be that it is better to force the user to make an informed decision. There are some other differences between the two situations but they seem inconsequential.

Maybe a moderator can close this issue (I am the one that opened it, but I do not want to overstep) and we can leave the commentary here for future Julia historians. It can be reopened in the future if there is someone with stronger opinion.