[Numpy-discussion] problem with vectorized difference equation

[Sameer Grover]

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]

whereas the code:

for i in range(0,n-1):
      y[i+1] = y[i] + u[i+1]

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
>>
>>