Trying to vectorize iterative calculation with numpy

Question

I am trying to make part of the code more efficient using vectorized form in numpy. Let me show you an example so that you know what I mean.

Given the following code:

a = np.zeros([4,4]) a[0] = [1., 2., 3., 4.] for i in range(len(a)-1): a[i+1] = 2*a[i] print a 

He outputs

 [[ 1. 2. 3. 4.] [ 2. 4. 6. 8.] [ 4. 8. 12. 16.] [ 8. 16. 24. 32.]] 

When I now try to vectorize the code as follows:

 a = np.zeros([4,4]) a[0] = [1., 2., 3., 4.] a[1:] = 2*a[0:-1] print a 

I get only the first iteration:

 [[ 1. 2. 3. 4.] [ 2. 4. 6. 8.] [ 0. 0. 0. 0.] [ 0. 0. 0. 0.]] 

Is it possible to write the code above efficiently in vectorized form (where the next iteration always refers to the previous iteration) or do I need to maintain a for loop?

Answer 1

A linear recursion such as this can be calculated using scipy.signal.lfilter :

 In [19]: from scipy.signal import lfilter In [20]: num = np.array([1.0]) In [21]: alpha = 2.0 In [22]: den = np.array([1.0, -alpha]) In [23]: a = np.zeros((4,4)) In [24]: a[0,:] = [1,2,3,4] In [25]: lfilter(num, den, a, axis=0) Out[25]: array([[ 1., 2., 3., 4.], [ 2., 4., 6., 8.], [ 4., 8., 12., 16.], [ 8., 16., 24., 32.]]) 

For more information, see the following: recursive python vectorization with timings , Recursive definitions in Pandas

---

Please note that using lfilter really makes sense if you are solving a heterogeneous problem, such as x[i+1] = alpha*x[i] + u[i] , where u is the given input array. For a simple repetition of a[i+1] = alpha*a[i] you can use the exact solution a[i] = a[0]*alpha**i . A solution for multiple seed values ​​can be vectorized using broadcast. For instance,

 In [271]: alpha = 2.0 In [272]: a0 = np.array([1, 2, 3, 4]) In [273]: n = 5 In [274]: a0 * (alpha**np.arange(n).reshape(-1, 1)) Out[274]: array([[ 1., 2., 3., 4.], [ 2., 4., 6., 8.], [ 4., 8., 12., 16.], [ 8., 16., 24., 32.], [ 16., 32., 48., 64.]]) 

Answer 2

Calculation of bill vectors acts on a vector, and not as a sequence of steps, so you need to vectorize the entire expression. For instance:

 np.multiply(np.arange(1,5), 2**np.arange(0,4)[np.newaxis].T) 

To solve the "final" question, yes, you have to keep the for loop if you want to do a consistent calculation. You can make it more efficient with map or [... for ...] , but it takes a lot of trial and error to optimize this path. The beauty of thinking in vector terms and using Numpy for implementation is that you get the result efficiently without all trial and error.

The cumsum and cumprod can do something similar to what you are asking for. Instead of 2**np.arange(...) you can get the same from

 np.multiply(np.arange(1,5), np.cumprod([1,2,2,2,])[np.newaxis].T) 

Answer 3

view your result matrix

 result = [[ 1, 2, 3, 4,] [ 2, 4, 6, 8,] [ 4, 8, 12, 16,] [ 8, 16, 24, 32,]] 

it can be deconstructed into a product (by type) of two matrices as

 a = [[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]] b = [[1, 1, 1, 1], [2, 2, 2, 2], [4, 4, 4, 4] [8, 8, 8, 8]] result = a * b 

you can calculate this type of operation using the meshgrid function

 aa, bb = np.meshgrid(np.array([1.0, 2.0, 3.0, 4.0]), np.array([1.0, 2.0, 4.0, 8.0])) result = aa * bb