`expm` missing method definitions, result type inconsistent

[jakebolewski]

There is some inconsistency with the output type for special matrix types.
expm(::Diagonal) -> Diagonal
expm(::Union(Symmetric,Hermitian)) -> Matrix

expm is currently undefined for the following AbstractMatrix subtypes:

Bidiagonal
UpperTriangular
LowerTriangular
SymTridiagonal

[andreasnoack]

Regarding the first part, I think it is covered by #18. For the second part, I'm not sure there is an issue. It is not obvious that we would ever want to define expm(SymTridiagonal).

[jiahao]

expm of a Triangular matrix is Triangular, so it would make sense to define that (and by extension, Bidiagonal too). But the others would in general produce dense matrices.

[tkelman]

I'm not aware of a good algorithm for sparse matrix exponential, the result will generally be dense as Jiahao said. Does SciPy implement these? I don't think we have a good general matrix-free operator abstraction in Base, but there are a few experiments out in packages. You probably don't want the full matrix exponential in these special cases, but an operator that computes the action of the matrix exponential applied to a vector would be useful and can be done efficiently. https://github.com/acroy/Expokit.jl ?

[argriffing]

I'm not aware of a good algorithm for sparse matrix exponential, the result will generally be dense as Jiahao said. Does SciPy implement these?

I've worked a bit on the scipy matrix exponential functions, and I agree there is not much point in a sparse matrix exponential except for certain structured matrices (e.g. (block) triangular or (block) diagonal).

I hope Julia gets a good matrix-free linear operator abstraction, if it doesn't have this already. The one in scipy is a bit of a mess in my opinion. For the combined expm(M)*v function, as opposed to just expm(M), it is very helpful to use the sparsity of M (presumably through the matrix-free linear operator abstraction).

If Julia has a nice linear operator abstraction and is able to work with the Higham lab and other historically MATLAB or NAG based labs in terms of licensing and implementations then I'll probably switch from Python to Julia for scientific computing.

[jiahao]

I hope Julia gets a good matrix-free linear operator abstraction,

There are several ideas floating around in IterativeSolvers.jl, KrylovMethods.jl, LinearMaps.jl, and LinearOperators.jl. I've been meaning to find time to sit down and look at all of these and find common ground.

[argriffing]

I've been meaning to find time to sit down and look at all of these and find common ground.

I would be excited to see this.