[SciPy-User] faster expm
Alex Leach
beamesleach@gmail....
Wed Nov 7 08:47:51 CST 2012
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-0001.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-0001.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-0001.py
More information about the SciPy-User
mailing list