[Numpy-discussion] problem with vectorized difference equation
eat
e.antero.tammi at gmail.com
Sat Apr 7 16:50:48 EDT 2012
Hi,
On Sat, Apr 7, 2012 at 1:27 AM, Sameer Grover <sameer.grover.1 at gmail.com>wrote:
> On Saturday 07 April 2012 02:51 AM, Francesco Barale wrote:
> > Hello Sameer,
> >
> > Thank you very much for your reply. My goal was to try to speed up the
> loop
> > describing the accumulator. In the (excellent) book I was mentioning in
> my
> > initial post I could find one example that seemed to match what I was
> trying
> > to do. Basically, it is said that a loop of the following kind:
> >
> > n = size(u)-1
> > for i in xrange(1,n,1):
> > u_new[i] = u[i-1] + u[i] + u[i+1]
> >
> > should be equivalent to:
> >
> > u[1:n] = u[0:n-1] + u[1:n] + u[i+1]
> This example is correct.
>
> What I was trying to point out was that the single expression y[1:n] =
> y[0:n-1] + u[1:n] will iterate over the array but will not accumulate.
> It will add y[0:n-1] to u[1:n] and assign the result to y[1:n].
>
> For example,
> y = [1,2,3,4]
> u = [5,6,7,8]
> Then y[0:n-1] = [1,2,3] and u[1:n]=[6,7,8]
>
> The statement y[1:n] = y[0:n-1] + u[1:n] implies
> y[1:n] = [6+1,7+2,8+3] = [7,9,11]
> yielding y = [1,7,9,11]
>
FWIIFO, if assignment in loop like this ever makes any sense (which I
doubt),
>
> whereas the code:
>
> for i in range(0,n-1):
> y[i+1] = y[i] + u[i+1]
>
then it can be captured in a function, like
In []: def s0(y, u):
..: yn= y.copy()
..: for k in xrange(y.size- 1):
..: yn[k+ 1]= yn[k]+ u[k+ 1]
..: return yn
..:
and now this function can be easily transformed to utilize cumsum, like
In []: def s1(y, u):
..: un= u.copy()
..: un[0]= y[0]
..: return cumsum(un)
..:
thus
In []: y, u= rand(1e5), rand(1e5)
In []: allclose(s0(y, u), s1(y, u))
Out[]: True
and definitely this transformation will outperform a plain python loop
In []: timeit s0(y, u)
10 loops, best of 3: 122 ms per loop
In []: timeit s1(y, u)
100 loops, best of 3: 2.16 ms per loop
In []: 122/ 2.16
Out[]: 56.48148148148148
My 2 cents,
-eat
>
> will accumulate and give y = [1,7,14,22]
>
> Sameer
> > Am I missing something?
> >
> > Regards,
> > Francesco
> >
> >
> > Sameer Grover wrote:
> >> On Saturday 07 April 2012 12:14 AM, francesco82 wrote:
> >>> Hello everyone,
> >>>
> >>> After reading the very good post
> >>>
> http://technicaldiscovery.blogspot.com/2011/06/speeding-up-python-numpy-cython-and.html
> >>> and the book by H. P. Langtangen 'Python scripting for computational
> >>> science' I was trying to speed up the execution of a loop on numpy
> arrays
> >>> being used to describe a simple difference equation.
> >>>
> >>> The actual code I am working on involves some more complicated
> equations,
> >>> but I am having the same exact behavior as described below. To test the
> >>> improvement in speed I wrote the following in vect_try.py:
> >>>
> >>> #!/usr/bin/python
> >>> import numpy as np
> >>> import matplotlib.pyplot as plt
> >>>
> >>> dt = 0.02 #time step
> >>> time = np.arange(0,2,dt) #time array
> >>> u = np.sin(2*np.pi*time) #input signal array
> >>>
> >>> def vect_int(u,y): #vectorized function
> >>> n = u.size
> >>> y[1:n] = y[0:n-1] + u[1:n]
> >>> return y
> >>>
> >>> def sc_int(u,y): #scalar function
> >>> y = y + u
> >>> return y
> >>>
> >>> def calc_vect(u, func=vect_int):
> >>> out = np.zeros(u.size)
> >>> for i in xrange(u.size):
> >>> out = func(u,out)
> >>> return out
> >>>
> >>> def calc_sc(u, func=sc_int):
> >>> out = np.zeros(u.size)
> >>> for i in xrange(u.size-1):
> >>> out[i+1] = sc_int(u[i+1],out[i])
> >>> return out
> >>>
> >>> To verify the execution time I've used the timeit function in Ipython:
> >>>
> >>> import vect_try as vt
> >>> timeit vt.calc_vect(vt.u) --> 1000 loops, best of 3: 494 us per loop
> >>> timeit vt.calc_sc(vt.u) -->10000 loops, best of 3: 92.8 us per loop
> >>>
> >>> As you can see the scalar implementation looping one item at the time
> >>> (calc_sc) is 494/92.8~5.3 times faster than the vectorized one
> >>> (calc_vect).
> >>>
> >>> My problem consists in the fact that I need to iterate the execution of
> >>> calc_vect in order for it to operate on all the elements of the input
> >>> array.
> >>> If I execute calc_vect only once, it will only operate on the first
> slice
> >>> of
> >>> the vectors leaving the remaining untouched. My understanding was that
> >>> the
> >>> vector expression y[1:n] = y[0:n-1] + u[1:n] was supposed to iterate
> over
> >>> all the array, but that's not happening for me. Can anyone tell me
> what I
> >>> am
> >>> doing wrong?
> >>>
> >>> Thanks!
> >>> Francesco
> >>>
> >> 1. By vectorizing, we mean replacing a loop with a single expression. In
> >> your program, both the scalar and vector implementations (calc_vect and
> >> calc_sc) have a loop each. This isn't going to make anything faster. The
> >> vectorized implementation is just a convoluted way of achieving the same
> >> result and is slower.
> >>
> >> 2. The expression y[1:n] = y[0:n-1] + u[1:n] is /not/ equivalent to the
> >> following loop
> >>
> >> for i in range(0,n-1):
> >> y[i+1] = y[i] + u[i+1]
> >>
> >> It is equivalent to something like
> >>
> >> z = np.zeros(n-1)
> >> for i in range(0,n-1):
> >> z[i] = y[i] + u[i+1]
> >> y[1:n] = z
> >>
> >> i.e., the RHS is computed in totality and then assigned to the LHS. This
> >> is how array operations work even in other languages such as Matlab.
> >>
> >> 3. I personally don't think there is a simple/obvious way to vectorize
> >> what you're trying to achieve.
> >>
> >> Sameer
> >>
> >> _______________________________________________
> >> NumPy-Discussion mailing list
> >> NumPy-Discussion at scipy.org
> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >>
> >>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20120407/f97015b1/attachment.html>
More information about the NumPy-Discussion
mailing list