<br><br><div class="gmail_quote">On Sat, Nov 12, 2011 at 11:16 AM,  <span dir="ltr"><<a href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">

On Sat, Nov 12, 2011 at 11:32 AM, Warren Weckesser<br>
<div><div></div><div class="h5"><<a href="mailto:warren.weckesser@enthought.com">warren.weckesser@enthought.com</a>> wrote:<br>
><br>
><br>
> On Sat, Nov 12, 2011 at 9:59 AM, <<a href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>> wrote:<br>
>><br>
>> On Sat, Nov 12, 2011 at 10:31 AM, Warren Weckesser<br>
>> <<a href="mailto:warren.weckesser@enthought.com">warren.weckesser@enthought.com</a>> wrote:<br>
>> ><br>
>> ><br>
>> > On Sat, Nov 12, 2011 at 6:43 AM, <<a href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>> wrote:<br>
>> >><br>
>> >> On Sat, Nov 12, 2011 at 3:36 AM, Geoffrey Zhu <<a href="mailto:zyzhu2000@gmail.com">zyzhu2000@gmail.com</a>><br>
>> >> wrote:<br>
>> >> > Hi,<br>
>> >> ><br>
>> >> > I am playing with multiple ways to speed up the following expression<br>
>> >> > (it is in the inner loop):<br>
>> >> ><br>
>> >> ><br>
>> >> > C[1:(M - 1)]=(a * C[2:] + b * C[1:(M-1)] + c * C[:(M-2)])<br>
>> >> ><br>
>> >> > where C is an array of about 200-300 elements, M=len(C), a, b, c are<br>
>> >> > scalars.<br>
>> >> ><br>
>> >> > I played with numexpr, but it was way slower than directly using<br>
>> >> > numpy. It would be nice if I could create a Mx3 matrix without<br>
>> >> > copying<br>
>> >> > memory and so I can use dot() to calculate the whole thing.<br>
>> >> ><br>
>> >> > Can anyone help with giving some advices to make this faster?<br>
>> >><br>
>> >> looks like a np.convolve(C, [a,b,c])  to me except for the boundary<br>
>> >> conditions.<br>
>> ><br>
>> ><br>
>> ><br>
>> > As Josef pointed out, this is a convolution.  There are (at least)<br>
>> > three convolution functions in numpy+scipy that you could use:<br>
>> > numpy.convolve, scipy.signal.convolve, and scipy.ndimage.convolve1d.<br>
>> > Of these, scipy.ndimage.convolve1d is the fastest.  However, it doesn't<br>
>> > beat the simple expression<br>
>> >    a*x[2:] + b*x[1:-1] + c*x[:-2]<br>
>> > Your idea of forming a matrix without copying memory can be done using<br>
>> > "stride tricks", and for arrays of the size you are interested in, it<br>
>> > computes the result faster than the simple expression (see below).<br>
>> ><br>
>> > Another fast alternative is to use one of the inline code generators.<br>
>> > This example is a easy to implement with scipy.weave.blitz, and it gives<br>
>> > a big speedup.<br>
>> ><br>
>> > Here's a test script:<br>
>> ><br>
>> > #----- convolve1dtest.py -----<br>
>> ><br>
>> ><br>
>> > import numpy as np<br>
>> > from numpy.lib.stride_tricks import as_strided<br>
>> > from scipy.ndimage import convolve1d<br>
>> > from scipy.weave import blitz<br>
>> ><br>
>> > # weighting coefficients<br>
>> > a = -0.5<br>
>> > b = 1.0<br>
>> > c = -0.25<br>
>> > w = np.array((a,b,c))<br>
>> > # Reversed w:<br>
>> > rw = w[::-1]<br>
>> ><br>
>> > # Length of C<br>
>> > n = 250<br>
>> ><br>
>> > # The original version of the calculation:<br>
>> > # Some dummy data<br>
>> > C = np.arange(float(n))<br>
>> > C[1:-1] = a * C[2:] + b * C[1:-1] + c * C[:-2]<br>
>> > # Save for comparison.<br>
>> > C0 = C.copy()<br>
>> ><br>
>> > # Do it again using a matrix multiplication.<br>
>> > C = np.arange(float(n))<br>
>> > # The "virtual" matrix view of C.<br>
>> > V = as_strided(C, shape=(C.size-2, 3), strides=(C.strides[0],<br>
>> > C.strides[0]))<br>
>> > C[1:-1] = np.dot(V, rw)<br>
>> > C1 = C.copy()<br>
>> ><br>
>> > # Again, with convolve1d this time.<br>
>> > C = np.arange(float(n))<br>
>> > C[1:-1] = convolve1d(C, w)[1:-1]<br>
>> > C2 = C.copy()<br>
>> ><br>
>> > # scipy.weave.blitz<br>
>> > C = np.arange(float(n))<br>
>> > # Must work with a copy, D, in the formula, because blitz does not use<br>
>> > # a temporary variable.<br>
>> > D = C.copy()<br>
>> > expr = "C[1:-1] = a * D[2:] + b * D[1:-1] + c * D[:-2]"<br>
>> > blitz(expr, check_size=0)<br>
>> > C3 = C.copy()<br>
>> ><br>
>> ><br>
>> > # Verify that all the methods give the same result.<br>
>> > print np.all(C0 == C1)<br>
>> > print np.all(C0 == C2)<br>
>> > print np.all(C0 == C3)<br>
>> ><br>
>> > #-----<br>
>> ><br>
>> > And here's a snippet from an ipython session:<br>
>> ><br>
>> > In [51]: run convolve1dtest.py<br>
>> > True<br>
>> > True<br>
>> > True<br>
>> ><br>
>> > In [52]: %timeit C[1:-1] = a * C[2:] + b * C[1:-1] + c * C[:-2]<br>
>> > 100000 loops, best of 3: 16.5 us per loop<br>
>> ><br>
>> > In [53]: %timeit C[1:-1] = np.dot(V, rw)<br>
>> > 100000 loops, best of 3: 9.84 us per loop<br>
>> ><br>
>> > In [54]: %timeit C[1:-1] = convolve1d(C, w)[1:-1]<br>
>> > 100000 loops, best of 3: 18.7 us per loop<br>
>> ><br>
>> > In [55]: %timeit D = C.copy(); blitz(expr, check_size=0)<br>
>> > 100000 loops, best of 3: 4.91 us per loop<br>
>> ><br>
>> ><br>
>> ><br>
>> > scipy.weave.blitz is fastest (but note that blitz has already been<br>
>> > called<br>
>> > once, so the time shown does not include the compilation required in<br>
>> > the first call).  You could also try scipy.weave.inline, cython.inline,<br>
>> > or np_inline (<a href="http://pages.cs.wisc.edu/%7Ejohnl/np_inline/" target="_blank">http://pages.cs.wisc.edu/~johnl/np_inline/</a>).<br>
>> ><br>
>> > Also note that C[-1:1] = np.dot(V, rw) is faster than either the simple<br>
>> > expression or convolve1d.  However, if you also have to set up V inside<br>
>> > your inner loop, the speed gain will be lost.  The relative speeds also<br>
>> > depend on the size of C.  For large C, the simple expression is faster<br>
>> > than the matrix multiplication by V (but blitz is still faster).  In<br>
>> > the following, I have changed n to 2500 before running<br>
>> > convolve1dtest.py:<br>
>> ><br>
>> > In [56]: run convolve1dtest.py<br>
>> > True<br>
>> > True<br>
>> > True<br>
>> ><br>
>> > In [57]: %timeit C[1:-1] = a * C[2:] + b * C[1:-1] + c * C[:-2]<br>
>> > 10000 loops, best of 3: 29.5 us per loop<br>
>> ><br>
>> > In [58]: %timeit C[1:-1] = np.dot(V, rw)<br>
>> > 10000 loops, best of 3: 56.4 us per loop<br>
>> ><br>
>> > In [59]: %timeit C[1:-1] = convolve1d(C, w)[1:-1]<br>
>> > 10000 loops, best of 3: 37.3 us per loop<br>
>> ><br>
>> > In [60]: %timeit D = C.copy(); blitz(expr, check_size=0)<br>
>> > 100000 loops, best of 3: 10.3 us per loop<br>
>> ><br>
>> ><br>
>> > blitz wins, the simple numpy expression is a distant second, and now<br>
>> > the matrix multiplication is slowest.<br>
>> ><br>
>> > I hope that helps--I know I learned quite a bit. :)<br>
>><br>
>> Interesting, two questions<br>
>><br>
>> does scipy.signal convolve have a similar overhead as np.convolve1d ?<br>
><br>
><br>
> Did you mean np.convolve?  There is no np.convolve1d.   Some of the tests<br>
> that I've done with convolution functions are here:<br>
>     <a href="http://www.scipy.org/Cookbook/ApplyFIRFilter" target="_blank">http://www.scipy.org/Cookbook/ApplyFIRFilter</a><br>
> I should add np.convolve to that page.  For the case considered here,<br>
> np.convolve is a bit slower than scipy.ndimage.convolve1d, but for larger<br>
> arrays, it looks like np.convolve can be much faster.<br>
><br>
><br>
>><br>
>> memory:<br>
>> the blitz code doesn't include the array copy (D), so the timing might<br>
>> be a bit misleading?<br>
><br>
><br>
> Look again at my %timeit calls in the ipython snippets.  :)<br>
><br>
<br>
</div></div>Sorry, I was reading way too fast or selectively (skipping the<br>
imports, namespaces for example).<br>
(I thought convolve1d was numpy not ndimage)<br>
<br>
I tried out your cookbook script a while ago, it's nice, I had only a<br>
rough idea about the timing before.<br>
Compared to the current case the x is much longer, (unless I don't<br>
remember or read it correctly.)<br></blockquote><div><br><br>I just updated the page with np.convolve:<br>    <a href="http://www.scipy.org/Cookbook/ApplyFIRFilter">http://www.scipy.org/Cookbook/ApplyFIRFilter</a><br>A big part of the overhead in using np.convolve in that script is that the function only accepts 1-d arrays, so a python loop is required to apply it along an axis of a 2-d array.<br>

<br>Warren<br><br></div></div>