CODEWITHSUNDEEP
pythonarraysnumpyvectorization
PyRsquared
PyRsquared

Reputation: 7338

How can I apply a vectorized function to the previous element of a numpy array?

I want to apply a function like this:

s[i] = a*x[i] + (1 - a)*s[i-1]

where s and x are both arrays of the same length.

I don't want to use a for loop as these arrays are very large (>50 mil). I have tried doing something like this

def f(a,x):
    s = [0]*len(x)
    s[i] = a*x[i] + (1 - a)*s[i-1]
    return s

but of course i isn't defined so this doesn't work.

Is there a way to do this using map or numpy.apply_along_axis or some other vectorized method?

I haven't come across a method that applies functions to current and previous elements of an array without using for loops, and that is really what I want to understand how to do here.

EDIT

To be unambiguous, here is the for loop implementation which works but that I want to avoid

s = [0]*len(x)
a=0.45
for i in range(len(x)):
    s[i] = a*x[i] + (1-a)*s[i-1]

s[0] = x[0] # reset value of s[0]

Upvotes: 2

Views: 1485

Answers (2)

javidcf
javidcf

Reputation: 59711

You can avoid the loop, although rather than vectorizing is more like "computing everything for each value".

import numpy as np

# Loop
def fun(x, a):
    s = [0] * len(x)
    for i in range(len(x)):
        s[i] = a * x[i] + (1 - a) * s[i - 1]
    s[0] = x[0]
    return s

# Vectorized
def fun_vec(x, a):
    x = np.asarray(x, dtype=np.float32)
    n = np.arange(len(x))
    p = a * (1 - a) ** n
    # Trick from here: https://stackoverflow.com/q/49532575/1782792
    pz = np.concatenate((np.zeros(len(p) - 1, dtype=p.dtype), p))
    pp = np.lib.stride_tricks.as_strided(
        pz[len(p) - 1:], (len(p), len(p)),
        (p.strides[0], -p.strides[0]), writeable=False)
    t = x[np.newaxis] * pp
    s = np.sum(t, axis=1)
    s[0] = x[0]
    return s


x = list(range(1, 11))
a = 0.45
print(np.allclose(fun(x, a), fun_vec(x, a)))
# True

This kind of strategy is takes O(n2) memory though, and it is more computation. Depending on the case it can be faster due to parallelism (I did something similar to eliminate a tf.while_loop in TensorFlow to great success), but in this case it is actually slower:

x = list(range(1, 101))
a = 0.45
%timeit fun(x, a)
# 31 µs ± 85.2 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit fun_vec(x, a)
# 147 µs ± 2.42 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

So, there can be a non-loop version, but it is more of a curiosity than anything else.

Upvotes: 0

Graipher
Graipher

Reputation: 7186

As I wrote in an answer to basically the same question, you can't:

There is no other way (in general) except for an explicit for loop. This is because there is no way to parallelize this task across the rows (since every row depends on some other row).

What makes this even harder is that you can easily generate chaotic behavior, for example with the seemingly innocent looking logistic map: x_{n+1} = r * x_n * (1 - x_{n-1}).

You can only find a way around this if you manage to find a closed form, essentially eliminating the recurrence relation. But this has to be done for each recurrence relation and I am pretty sure you are not even guaranteed that a closed form exists...

Upvotes: 1

Related Questions