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?

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

Improve this answer
\$\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\$
    – Delgan
    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...

Improve this answer
\$\endgroup\$
2
  • \$\begingroup\$ Indeed, I thought that there might be an optimized way to do that in Numpy, but it seems not. \$\endgroup\$
    – Delgan
    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\$
    – bobflux
    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
Improve this answer
\$\endgroup\$
7
  • \$\begingroup\$ I am not sure about your updated populate() function, it does not add the calculated values to s. \$\endgroup\$
    – Delgan
    Dec 7, 2017 at 17:25
  • \$\begingroup\$ I thought it did, have you run it? \$\endgroup\$
    – 13ros27
    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\$
    – Delgan
    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\$
    – Delgan
    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\$
    – Delgan
    Dec 7, 2017 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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