[sage-support] Re: speed question, numpy vs CDF
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? ** 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 rober...@math.washington.edu wrote: On Mon, Feb 13, 2012 at 10:06 AM, William Steinwst...@gmail.com wrote: On Mon, Feb 13, 2012 at 9:59 AM, Pierrepierre.guil...@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. Seehttp://sagenb.org/home/pub/4303/ Nice. Thanks, Jason -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 Washingtonhttp://wstein.org -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
[sage-support] Re: speed question, numpy vs CDF
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() ? On 14 fév, 18:11, Pierre pierre.guil...@gmail.com wrote: 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? ** 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 rober...@math.washington.edu wrote: On Mon, Feb 13, 2012 at 10:06 AM, William Steinwst...@gmail.com wrote: On Mon, Feb 13, 2012 at 9:59 AM, Pierrepierre.guil...@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. Seehttp://sagenb.org/home/pub/4303/ Nice. Thanks, Jason -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 Washingtonhttp://wstein.org -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
On Tue, Feb 14, 2012 at 9:11 AM, Pierre pierre.guil...@gmail.com wrote: 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 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 rober...@math.washington.edu wrote: On Mon, Feb 13, 2012 at 10:06 AM, William Steinwst...@gmail.com wrote: On Mon, Feb 13, 2012 at 9:59 AM, Pierrepierre.guil...@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. Seehttp://sagenb.org/home/pub/4303/ Nice. Thanks, Jason -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 Washingtonhttp://wstein.org -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
On Tue, Feb 14, 2012 at 9:14 AM, Pierre pierre.guil...@gmail.com wrote: 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 On 14 fév, 18:11, Pierre pierre.guil...@gmail.com wrote: 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? ** 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 rober...@math.washington.edu wrote: On Mon, Feb 13, 2012 at 10:06 AM, William Steinwst...@gmail.com wrote: On Mon, Feb 13, 2012 at 9:59 AM, Pierrepierre.guil...@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. Seehttp://sagenb.org/home/pub/4303/ Nice. Thanks, Jason -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 Washingtonhttp://wstein.org -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
[sage-support] Re: speed question, numpy vs CDF
Nope! If you just want to use Python's complex data type you could rewrite the code to not use any imports. just to be clear: you mean that the first two lines of Jason's code 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 this reminds me that ages ago, you mentioned a book about Sage that 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 !!!) -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
[sage-support] Re: speed question, numpy vs CDF
On 2/14/12 11:20 AM, William Stein wrote: ** 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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
[sage-support] Re: speed question, numpy vs CDF
On Feb 14, 9:20 am, William Stein wst...@gmail.com wrote: PY_NEW is a macro that creates the class without calling __init__. This is an optimization. Hasn't that construction been made obsolete by 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. -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
[sage-support] Re: speed question, numpy vs CDF
Hitting the nail one the head! sounds like your solution (I mean numpy arrays with CDF elements, which I didn't know was possible) is going to be perfect for me. Thank you so much! On 13 fév, 01:30, Nils Bruin nbr...@sfu.ca wrote: On Feb 12, 1:39 pm, Pierre pierre.guil...@gmail.com wrote: 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. No, that's not the problem. It's simply that numpy's default complex 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 625 loops, best of 3: 11.8 µs per loop sage: import numpy sage: AA= numpy.array(A); BB= numpy.array(B) sage: %timeit AA.dot(BB) 625 loops, best of 3: 1.28 µs per loop 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' -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
On Sun, Feb 12, 2012 at 4:30 PM, Nils Bruin nbr...@sfu.ca wrote: On Feb 12, 1:39 pm, Pierre pierre.guil...@gmail.com wrote: 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. No, that's not the problem. It's simply that numpy's default complex 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 625 loops, best of 3: 11.8 µs per loop sage: import numpy sage: AA= numpy.array(A); BB= numpy.array(B) sage: %timeit AA.dot(BB) 625 loops, best of 3: 1.28 µs per loop 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' 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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
[sage-support] Re: speed question, numpy vs CDF
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 On 13 fév, 18:32, Robert Bradshaw rober...@math.washington.edu wrote: On Sun, Feb 12, 2012 at 4:30 PM, Nils Bruin nbr...@sfu.ca wrote: On Feb 12, 1:39 pm, Pierre pierre.guil...@gmail.com wrote: 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. No, that's not the problem. It's simply that numpy's default complex 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 625 loops, best of 3: 11.8 µs per loop sage: import numpy sage: AA= numpy.array(A); BB= numpy.array(B) sage: %timeit AA.dot(BB) 625 loops, best of 3: 1.28 µs per loop 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' 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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
On Mon, Feb 13, 2012 at 9:59 AM, Pierre pierre.guil...@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. -- William On 13 fév, 18:32, Robert Bradshaw rober...@math.washington.edu wrote: On Sun, Feb 12, 2012 at 4:30 PM, Nils Bruin nbr...@sfu.ca wrote: On Feb 12, 1:39 pm, Pierre pierre.guil...@gmail.com wrote: 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. No, that's not the problem. It's simply that numpy's default complex 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 625 loops, best of 3: 11.8 µs per loop sage: import numpy sage: AA= numpy.array(A); BB= numpy.array(B) sage: %timeit AA.dot(BB) 625 loops, best of 3: 1.28 µs per loop 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' 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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
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.guil...@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. On 13 fév, 18:32, Robert Bradshaw rober...@math.washington.edu wrote: On Sun, Feb 12, 2012 at 4:30 PM, Nils Bruin nbr...@sfu.ca wrote: On Feb 12, 1:39 pm, Pierre pierre.guil...@gmail.com wrote: 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. No, that's not the problem. It's simply that numpy's default complex 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 625 loops, best of 3: 11.8 µs per loop sage: import numpy sage: AA= numpy.array(A); BB= numpy.array(B) sage: %timeit AA.dot(BB) 625 loops, best of 3: 1.28 µs per loop 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' 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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
Anyone using numpy from Sage should beware of the following: 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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
Re: [sage-support] Re: speed question, numpy vs CDF
On Mon, Feb 13, 2012 at 10:29 AM, Robert Bradshaw
rober...@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.guil...@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/
+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.
On 13 fév, 18:32, Robert Bradshaw rober...@math.washington.edu
wrote:
On Sun, Feb 12, 2012 at 4:30 PM, Nils Bruin nbr...@sfu.ca wrote:
On Feb 12, 1:39 pm, Pierre pierre.guil...@gmail.com wrote:
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.
No, that's not the problem. It's simply that numpy's default complex
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
625 loops, best of 3: 11.8 µs per loop
sage: import numpy
sage: AA= numpy.array(A); BB= numpy.array(B)
sage: %timeit AA.dot(BB)
625 loops, best of 3: 1.28 µs per loop
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'
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
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support+unsubscr...@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
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support+unsubscr...@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
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
Re: [sage-support] Re: speed question, numpy vs CDF
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 rober...@math.washington.edu wrote: On Mon, Feb 13, 2012 at 10:06 AM, William Steinwst...@gmail.com wrote: On Mon, Feb 13, 2012 at 9:59 AM, Pierrepierre.guil...@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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@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 -- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
