[Numpy-discussion] problem with vectorized difference equation

Francesco Barale francesco.barale@gmail....
Fri Apr 6 16:21:48 CDT 2012


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]

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@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> 

-- 
View this message in context: http://old.nabble.com/problem-with-vectorized-difference-equation-tp33641688p33645699.html
Sent from the Numpy-discussion mailing list archive at Nabble.com.



More information about the NumPy-Discussion mailing list