Roadmap

[mforets]

as a step forward what do you think about:

• add more functionality:
  • phimv Add phimv #11
  • padm (a preliminary version is here, and some results here) Add padm.jl #14
  • chbv Add chbv #15
• make docs
• register the package
• benchmarks Benchmarks #21

CC: @acroy, @garrison

[acroy]

Sounds great! I would probably add some benchmarks, but this could also be done later on.

[mforets]

cool. let me add those new issues.

[Jutho]

How about giving these functions a bit more of a descriptive name, in line with Julia's style conventions? E.g. using a keyword argument for the algorithm or so.

[mforets]

+1. suggestion:

• padm(A) -> expm(A, algorithm="pade") (by extending Base.expm)
• expmv(A, v) -> expm_vector(A, v)

[acroy]

I am not sure about that; expmv, phimv and padm are calculating different things (i.e. solve different ODEs) in different ways. At least expmv (or expv) is used in other languages as well and thus people "know what they get". Extending Base.expm is probably reasonable.

How about having a "matrix exponential" type? Has that been discussed already?

[Jutho]

I agree that a suitable design for the interface is non-trivial. You cannot extend Base.expm by just a keyword argument.

A matrix exponential type sounds like a nice Julian way to go. You can then define A_mul_B!(w, expA, v; alg = ...) for algorithms that only act on a vector, and functions like full(expA) for algorithms that compute the full exponential.

[acroy]

Ah, you are right of course about Base.expm (apparently, I got a bit rusty).

The behaviour you describe is exactly what I had in mind. I guess the new type (say MatrixExp) would be a subtype of AbstractMatrix or is there any catch (partly having in mind LinearMaps.jl)?

[Jutho]

Good question. I guess an AbstractMatrix does require that you can access individual entries more or less efficiently, though there is no strict definition I think. It could be implemented similarly as a LinearMap, or even as a specific subtype of it (but then you need to add a dependency, and that's something I am also always rather hesitant about).

[acroy]

I would like to avoid dependencies at this stage. Later on one could still tie into your package or whatever.

Access to individual entries could be slow. (One can get an entire column via exp(A)*[...,0,1,0,...] etc.) In the manual I found the statement

The AbstractArray type includes anything vaguely array-like, and implementations of it might be quite different from conventional arrays. For example, elements might be computed on request rather than stored. [...] It is recommended that these operations [getindex and setindex] have nearly constant time complexity, or technically Õ(1) complexity, as otherwise some array functions may be unexpectedly slow.

I am not sure if the latter can be guaranteed in our case.

[Jutho]

That's what I am also thinking; they shouldn't be subtypes of AbstractArray.

[ChrisRackauckas]

Thanks for reviving this project! We are creating some exponential integrators and hit a wall since this area needs some Julia development.

expmv is a sufficiently well-known name (it's the name of the Higham algorithm at least) that it shouldn't be avoided. phimv is related and sufficiently clear that it's the phi matrix vector product. I would avoid the unnecessary allocation though and make the Pade part algorithm=:pade.

Like matrix exponentials and all of that, it should be a function and the matrix exponential type which has A_mul_B! call expmv should be added separately so that way users could easily add this kind of overload to their own abstract linear maps.

Does anyone see a license at all for ExpoKit? I can't find it, so I'd be very careful and not look at their code.

It would be nice if the Higham algorithm was offered as well (http://epubs.siam.org/doi/abs/10.1137/100788860). I would like it if there was sufficient ability in one package to pull off this kind of benchmarking (http://profs.sci.univr.it/~caliari/pdf/preCKOR13.pdf).

Anyways, if this project is revived (looks like it just was a few days ago!) we will likely abandon our efforts over at SciML/Expmv.jl#1 and use this.

[acroy]

@ChrisRackauckas : I was also very happy that @mforets woke me from my "hibernation" ...

Regarding the license: Before I started this package few years ago I was in contact with R.B. Sidje and asked him about it. It seems there is no strict license and he is Ok with us using his code as a starting point.

The Higham algorithm was actually implemented by @marcusps in his package Expmv.jl. He also did some benchmarking, but it seems not much happened recently. There is also #1 on having the algorithm in this package.

[ChrisRackauckas]

Expmv.jl allocates a ton (no in place form either) and doesn't have the norm 1 estimation function, and it has seemed that development stopped, so I was hoping you guys would pick it up. Looks like that is on the list as #1

[mforets]

that sounds goood to me :)
within 4 weeks or so i propose to focus on adding functionality and benchmarks. the interface for inclusion in JuliaDiffEq can be revised afterwards.

[ChrisRackauckas]

Well since you gave me a timeframe I'll hold you to it 😄