Julia: `exp` does not work on many sparse `LinearAlgebra` matrices
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.
Krastanov
All 22 comments
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.
mbauman
on 14 Aug 2019
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.
Krastanov
on 15 Aug 2019
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.
mbauman
on 15 Aug 2019
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
on 15 Aug 2019
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.
Krastanov
on 15 Aug 2019
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...
antoine-levitt
on 15 Aug 2019
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.
Krastanov
on 15 Aug 2019
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.
antoine-levitt
on 15 Aug 2019
@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.
Krastanov
on 15 Aug 2019
FWIW, @Krastanov, @antoine-levitt, my research group has encountered a similar problem related to quantum mechanics involving sparse matrices.
ryanengle
on 8 Oct 2019
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.
fab4ap
on 26 Mar 2020
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
on 26 Mar 2020
If exp(A') is the same as exp(A)' that's an easier rule to add that may work in more cases.
ViralBShah
on 26 Mar 2020
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.
fredrikekre
on 26 Mar 2020
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
.
fab4ap
on 27 Mar 2020
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.
ViralBShah
on 27 Mar 2020
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
.
fab4ap
on 27 Mar 2020
If you want to do elementwise exp then you should write exp.(A) instead of exp(A).
StefanKarpinski
on 27 Mar 2020
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
.
fab4ap
on 27 Mar 2020
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.
StefanKarpinski
on 27 Mar 2020
@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.
mbauman
on 27 Mar 2020
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
.
fab4ap
on 31 Mar 2020
Related issues
tkoolen
·
3Comments
wilburtownsend
·
3Comments
m-j-w
·
3Comments
StefanKarpinski
·
3Comments
StefanKarpinski
·
3Comments
Most helpful comment
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
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: