[Numpy-discussion] problem with vectorized difference equation

[Sameer Grover]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20120407/6f2a9082/attachment.html>