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

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

\$\endgroup\$
  • 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\$ – Delgan Dec 7 '17 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\$
  • \$\begingroup\$ Indeed, I thought that there might be an optimized way to do that in Numpy, but it seems not. \$\endgroup\$ – Delgan Dec 7 '17 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\$ – peufeu Dec 7 '17 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\$
  • \$\begingroup\$ I am not sure about your updated populate() function, it does not add the calculated values to s. \$\endgroup\$ – Delgan Dec 7 '17 at 17:25
  • \$\begingroup\$ I thought it did, have you run it? \$\endgroup\$ – 13ros27 Dec 7 '17 at 17:26
  • 1
    \$\begingroup\$ Yes, it returns an array full of 0. Look at the code: it never modifies s. \$\endgroup\$ – Delgan Dec 7 '17 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\$ – Delgan Dec 7 '17 at 17:31
  • 2
    \$\begingroup\$ Yes, I think you should edit, because currently it does not seem to work as the initial populate(). \$\endgroup\$ – Delgan Dec 7 '17 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.