Coverage remained the same when pulling df566b3ff30f589a2745371a4d376ec955b87389 on argriffing:alex-leach-expokit-wrapper into 70fd6ff on scipy:master.

The interface should IMHO be expm(A, method="expokit") and expm(A, method="higham") rather than adding a submodule for one new function.

Cleaned up + rebased branch here: https://github.com/pv/scipy-work/tree/pr-2974
The following commands update the current branch to that point

git remote add pv https://github.com/pv/scipy-work.git
git fetch pv
git reset --hard pv/pr-2974

With the other PR with the exact norm estimation, your implementation of Higham's algorithms seem to be within 20% of Expokit speed. Not bad :)

The Krylov estimate part still needs some love; it shouldn't call todense(), and this seems to be a bit light on tests.

[EDIT: fixed the git stuff]

I'm probably not up for writing sparse matrix Fortran<->SciPy code if that's what would be required to make the EXPOKIT sparse expm_multiply work, but I might be up for a more detailed investigation of the accuracies of the expm implementations, which Alex Leach has begun in his PR thread.

Regarding the organization of the code, I would be OK with either expm(A, method=foo) or with a separate expokit submodule. My reasoning for tolerating the submodule is that EXPOKIT has an interface that potentially provides many functions related to expm, ode, exponential integrators, etc., if all of these were to be wrapped. It also seems somewhat standard, as it is prominent in Burkardt's comparison of expm implementations. Furthermore, the EXPOKIT author Roger Sidje seems to be not only alive, but also actively communicating with and helping the scipy devs and is using his EXPOKIT in new publications. Also, I believe that EXPOKIT has expm functions specialized for rate matrices (square matrices with non-negative off-diagonals, and rows summing to zero) which are the usual matrices that I personally use for expm in my work, but I think these more specialized functions are not yet wrapped for scipy.

@argriffing: use git push -f. "non-fast-forward" means you have an all-new patch set, which is not derived from the previous ones.

@pv thanks, your commits are now added into this PR.

Coverage remained the same when pulling 0447188 on argriffing:alex-leach-expokit-wrapper into 57ac922 on scipy:master.

Looking more closely into the accuracy differences, I see that for medium-sized (n=20) complex matrices constructed by taking the logm of matrices whose entries are independently standard normally distributed, the higham expm always seems to give closer (in norm) answers than does the expokit expm. For other matrix distributions and other ways of comparing accuracy, the higham expm still seems to be usually better, but not for every matrix. These comparisons still have the caveat that the accuracy calculations depend on the accuracy of logm, and the comparisons extend the earlier comparisons by @alexleach.

expm accuracy testing code
https://gist.github.com/argriffing/6978763

linalg.expm vs expokit expm for logm of entrywise standard random normal matrix
https://gist.github.com/argriffing/6978695

linalg.expm vs expokit expm for logm of entrywise standard uniform random matrix
https://gist.github.com/argriffing/6978754

I think the things that need to be decided before merging are:

• What to do with the Krylov stuff. We certainly want tol be tunable by the user. This part needs tests, too --- I'm not 100% sure the code currently there actually works...
• Taking care of the t argument and other extra functionality of expokit. Adding t= and structure= keywords to expm could be one option. For the higham mode, t would just pre-multiply, and structure would have no effect.

The second option to stuffing expm full of keyword args is to add more specialized functions for each special task.

closing https://mail.scipy.org/pipermail/numpy-discussion/2016-January/074810.html

maybe tag these with something like might-resurrect or something to this effect?

Chuck suggested inactive --- added that.