Solving system of linear equations using Cython

[Shriramana Sharma]

Hello. I wrote a Python function which involves solving a system of
linear equations. Currently I'm using SymPy's Matrix LUsolve, but
SymPy is pure Python. I was thinking of porting the function (or even
module) to Cython for speed. Do I understand that for solving a system
of linear equations in Cython, it is best to use Cython's link with
NumPy? http://docs.cython.org/src/tutorial/numpy.html hints at this,
but doesn't provide an example exactly with linear equations, hence my
question.

Thank you.

--
Shriramana Sharma

[Shriramana Sharma]

On Wed, Nov 14, 2012 at 9:15 PM, Shriramana Sharma <sam...@gmail.com> wrote:
> http://docs.cython.org/src/tutorial/numpy.html hints at this,
> but doesn't provide an example exactly with linear equations, hence my
> question.

Sorry for the additional mail -- I didn't voice my requirement very
clearly. I'd like to see a minimal example of something like the
following done in Cython:

from numpy import matrix,linalg
A = matrix([[2,2,3],[11,24,13],[21,22,46]])
B = matrix([[15],[98],[203]])
X = linalg.solve(A,B)

Thanks!

--
Shriramana Sharma

[Bradley Froehle]

Hi Shriramana:

There isn't much benefit to doing this in Cython.  NumPy will (hopefully) use an optimized BLAS/LAPACK library to do the matrix solve, and that is where you will likely spend all of your time anyway.

-Brad

[Shriramana Sharma]

On Wed, Nov 14, 2012 at 10:16 PM, Bradley Froehle
<brad.f...@gmail.com> wrote:
>
> There isn't much benefit to doing this in Cython. NumPy will (hopefully)
> use an optimized BLAS/LAPACK library to do the matrix solve, and that is
> where you will likely spend all of your time anyway.

Oh -- thanks for that! But I'll probably use Cython for some other parts
of my project as well. Will get back later.

--
Shriramana Sharma

[Dag Sverre Seljebotn]

On 11/14/2012 05:46 PM, Bradley Froehle wrote:
> Hi Shriramana:
>
> There isn't much benefit to doing this in Cython. NumPy will
> (hopefully) use an optimized BLAS/LAPACK library to do the matrix solve,
> and that is where you will likely spend all of your time anyway.

It's worth stressing that this is actually an understatement -- unless
you know *a lot* about your CPU and write code tuned specifically for
it, anything one implements oneself in Cython will be 10x to 100x slower.

Dag Sverre

[David Warde-Farley]

To add to what Dag said, if Python overhead is a significant problem (i.e. you have lots of 3x3 or 4x4 matrices, and a lot of time is spent in stupid error-checking and object unboxing by python/numpy)  it is possible to call BLAS and LAPACK directly from Cython.

numpy.show_config() should be able to tell you what BLAS and LAPACK libraries NumPy is using on your system, so that you can link against them directly.

On the other hand, if your matrices are sufficiently large that time inside spent inside LAPACK functions absolutely dwarfs everything else, Cython will not buy you anything. You'd be better off making sure that your BLAS library is a high-quality multithreaded implementation such as a tuned ATLAS or the Intel Math Kernel Library.

These materials from a tutorial Dag gave might help if you have the "lots of small matrices" problem:

http://sage.math.washington.edu/home/wstein/www/home/dagss/cython-notur09/

[Robert Bradshaw]

On Wed, Nov 14, 2012 at 1:06 PM, David Warde-Farley
<d.warde...@gmail.com> wrote:
> On Wed, Nov 14, 2012 at 2:33 PM, Dag Sverre Seljebotn
> <d.s.se...@astro.uio.no> wrote:
>>
>> On 11/14/2012 05:46 PM, Bradley Froehle wrote:
>>>
>>> Hi Shriramana:
>>>
>>> There isn't much benefit to doing this in Cython. NumPy will
>>> (hopefully) use an optimized BLAS/LAPACK library to do the matrix solve,
>>> and that is where you will likely spend all of your time anyway.
>>
>>
>> It's worth stressing that this is actually an understatement -- unless you
>> know *a lot* about your CPU and write code tuned specifically for it,
>> anything one implements oneself in Cython will be 10x to 100x slower.
>
>
> To add to what Dag said, if Python overhead is a significant problem (i.e.
> you have lots of 3x3 or 4x4 matrices, and a lot of time is spent in stupid
> error-checking and object unboxing by python/numpy) it is possible to call
> BLAS and LAPACK directly from Cython.

Those libraries are primarily optimized for larger matrix
computations; if you really have 3x3 matrices I wouldn't be surprised
if you could write a faster specialized implementation. (I could be
wrong though.)

> numpy.show_config() should be able to tell you what BLAS and LAPACK
> libraries NumPy is using on your system, so that you can link against them
> directly.
>
> On the other hand, if your matrices are sufficiently large that time inside
> spent inside LAPACK functions absolutely dwarfs everything else, Cython will
> not buy you anything. You'd be better off making sure that your BLAS library
> is a high-quality multithreaded implementation such as a tuned ATLAS or the
> Intel Math Kernel Library.

+1

> These materials from a tutorial Dag gave might help if you have the "lots of
> small matrices" problem:
>
> http://sage.math.washington.edu/home/wstein/www/home/dagss/cython-notur09/

Note that a lot has changed since those slides were written, though
they're still a good resource.

- Robert

[Dag Sverre Seljebotn]

Or OpenBLAS!:

https://github.com/xianyi/OpenBLAS

A lot simpler to compile than ATLAS, and it can outperform MKL in some
cases.

Dag Sverre

[Sturla Molden]

On 14.11.2012 22:24, Dag Sverre Seljebotn wrote:

> Or OpenBLAS!:
>
> https://github.com/xianyi/OpenBLAS
>
> A lot simpler to compile than ATLAS, and it can outperform MKL in some
> cases.

Yes, ATLAS is a PITA to build, and slower than OpenBLAS.

After compiling OpenBLAS then compile against LAPACK from Netlib. It is
just a makefile that needs no configuration (except where to locate BLAS).

ACML is good too, and no compilation is needed.

To avoid dependency on Fortran ABI (and get rid of pesky Fortran 77 work
arrays, etc.), we can use LAPACKE from Cython instead of calling LAPACK
directly.

http://www.netlib.org/lapack/lapacke.html


Sturla