John Cremona posted the following message to the sage development list about matrix exponentials. I'm copying it to here since it asks about the scipy matrix exponential method (we say numpy below, but we really mean scipy...) John Cremona wrote:
I have just been to a colloquium talk by numerical analyst Nick Higham (Manchester) called "How to compute and not to compute a matrix exponential". He has new methods which are now in mathematica, matlab and NAG but (apparantly) nowhere else. He only seemed interested in getting good speed & precision to 16 decimals but (when I asked) confirmed that the methods should apply to give arbitrary precision.
I just checked and see that Sage's matrix exp() uses something stupid except over RDF/CDF where it uses a pade approximation method via numpy. The method of the talk was a variant of that, the main trick being to use exactly the right order of Pade approx. so maximise precision and speed.
I would like to know how good the numpy method is, and whether it can be improved to this "state of the art" version at least for RDF. Then it could be another selling point for Sage.
Could you CC the numpy devlist as well on this? It sounds exciting!
I will if you give me the address (or you can perhaps?). It might be worth including Higham's URL: http://www.maths.manchester.ac.uk/~higham/ as he has lots of his talks up there including some which are similar to the one I heard.