1

I have a recurrent relation that looks like the following:

f_0 = s_0
f_1 = alpha * s_1 + (1 - alpha) * f_0
f_2 = alpha * s_2 + (1 - alpha) * f_1
...
f_n = alpha * s_n + (1 - alpha) * f_n-1

this is easy enough to code up as a for-loop (example code)

import numpy as np
## values for example only ##
s = np.linspace(0, 10, 3000)
alpha = 1 / 5
f = s.copy()
for i in range(1, len(s)):
    f[i] = alpha * s[i] + (1 - alpha) * f[i-1]

my struggle is finding a way to vectorize this relation. I've tried being clever with slice indexing, but something like

f[1:] = alpha * s[1:] + (1 - alpha) * f[:-1]

is clearly not right, as it doesn't update f.

Is there an obvious way to vectorize this recurrence relation? I'd rather avoid this for-loop if I can, as my code will have to do this thousands of times as part of an optimization routine.

Thanks!

Update

Thanks for the feedback. I tried a naive numba implementation @jit and it actually made the loop 50x faster?! Being new to numba, is there any way to further tune these results?

not link
  • 858
  • 1
  • 16
  • 29

1 Answers1

0

You could change this into a matrix equation. If you write f and s as numpy arrays (1d) and write a 2d numpy matrix M, such that:

F = numpy.dot(M,S)

then you only have to create the matrix once (it only depends on alpha), and can keep using it on different S. On first inspection the matrix looks something like:

|1,0,0,0|
|(1-alpha),alpha,0,0|
|(1-alpha)^2,(1-alpha),alpha,0|
|(1-alpha)^3,(1-alpha)^2,(1-alpha),alpha|
ergoregion
  • 1