[SciPy-User] faster expm

[Jonathan Helmus]

Alex this is some very nice work.  Unfortunately I don't think the 
Expokit licensed (http://www.maths.uq.edu.au/expokit/copyright) is 
compatible with the BSD license SciPy uses. Specifically approval must 
be obtain before commercial use.  Including the package in the Topical 
Software list and/or the SciPy Cookbook/SciPy Central might be prudent.

     - Jonathan Helmus

On 11/07/2012 09:47 AM, Alex Leach wrote:
>
> 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
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20121107/09583dff/attachment.html>