5
\$\begingroup\$

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?

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

3 Answers 3

Reset to default
4
\$\begingroup\$

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

\$\endgroup\$
1
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Commented Dec 7, 2017 at 15:14
4
\$\begingroup\$

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

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

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

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.