Building a numpy array as a function of previous element

Question

I would like to create a numpy array where the first element is a defined constant, and every next element is defined as the function of the previous element in the following way:

import numpy as np

def build_array_recursively(length, V_0, function):
    returnList = np.empty(length)
    returnList[0] = V_0
    for i in range(1,length):
        returnList[i] = function(returnList[i-1])
    return returnList

d_t = 0.05
print(build_array_recursively(20, 0.3, lambda x: x-x*d_t+x*x/2*d_t*d_t-x*x*x/6*d_t*d_t*d_t))

The print method above outputs

[0.3 0.28511194 0.27095747 0.25750095 0.24470843 0.23254756 0.22098752
 0.20999896 0.19955394 0.18962586 0.18018937 0.17122037 0.16269589
 0.15459409 0.14689418 0.13957638 0.13262186 0.1260127 0.11973187 0.11376316]

Is there a fast way of doing this in numpy without a for loop? If so is there a way to handle two elements before the current one, e.g. can a Fibonacci array be constructed similarly? I found a similar question here

Is it possible to vectorize recursive calculation of a NumPy array where each element depends on the previous one?

but was not answered in general. In my example, the difference equation is difficult to solve manually.

Answer 1

This is faster for what you want to do. You don't have to use recursion for the function.

Calculate the element based on previous element. Append calculated element to a list, and then change the list to numpy.

def method2(length, V_0, d_t):
    k = [V_0]
    x = V_0
    for i in range(1, length):
        x = x - x * d_t + x * x / 2 * d_t * d_t - x * x * x / 6 * d_t * d_t * d_t
        k.append(x)
    return np.asarray(k)

print(method2(20,0.3, 0.05))

Running you existing method 10000 times takes 0.438 seconds, while method2 takes 0.097 seconds.

Answer 2

Using a function to make the code clearer (instead of the inline lambda):

def fn(x):
    return x-x*d_t+x*x/2*d_t*d_t-x*x*x/6*d_t*d_t*d_t

And a function that combines elements of build_array_recursively and method2:

def foo1(length, V_0, function):
    returnList = np.empty(length)
    returnList[0] = x = V_0
    for i in range(1,length):
        returnList[i] = x = function(x)
    return returnList

In [887]: timeit build_array_recursively(20,0.3, fn);                                                                     
61.4 µs ± 63 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [888]: timeit method2(20,0.3, fn);                                                                  
16.9 µs ± 103 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [889]: timeit foo1(20,0.3, fn);                                                                     
13 µs ± 29.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

The main time saver in method2 and foo2 is carrying over x, the last value, from one iteration to the next, rather than indexing with returnList[i-1].

The accumulation method, assigning to a preallocated array, or list append, is less important. Performance is usually similar.

Here the calculation is simple enough that details of what you do in the loop makes a big difference in the overall time.

All of these are loops. Some ufunc have a reduce (and accumulate) method, that can apply the function repeatedly to a elements of the input array. np.sum, np.cumsum, etc make use of this. But you can't do that with a general Python function.

You have to use some sort of compilation tool like numba to perform this sort of loop much faster.