faster expm

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

faster expm

Bugzilla from beamesleach@gmail.com

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
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

expokit.f (162K) Download Attachment
time_expm.py (2K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 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
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Bugzilla from beamesleach@gmail.com

On Wednesday 07 Nov 2012 10:04:55 Jonathan Helmus wrote:

> 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

 

(sorry, replied directly before)

 

Hi Jonathan,

 

Good spot about the license issue. That's a shame, but seems like the author might give SciPy permission to use it, if asked..

 

Including it as a package, documented build method, or otherwise would be cool. If someone wants to go ahead and do that, that's fine; I don't think I have the access rights (or the time, really) to do any of that myself.

 

One other thing; just did a diff to see if I'd made any other changes. Looks like I also had to edit some calls to the matvec functions. Converted C code fails during compilation without them. Universal diff attached.

 

Cheers,

Alex


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

expokit.udiff (11K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Alan G Isaac-2
On 11/7/2012 10:57 AM, Alex Leach wrote:
> seems like the author might give SciPy permission to use it, if asked.


So, will you ask?

Just fyi, I have had good luck with such requests,
when I ask if BSD/MIT would be possible and explain
why it is needed for SciPy.

Alan Isaac

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Nathaniel Smith
In reply to this post by Bugzilla from beamesleach@gmail.com

On 7 Nov 2012 15:58, "Alex Leach" <[hidden email]> wrote:
>
> On Wednesday 07 Nov 2012 10:04:55 Jonathan Helmus wrote:
>
> > 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
>
>  
>
> (sorry, replied directly before)
>
>  
>
> Hi Jonathan,
>
>  
>
> Good spot about the license issue. That's a shame, but seems like the author might give SciPy permission to use it, if asked..

Note that just asking for SciPy to have permission is insufficient. We actually need it to allowed for use under bsd-like terms. This is because we need to make sure that everyone who uses scipy also can use it. The current license makes it impossible to use expokit and GPL code in the same program, eg.

-n

>  
>
> Including it as a package, documented build method, or otherwise would be cool. If someone wants to go ahead and do that, that's fine; I don't think I have the access rights (or the time, really) to do any of that myself.
>
>  
>
> One other thing; just did a diff to see if I'd made any other changes. Looks like I also had to edit some calls to the matvec functions. Converted C code fails during compilation without them. Universal diff attached.
>
>  
>
> Cheers,
>
> Alex
>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Bugzilla from beamesleach@gmail.com

On Thursday 08 Nov 2012 07:19:51 Nathaniel Smith wrote:

>

> Note that just asking for SciPy to have permission is insufficient. We

> actually need it to allowed for use under bsd-like terms. This is because

> we need to make sure that everyone who uses scipy also can use it. The

> current license makes it impossible to use expokit and GPL code in the same

> program, eg.

>

 

Kind email sent. Fingers crossed.

 

Alex


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Bugzilla from beamesleach@gmail.com
In reply to this post by Nathaniel Smith

Great news guys. The Expokit author, Professor Roger Sidje, is perfectly happy for Expokit to be used in SciPy. He just likes to keep a tab on things, for appraisals and such like.

 

His full response, and my original message, is below:-

 

Alex,

 

Thanks for your interest.

 

Please note that the wording of the license is merely a way to keep

track of where it is used rather than an attempt to restrict its

usage.

 

The license is simple -- any other open/closed source can use the

software any way they want.

 

All the license does is to allow me to know/have a bullet point in

those end-of-year activity reports for the powers that be, telling

that the software is found useful.

 

I have had similar requests in the past re:GPL. All I say is, feel

free to use/modify the software any way you want. This e-mail of yours

is exactly the kind of expression-of-interest indicated (even for use

in a commercial setting).

 

Because I have moved from UQ (in Australia) to UA (in the USA), making

changes to Expokit won't be convenient for a while. Please consider

this response as a plain acknowledgment that Expokit can be

embedded/released under GPL. You can add the following in the meantime:

 

"Alternatively, EXPOKIT may be used under the terms of the GPL."

 

By so doing, it would mean that you have chosen to branch on the

GPL-alternative, while allowing others to do what they want on the

other alternatives (and it is understood that, per GPL, if someone

ever gets the GPL alternative or its derivatives, they will abide by

GPL).

 

In particular, possible future versions of Expokit not derived from

your copy will remain free to be licensed any way (which is an option

that commercial packages want -- even though they use EXPOKIT without

royalties).

 

Thanks

---

RBS

 

On Thu, Nov 8, 2012 at 5:44 AM, Alex Leach <[hidden email]> wrote:

> Dear Professor Sidje,

>

> I hope this email finds you well.

>

> As way of introduction, I am a PhD student studying Computational Biology at

> the University of York. I have recently used Expokit for modelling

> polymerisation in a specific chemical reaction, and had some good

> experiences with it.

>

> I thought I'd share my joy with the SciPy community, who have implemented

> similar algorithms for computing matrix exponentials, but in Python. After

> writing a small wrapper function (and making some very minor changes to the

> Expokit source code), Expokit's dgpadm worked as a drop-in replacement for

> their `expm` implementation, but was about twice as fast.

>

> The exact changes I made can be seen in expokit.udiff, which I sent as an

> attachment with this email to the SciPy list. (The SciPy mailing system made

> the nasty file name: attachment-0001.bin.)

>

> Impressed with the speed improvements, the SciPy community would be

> interested in using Expokit in precedence of their original implementation.

>

> However, there is an issue with licensing... Python is GPL source code, and

> SciPy uses a BSD license. I'm sure you know more about these licenses than I

> do, but from what I gather, the Expokit Copyright is restrictive for

> corporations who would wish to use the software for commercial purposes,

> whereas the BSD license is non-restrictive in this sense.

>

> So, as you've probably guessed, I'm writing to ask your kind permission, on

> behalf of the SciPy community, to include and distribute the Expokit source

> code, as a part of the SciPy source code distribution.

>

> With your permission, expokit.f would be included in future SciPy releases

> and built as a Python shared module, which is compiled from intermediate C

> code that includes the relevant Python headers, etc. The first link in this

> message, above, shows the performance improvements I achieved, build

> instructions and python test module, detailing how I did this.

>

> If you have any questions or requirements with regard to the usage of

> expokit, please don't hesitate to ask.

>

> There would be some work to merge Expokit into SciPy, but there seems to be

> a lot of interest in doing this, so if you give the green light, the

> community would get this done in no time, I reckon. It's impossible to tell

> how many tools would benefit from this change downstream, but I know that

> some popular protein visualisation tools use SciPy internally, so I think

> this would benefit a massive user-base.

>

> I look forward to hearing back from you.

> Kind regards,

> Alex

>

>

> --

> Alex Leach BSc. MRes.

> Department of Biology

> University of York

> York YO10 5DD

> United Kingdom

> www: http://bioltfws1.york.ac.uk/~albl500

> EMAIL DISCLAIMER: http://www.york.ac.uk/docs/disclaimer/email.htm


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Thomas Kluyver-2
On 8 November 2012 23:19, Alex Leach - Gmail <[hidden email]> wrote:
Great news guys. The Expokit author, Professor Roger Sidje, is perfectly happy for Expokit to be used in SciPy. He just likes to keep a tab on things, for appraisals and such like.

It certainly sounds like he's happy to see the code reused, but I'm not convinced that his e-mail technically provides the necessary permission:

> All I say is, feel free to use/modify the software any way you want.

This statement appears to grant a very permissive license, but it doesn't mention any right to redistribute the code, which is key.

> Please consider this response as a plain acknowledgment that Expokit can be embedded/released under GPL.

Unfortunately, we can't integrate GPL licensed code into Scipy, as Scipy is BSD licensed. To integrate that code, we'd need to get a similarly permissive license. I've just found this handy page: http://www.scipy.org/License_Compatibility

He might not be willing to give us that - he wants to ensure he knows about who's using his software, which is a requirement you can't embed in a permissive license. I think it's worth a try, though. A couple of points you could use:
- We're also keen to know who's using our software, although we don't make it compulsory for commercial users to tell us. E.g. IPython maintains a list of projects using IPython, which we use to demonstrate impact.
- He could grant a BSD license to just the bit of the code we're interested in, rather than all of Expokit.

Kudos for working on this,

Thomas

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Bugzilla from beamesleach@gmail.com

On 9 Nov 2012, at 00:57, Thomas Kluyver <[hidden email]> wrote:

Unfortunately, we can't integrate GPL licensed code into Scipy, as Scipy is BSD licensed. To integrate that code, we'd need to get a similarly permissive license. I've just found this handy page: http://www.scipy.org/License_Compatibility


All good.
Just got confirmation back; below.


On 9th Nov 2012, at 14:44, Roger B. Sidje wrote:

As I indicated the license is very permissive indeed.

Just substitute GPL with BSD in my earlier post.

You can branch with your BSD copy, make change and redistribute SciPy
with those changes any way you guys want, while allowing others to do
what they want with the other alternatives not derived from your
branch.
---
RBS

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: faster expm

Thomas Kluyver-2
On 9 November 2012 15:12, Alex Leach <[hidden email]> wrote:
As I indicated the license is very permissive indeed.

Just substitute GPL with BSD in my earlier post.

That sounds OK to me :-)

Thomas

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user