Add Julia matrix exponential?

[sdewaele]

I needed a pure Julia implementation of matrix exponentiation, so I updated an old 2012 gist to Julia 1.0, see here. I am happy to contribute it to this repo if desired. Probably a better implementation can be made that is closer to the current Julia implementation. After all, it only uses LAPACK.gebal! and LAPACK.gesv! as far as I can tell.

[andreasnoack]

That would be great. There is absolutely no reason why the base implementation uses gesv! directly. Instead of LAPACK.gesv!(V-U, X) is should just do ldiv!(lu!(V-U), X). It would be great if you could prepare such a PR for for the julia repo. I'd also be happy to accept a PR with a generic gebal implementation in this repo.

[sdewaele]

Thanks for the pointer for the replacement of LAPACK.gesv, that is a good start. I have made this replacement and did some performance testing on it. I do find some increase in execution time, typically around 1% or so. I get occasional larger increases, up to ~30%, which seem to occur when there is a competing process running. This would all be fine if it would give us a fully generic matrix exponential, however, there is still the call to LAPACK.gebal!. In all, I don't think it is worthwhile submitting a PR for the Julia repo. It appears that creating a generic version for gebal is not that simple.

[sdewaele]

I created a generic matrix exponential based on the current Julia matrix exponential, but without matrix balancing, see here for the added code. I think it is valuable to have this as a starting point. Please let me know if you would like to include it in this repository.

[sdewaele]

Is there an interest in adding the generic matrix exponential?

[andreasnoack]

Yes. Did you get around looking into a generic gebal or is your version still not balancing?

[timholy]

We are quite interested in this. @sdewaele, your link for your code is broken, is there some place we can see it and perhaps help get this done?

[sdewaele]

I am sorry to say I have deleted my fork of the repo.

Anyway, not too much is lost: all I did compared the the current implementation is replace LAPACK.gesv!(V-U, X) by ldiv!(lu!(V-U), X). and skipped balancing.