Vectorizing a sequence \$x_{i+1} = f(x_i)\$ with NumPy

Question

I would like to populate a NumPy array with "accumulated" results from a custom function.

Currently, my code is:

import numpy as np
def f(x, mu):
 return mu * x * (1 - x)
def populate(x0, mu, n):
 s = np.zeros(n)
 x = x0
 for i in range(n):
 s[i], x = x, f(x, mu)
 return s

It does not take advantage of the vectorization performance of NumPy.

Is there any way to improve the speed of creating arrays like this?

• Delgan: @Graipher Are you sere about that? I tested it and s[0] == x0.
• Graipher: You are right, I misread the swap.
• Snowbody: You might be able to gain a tiny bit of speed by calculating the composition of f with itself say 4 times, then call that on four values at the same time. It wastes a bunch of work but might end up being faster. f(f(x)) = r(rx(1-x))(1-rx(1-x)), do the algebra, come up with a polynomial. Sort of like unrolling the loop.
• Acccumulation: @Snowbody that's an interesting idea. You might get different answers based on rounding errors, thought.
• Snowbody: I did the algebra. There's 69 terms (constant * power of mu * power of x) that have to be added up for the 4x-composition. There's no way this is faster, even if you get the full speedup. I conclude this is not parallelizable.

Answer 1

You can't vectorize this. Repeated application of \$x_i = \mu x_{i-1} (1-x_{i-1})\$ is chaotic.

• Delgan: Yes, I realize that trying to "vectorize" this kind of function makes no sense. But I wonder if there is a way to avoid the Python loop and use only Numpy features.

Answer 2

As the others said it is not vectorizable via numpy since each iteration depends on the previous application of the function.

However, if you need speed, you can always use numba to JIT-compile the whole python loop into machine code...

Time to run populate(0.2,0.3,1000000), for three iterations:

Without numba

0.749
0.766
0.727

With numba

0.356
0.010
0.007

The first iteration is not much faster than original python, as it includes the compilation to machine code (I presume some libraries have to be loaded, etc). The others are a lot faster... 10ms versus 750ms in pure python...

• Delgan: Indeed, I thought that there might be an optimized way to do that in Numpy, but it seems not.
• bobflux: I made a little benchmark, considering adding numba to your program takes 2 lines of code, a 75x speedup on this function aint too bad!

Answer 3

It is not a bad thing to have your code take up more space!

s[i] = x
x = f(x, mu)

is a lot easier to read than

s[i], x = x, f(x, mu)

So overall that would be

import numpy as np 
def f(x, mu):
 return mu * x * (1 - x)
def populate(x0, mu, n):
 s = np.zeros(n)
 x = x0
 for i in range(n):
 s[i] = x
 x = f(x, mu)
 return s

• Delgan: I am not sure about your updated populate() function, it does not add the calculated values to s.
• 13ros27: I thought it did, have you run it?
• Delgan: Yes, it returns an array full of 0. Look at the code: it never modifies s.
• Delgan: But num is just a float, so "modifying" it is actually just assigning a new value to the variable num.
• Delgan: Yes, I think you should edit, because currently it does not seem to work as the initial populate().