[SciPy-User] faster expm

[Alex Leach]

Sorry to pick up on this old thread, but I too was looking for a faster expm 
implementation, and had some joy. Thought I'd share...

The expokit library seems to be a highly regarded matrix exponentiation 
library (that Matlab probably uses), which is written in Fortran. I downloaded 
it a while ago and compiled it with f2py; first time I've used f2py and it 
worked pretty much out of the box!

The Expokit source code is available from:-
http://www.maths.uq.edu.au/expokit/download.html

Expokit just needs to be compiled and linked against LAPACK and BLAS 
libraries, so I've got two versions: one linked against NVIDIA's libcublas and 
another linked against Intel's libmkl (most people would link against blas and 
lapack libraries from ATLAS, I guess). The library can be specified manually, 
but it's probably easiest to just let numpy choose (like the command below 
does).

I've had it kicking around for a while, but just checked, and a compile 
command that just worked:-

f2py --fcompiler=intelem -c expokit.f -m expokit --link-blas_opt --link-
lapack_opt

Unless you use Intel's 64bit Fortran compiler, you'll want to remove or change 
the --fcompiler option.

This creates the Python shared module: expokit.so
This assumes a working numpy installation, and fortran compiler.

I've only needed to use Padé approximation myself, which is done with the 
subroutines [ZD]GPADM. These should give equivalent results to 
scipy.linalg.expm, but they take a lot more function arguments.

time_expm.py, attached, has an example wrapper function and runs some basic 
timings against scipy.linalg.expm.

The only change I made to the Fortran code was to add f2py style comments, so 
the output variables do return changed, and can be seen when inspecting with 
the numpy.info function. e.g.

>>> import numpy as np
>>> import expokit
>>> np.info(expokit.dgpadm)
dgpadm - Function signature:
  dgpadm(ideg,t,h,wsp,ipiv,iexph,ns,iflag,[m,ldh,lwsp])
Required arguments:
  ideg : input int
  t : input float
  h : input rank-2 array('d') with bounds (ldh,m)
  wsp : input rank-1 array('d') with bounds (lwsp)
  ipiv : input rank-1 array('i') with bounds (m)
  iexph : in/output rank-0 array(int,'i')
  ns : in/output rank-0 array(int,'i')
  iflag : in/output rank-0 array(int,'i')
Optional arguments:
  m := shape(h,1) input int
  ldh := shape(h,0) input int
  lwsp := len(wsp) in/output rank-0 array(int,'i')


For further information on the subroutine names and arguments, see the Expokit 
download link, above, which describes each subroutine, and provides separate 
links to each subroutine's source code. The Fortran source code is where to 
find the best documentation of each function and its arguments.

In terms of merging this into scipy, some work would need to be done...

	1) Wrapper functions would need to be created to conform to scipy's expm 
API; time_expm.py, attached, demonstrates python wrappers for dgpadm and 
zgpadm. Using PyFort to automatially write and compile C wrappers should give 
further, marginal performance improvements.
	2) Unit-tests on the wrapper functions.  I think scipy's expm unit-tests 
could be used here.
	3) The scipy build systems would need to be informed about it.
	4) scipy.linalg.__init__.py should probably try and import expokit wrapper 
functions over the expm, expm2 and expm3 functions. If that fails, fall back 
to the original scipy versions.

Performance improvements?
These are the results I got when just running it, averaging the times over 3 
runs. Note, both scipy and expokit use the same BLAS dot product functions, 
which does a lot of the heavy lifting; over multiple processor cores, I might 
add.

                                  ----- Mean time +- S.D. -----
Array size                Expokit                   scipy..expm              
25                        0.149 += 0.008            0.410 += 0.025           
50                        0.341 += 0.011            0.678 += 0.019           
75                        0.650 += 0.049            1.282 += 0.030           
100                       1.372 += 0.087            1.995 += 0.032           
125                       1.972 += 0.029            3.007 += 0.145           
150                       6.038 += 0.062            7.933 += 0.217           
175                       9.169 += 0.050            11.785 += 0.593          
200                       12.049 += 0.190           16.281 += 0.293

Hope this might be of help to someone.
Kind regards,
Alex
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20121107/8eea6928/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: expokit.f
Type: text/x-fortran
Size: 161471 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20121107/8eea6928/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: time_expm.py
Type: text/x-python
Size: 2935 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20121107/8eea6928/attachment.py>