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

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?

• @Graipher Are you sere about that? I tested it and s == x0. Dec 7, 2017 at 17:12
• You are right, I misread the swap. Dec 7, 2017 at 17:18
• 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. Dec 7, 2017 at 18:14
• @Snowbody that's an interesting idea. You might get different answers based on rounding errors, thought. Dec 7, 2017 at 20:28
• 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. Dec 11, 2017 at 21:12

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

• 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. Dec 7, 2017 at 15:14

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

• Indeed, I thought that there might be an optimized way to do that in Numpy, but it seems not. Dec 7, 2017 at 17:29
• 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! Dec 7, 2017 at 18:20

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

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