1

I am learning to vectorize my code, which is currently suffering from severe time issues. My goal is to get rid of all the loops that slow down the process. This has worked fine whenever I did not need to know the result of step i to calculate step i+1.

At one point, this does not hold true anymore. I tried to break down the problem to this simplified snippet:

x2 = np.arange(10)
f1, sumint = 0.0, 0.0

for i in x2:
    f2 = np.exp(x2[i]) + f1
    sumint += f2
    f1 = f2

The term np.exp(x2[i]) can of course be calculated for the whole matrix. But f1 is a term which carries the result of i-1 which I do not know in advance (or, do I?)

I am looking for a numpy-solution to this problem, because the real code involves more complex calculations and more of the variables like f1. I am aware of numpy.vectorize but the docs say that this does not speed up the code, but rather assists readibility.

Edit: Maybe this example is more suitable to show my actual case...

x2 = np.random.rand(10)
y2 = np.random.rand(10)
f2 = np.random.rand(10)
sumint, f1, x1, y1 = 0, 0, 0, 0

for i in range(10):
    sumint += (f2[i] - f1) * (x2[i] - x1) / (y2[i] - y1)
    x1 = x2[i]
    y1 = y2[i]
    f1 = f2[i]

Edit 2: Sorry for that, I found the solution. By breaking down the problem to its core I found that what I thought of a recursion was just a SHIFT. Solution to the problem was:

f1 = np.zeros(10)
x1 = np.zeros(10)
y1 = np.zeros(10)

x1[1:] = x2[:-1]
y1[1:] = y2[:-1]
f1[1:] = f2[:-1]

sumint = (f2 - f1) * (x2 - x1) / (y2 - y1)

which I can now apply to the real case.

offeltoffel
  • 2,691
  • 2
  • 21
  • 35
  • 1
    Does this answer your question? [Is it possible to vectorize recursive calculation of a NumPy array where each element depends on the previous one?](https://stackoverflow.com/questions/4407984/is-it-possible-to-vectorize-recursive-calculation-of-a-numpy-array-where-each-el) – Iguananaut Dec 23 '19 at 09:18
  • Although looking more closely at your case, I can't see how it's any different from just a sum of exponentials. Unless you've simplified your actual case, I'm missing the need for a recursive definition here... – Iguananaut Dec 23 '19 at 09:23
  • 1
    Indeed, the code you've written is exactly equivalent to `s = np.exp(np.arange(10)); f1 = np.sum(s); sumint = np.sum(np.cumsum(s))`. But in a general case you can't vectorize a recursively defined function. Here you've just taken a trivial summation and made it look more complicated than it is. But maybe your real problem is more complicated. – Iguananaut Dec 23 '19 at 09:27
  • I am kind of confused right now and may have asked something quite dumb ... I totally see how my example does not really represent a complicated recursive calculation. I found a more accurate example for my case. Would you have a look at it? Thanks for the hint that recursive calculations cannot be vectorized with NumPy! – offeltoffel Dec 23 '19 at 09:45
  • 1
    The question you should be asking yourself is: Is there a "closed form" expression for your sum. If you can find that then it can easily be vectorized. It may be quite easy to find with a little bit of algebra, or in some cases it might be more difficult. – Iguananaut Dec 23 '19 at 09:54
  • I think I just found the solution. I don't know what "closed form" is, but I can figure what it means. And had the feeling it should somehow be manageable. Thanks for your hint! – offeltoffel Dec 23 '19 at 09:56
  • 1
    Yup, I was about to tell you the same, that in this case you just need to sum a shift. The reason you can do it in this case as because you only depend trivially on the previous values in the array. If it was something possibly non-linear like `x_3 = f(x_2) - f(x_1)` then this might not be possible. – Iguananaut Dec 23 '19 at 10:01

1 Answers1

1

You can do this using 1 calls to np.cumsum and 1 call to np.exp:

sumint = np.sum(np.cumsum(np.exp(x2)))

If you also want the values of f1 or f2:

f1 = f2 = sum(np.exp(x2))

This gives the same values as your code, but is vectorized.

Rohan Saxena
  • 3,133
  • 2
  • 16
  • 34