speed question, numpy vs CDF
Pierre
I'm doing some computations with complex numbers and speed is
critical. I'm puzzled by the following.
sage: A= MatrixSpace(CDF, 2).random_element()
sage: B= MatrixSpace(CDF, 2).random_element()
sage: %timeit A*B
625 loops, best of 3: 20 µs per loop
sage: AA= numpy.array(A); BB= numpy.array(B)
sage: %timeit AA.dot(BB)
625 loops, best of 3: 2.24 µs per loop
So numpy seems much faster for this. But on the other hand
sage: z= A[0,0]
sage: %timeit z*z
625 loops, best of 3: 241 ns per loop
sage: zz= AA[0,0]
sage: %timeit zz*zz
625 loops, best of 3: 521 ns per loop
so now numpy is much slower!!
what is going on? is there a way to get the best of both worlds
simply, both for matrix operations and simple arithmetic?
i think zz above might still be considered as a 1 x 1 matrix instead
of a complex number, somehow, and this may be slowing things down.
any help appreciated!
thanks
pierre
Nils Bruin
> i think zz above might still be considered as a 1 x 1 matrix instead
> of a complex number, somehow, and this may be slowing things down.
number type is apparently a bit slower for individual element
arithmetic. It may well be that you're mainly measuring overhead,
though, so you should really test in a more representative situation
before committing to a particular implementation choice. numpy does
allow arbitrary types in its arrays. I doubt they're as optimized as
its own types, but you can try:
sage: A= MatrixSpace(CDF, 2).random_element()
sage: B= MatrixSpace(CDF, 2).random_element()
sage: %timeit A*B
sage: import numpy
sage: %timeit AA.dot(BB)
sage: AAA= numpy.array(A,dtype=type(A[0,0])); BBB=
numpy.array(B,dtype=type(B[0,0]))
sage: %timeit AAA.dot(BBB)
625 loops, best of 3: 2.33 µs per loop
sage: z=A[0,0]
sage: %timeit z*z
625 loops, best of 3: 101 ns per loop
sage: zz=AA[0,0]
sage: %timeit zz*zz
625 loops, best of 3: 253 ns per loop
sage: zzz=AAA[0,0]
sage: %timeit zzz*zzz
625 loops, best of 3: 107 ns per loop
sage: type(z); type(zz); type(zzz)
<type 'sage.rings.complex_double.ComplexDoubleElement'>
<type 'numpy.complex128'>
<type 'sage.rings.complex_double.ComplexDoubleElement'>
Pierre
arrays with CDF elements, which I didn't know was possible) is going
to be perfect for me.
Thank you so much!
Robert Bradshaw
With such small matrices (and elements), you're essentially measuring
overhead rather than arithmetic here. Of course if you have lots of
small matrices, that may be a relavant thing to measure. As the matrix
size grows, they should be the same, as multiplying CDF matrices
simply defers to multiplying numpy matrices.
sage: A= MatrixSpace(CDF, 200).random_element()
sage: B= MatrixSpace(CDF, 200).random_element()
sage: %timeit A*B
125 loops, best of 3: 7.31 ms per loop
sage: AA= numpy.array(A); BB= numpy.array(B)
sage: %timeit AA.dot(BB)
125 loops, best of 3: 7.34 ms per loop
- Robert
Pierre
i'm better off storing them as numpy matrices throughout... thanks for
your explanations though.
Pierre
On 13 fév, 18:32, Robert Bradshaw <rober...@math.washington.edu>
wrote:
William Stein
> I see. Well I *do* have hundreds of 2x2 matrices to multiply out so
> i'm better off storing them as numpy matrices throughout... thanks for
> your explanations though.
>
> Pierre
You might consider using Cython and writing a custom 2x2 matrix class.
It wouldn't be difficult... so I'll write one right now and respond
with the benchmarks.
-- William
> --
> To post to this group, send email to sage-s...@googlegroups.com
> To unsubscribe from this group, send email to sage-support...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-support
> URL: http://www.sagemath.org
--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org
Robert Bradshaw
> On Mon, Feb 13, 2012 at 9:59 AM, Pierre <pierre....@gmail.com> wrote:
>> I see. Well I *do* have hundreds of 2x2 matrices to multiply out so
>> i'm better off storing them as numpy matrices throughout... thanks for
>> your explanations though.
>>
>> Pierre
>
> You might consider using Cython and writing a custom 2x2 matrix class.
> It wouldn't be difficult... so I'll write one right now and respond
> with the benchmarks.
+1, exactly what I was going to consider. Depending on how naturally
it fits your problem, you could also consider packing your 2x2
matrices into larger arrays (e.g. representing n 2x2 matrices by a 4 x
n matrix and manually doing the multiplication) so you can your
computations in a "vectorized" form.
D. S. McNeil
sage: import numpy
sage: m = numpy.matrix([[1,2],[3,4]])
sage: m[:,0]
matrix([[1, 3]])
sage: m[:,int(0)]
matrix([[1],
[3]])
That is, if you use a Sage integer to index a numpy matrix, you don't
get the expected shape back. I know exactly why this happens --
http://trac.sagemath.org/sage_trac/ticket/10928 -- but ideally it
should be patched upstream.
Doug
William Stein
<robe...@math.washington.edu> wrote:
> On Mon, Feb 13, 2012 at 10:06 AM, William Stein <wst...@gmail.com> wrote:
>> On Mon, Feb 13, 2012 at 9:59 AM, Pierre <pierre....@gmail.com> wrote:
>>> I see. Well I *do* have hundreds of 2x2 matrices to multiply out so
>>> i'm better off storing them as numpy matrices throughout... thanks for
>>> your explanations though.
>>>
>>> Pierre
>>
>> You might consider using Cython and writing a custom 2x2 matrix class.
>> It wouldn't be difficult... so I'll write one right now and respond
>> with the benchmarks.
Here it is: http://480.sagenb.org/home/pub/97/
I ended up using GSL's complex matrix data type and the BLAS level 3
routine to do the multiplication. I did not add any other
convenience functions to the class, so some more will probably be
needed for your application.
The results:
sage: timeit('a*a2')
625 loops, best of 3: 643 ns per loop
sage: timeit('b*b2')
625 loops, best of 3: 25.8 µs per loop
sage: timeit('c.dot(c2)')
625 loops, best of 3: 2.52 µs per loop
sage: 2.52/.643
3.91912908242613
sage: 25.8/.643
40.1244167962675
So, my one off GSL based class is 4 times faster than numpy and 40
times faster than Sage at matrix squaring.
http://480.sagenb.org/home/pub/97/
Jason Grout
> On Mon, Feb 13, 2012 at 10:29 AM, Robert Bradshaw
> <robe...@math.washington.edu> wrote:
>> On Mon, Feb 13, 2012 at 10:06 AM, William Stein<wst...@gmail.com> wrote:
>>> On Mon, Feb 13, 2012 at 9:59 AM, Pierre<pierre....@gmail.com> wrote:
>>>> I see. Well I *do* have hundreds of 2x2 matrices to multiply out so
>>>> i'm better off storing them as numpy matrices throughout... thanks for
>>>> your explanations though.
>>>>
>>>> Pierre
>>>
>>> You might consider using Cython and writing a custom 2x2 matrix class.
>>> It wouldn't be difficult... so I'll write one right now and respond
>>> with the benchmarks.
>
> Here it is: http://480.sagenb.org/home/pub/97/
>
> I ended up using GSL's complex matrix data type and the BLAS level 3
> routine to do the multiplication. I did not add any other
> convenience functions to the class, so some more will probably be
> needed for your application.
I'm curious why you didn't just store the 4 complex numbers in C. I
tried it and got a much bigger speedup: 17x faster than numpy and 150x
faster than Sage. See http://sagenb.org/home/pub/4303/
Thanks,
Jason
William Stein
<jason...@creativetrax.com> wrote:
> On 2/13/12 12:41 PM, William Stein wrote:
>>
>> On Mon, Feb 13, 2012 at 10:29 AM, Robert Bradshaw
>> <robe...@math.washington.edu> wrote:
>>>
>>> On Mon, Feb 13, 2012 at 10:06 AM, William Stein<wst...@gmail.com> wrote:
>>>>
>>>> On Mon, Feb 13, 2012 at 9:59 AM, Pierre<pierre....@gmail.com>
>>>> wrote:
>>>>>
>>>>> I see. Well I *do* have hundreds of 2x2 matrices to multiply out so
>>>>> i'm better off storing them as numpy matrices throughout... thanks for
>>>>> your explanations though.
>>>>>
>>>>> Pierre
>>>>
>>>>
>>>> You might consider using Cython and writing a custom 2x2 matrix class.
>>>> It wouldn't be difficult... so I'll write one right now and respond
>>>> with the benchmarks.
>>
>>
>> Here it is: http://480.sagenb.org/home/pub/97/
>>
>> I ended up using GSL's complex matrix data type and the BLAS level 3
>> routine to do the multiplication. I did not add any other
>> convenience functions to the class, so some more will probably be
>> needed for your application.
>
>
> I'm curious why you didn't just store the 4 complex numbers in C.
I was concerned about "numerical stability" and figured BLAS would
deal with that. But maybe that doesn't matter.
Most importantly, I forgot about Cython's complex type, which I now
remember Robert Bradshaw wrote for his thesis work. I wasn't looking
forward to writing your line:
res.m00=self.m00*right.m00+self.m01*right.m10
but against the gsl library. But using the cython type in complex,
you get the above at C speed with easy notation, which is pretty
awesome.
> I tried
> it and got a much bigger speedup: 17x faster than numpy and 150x faster than
> Sage. See http://sagenb.org/home/pub/4303/
Nice.
>
> Thanks,
>
> Jason
Pierre
tend to restrict myself to basic C types and simple functions which do
not import anything -- it seems I can learn an awful lot by studying
your examples!
Can i ask the following questions about Jason's code?
** why so you 'import' CDF while you 'cimport' ComplexDoubleElement?
what is the difference?
** and why do these imports anyway, are they necessary for "cdef
complex..." to work?
** I didn't know PY_NEW. What is the difference between res=
PY_NEW(Matrix2) and simply res= Matrix2(0,0,0,0) (say) ? and what
about using PY_NEW with a class that *requires* parameters to its
__init__?
** more generally, is there a Cython tutorial that would be (just)
sufficiently advanced to cover such an example?
thanks!
Pierre
On 13 fév, 20:13, William Stein <wst...@gmail.com> wrote:
> On Mon, Feb 13, 2012 at 11:07 AM, Jason Grout
>
>
>
>
>
>
>
>
>
> <jason-s...@creativetrax.com> wrote:
> > On 2/13/12 12:41 PM, William Stein wrote:
>
> >> On Mon, Feb 13, 2012 at 10:29 AM, Robert Bradshaw
> >>>> On Mon, Feb 13, 2012 at 9:59 AM, Pierre<pierre.guil...@gmail.com>
Pierre
if i put all this in a foo.spyx file, i usually do "load foo.spyx". It
works, but compiles every time! how do you just load the compiled
version? and what about import it like a module, using
foo.function() ?
William Stein
> Hey, thanks for all the responses. I use Cython occasionally, but i
> tend to restrict myself to basic C types and simple functions which do
> not import anything -- it seems I can learn an awful lot by studying
> your examples!
>
> Can i ask the following questions about Jason's code?
>
> ** why so you 'import' CDF while you 'cimport' ComplexDoubleElement?
> what is the difference?
cimport gives access to the underlying C-level data structure.
> ** and why do these imports anyway, are they necessary for "cdef
> complex..." to work?
Nope! If you just want to use Python's complex data type you could
rewrite the code
to not use any imports.
> ** I didn't know PY_NEW. What is the difference between res=
> PY_NEW(Matrix2) and simply res= Matrix2(0,0,0,0) (say) ? and what
> about using PY_NEW with a class that *requires* parameters to its
> __init__?
PY_NEW is a macro that creates the class without calling __init__.
This is an optimization.
> ** more generally, is there a Cython tutorial that would be (just)
> sufficiently advanced to cover such an example?
You may find chapter 3 of this book I'm writing helpful:
http://wstein.org/books/sagebook/sagebook.pdf
William Stein
> oh and a related question -- been wondering about this for a while :
> if i put all this in a foo.spyx file, i usually do "load foo.spyx". It
> works, but compiles every time! how do you just load the compiled
> version?
> and what about import it like a module, using
> foo.function() ?
You can do all that -- it just costs you about 2 hours of background
education and practice. You have to use a .pyx file, probably a
setup.py file, and follow the standard instructions at cython.org:
http://docs.cython.org/src/quickstart/index.html
-- William
Pierre
> rewrite the code
> to not use any imports.
can be deleted, the class Matrix2 would still work as such?
> You may find chapter 3 of this book I'm writing helpful:
>
> http://wstein.org/books/sagebook/sagebook.pdf
you were writing, and I agreed to translate it into French (remember?
someone else was in charge of spanish i think). Is this the one? (I
guess not, since you didn't ask me)
(gmail tells me we had this conversation in august 2009 !!!)
Jason Grout
>> ** and why do these imports anyway, are they necessary for "cdef
>> > complex..." to work?
> Nope! If you just want to use Python's complex data type you could
> rewrite the code
> to not use any imports.
>
I had them because I was copying and modifying William's code and didn't
ever delete them.
Thanks,
Jason
Nils Bruin
> PY_NEW is a macro that creates the class without calling __init__.
> This is an optimization.
"Matrix2.__new__(Matrix2)" ?
See http://trac.cython.org/cython_trac/ticket/443.
I think I used it for that purpose in EclObject. It looks cleaner and
you can use it without an extra include. If you're writing this in a
book, you may want to check with the cython folk what the recommended
construct is. Plus I'd be interested to know if there is still a
reason to prefer PY_NEW.