<br><br><div class="gmail_quote">On Sat, Nov 12, 2011 at 6:43 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;">
<div class="im">On Sat, Nov 12, 2011 at 3:36 AM, Geoffrey Zhu <<a href="mailto:zyzhu2000@gmail.com">zyzhu2000@gmail.com</a>> 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 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 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>
</div>looks like a np.convolve(C, [a,b,c]) to me except for the boundary conditions.<br></blockquote><div><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><span style="font-family: courier new,monospace;">#----- convolve1dtest.py -----</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">import numpy as np</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">from numpy.lib.stride_tricks import as_strided</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">from scipy.ndimage import convolve1d</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">from scipy.weave import blitz</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># weighting coefficients</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">a = -0.5</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">b = 1.0</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">c = -0.25</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">w = np.array((a,b,c))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Reversed w:</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">rw = w[::-1]</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Length of C</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">n = 250</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># The original version of the calculation:</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"># Some dummy data</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">C = np.arange(float(n))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C[1:-1] = a * C[2:] + b * C[1:-1] + c * C[:-2]</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Save for comparison.</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C0 = C.copy()</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Do it again using a matrix multiplication.</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C = np.arange(float(n))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># The "virtual" matrix view of C.</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">V = as_strided(C, shape=(C.size-2, 3), strides=(C.strides[0], C.strides[0]))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">C[1:-1] = np.dot(V, rw)</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C1 = C.copy()</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Again, with convolve1d this time.</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C = np.arange(float(n))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">C[1:-1] = convolve1d(C, w)[1:-1]</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C2 = C.copy()</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># scipy.weave.blitz</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C = np.arange(float(n))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Must work with a copy, D, in the formula, because blitz does not use</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"># a temporary variable.</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">D = C.copy()</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">expr = "C[1:-1] = a * D[2:] + b * D[1:-1] + c * D[:-2]"</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">blitz(expr, check_size=0)</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">C3 = C.copy()</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Verify that all the methods give the same result.</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">print np.all(C0 == C1)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">print np.all(C0 == C2)</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">print np.all(C0 == C3)</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">#-----</span><br>
<br>And here's a snippet from an ipython session:<br><br><span style="font-family: courier new,monospace;">In [51]: run convolve1dtest.py</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">True</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">True</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">True</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">In [52]: %timeit C[1:-1] = a * C[2:] + b * C[1:-1] + c * C[:-2]</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">100000 loops, best of 3: 16.5 us per loop</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [53]: %timeit C[1:-1] = np.dot(V, rw)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">100000 loops, best of 3: 9.84 us per loop</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [54]: %timeit C[1:-1] = convolve1d(C, w)[1:-1]</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">100000 loops, best of 3: 18.7 us per loop</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [55]: %timeit D = C.copy(); blitz(expr, check_size=0)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">100000 loops, best of 3: 4.91 us per loop</span><br style="font-family: courier new,monospace;">
<br><br><br>scipy.weave.blitz is fastest (but note that blitz has already been 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/~johnl/np_inline/">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 convolve1dtest.py:<br><br><span style="font-family: courier new,monospace;">In [56]: run convolve1dtest.py</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">True</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">True</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">True</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">In [57]: %timeit C[1:-1] = a * C[2:] + b * C[1:-1] + c * C[:-2]</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">10000 loops, best of 3: 29.5 us per loop</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [58]: %timeit C[1:-1] = np.dot(V, rw)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">10000 loops, best of 3: 56.4 us per loop</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [59]: %timeit C[1:-1] = convolve1d(C, w)[1:-1]</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">10000 loops, best of 3: 37.3 us per loop</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [60]: %timeit D = C.copy(); blitz(expr, check_size=0)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">100000 loops, best of 3: 10.3 us per loop</span><br style="font-family: courier new,monospace;">
<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><br>Warren<br> <br><br></div><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<br>
Josef<br>
<div><div></div><div class="h5"><br>
<br>
><br>
> Thanks,<br>
> G<br>
> _______________________________________________<br>
> NumPy-Discussion mailing list<br>
> <a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
> <a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
><br>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</div></div></blockquote></div><br>