Numpy: calculate based on previous element?

Question

Say that I have array x and y:

x = numpy.array([1,2,3,4,5,6,7,8,9,10])  # actual content is the a result of another calculation step

There's a formula for y, and each element is based on the previous element, let i denote the index of y, each element is:

y[i] = y[i-1] * 2 + x[i]

When calculating the first element, let y[i-1] = 50. In other words, y should be:

[101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236]

How do I calculate y with numpy?

This is a non-homogeneous linear recurrence equation. Have you considered finding the analytic solution first, then seeing if it could be efficiently coded? — Ami Tavory, May 06 '15 at 05:43

score 11 · Accepted Answer · answered May 06 '15 at 06:42

Lets build a few of the items in your sequence:

y[0] = 2*y[-1] + x[0]
y[1] = 2*y[0] + x[1] = 4*y[-1] + 2*x[0] + x[1]
y[2] = 2*y[1] + x[2] = 8*y[-1] + 4*x[0] + 2*x[1] + x[2]
...
y[n] = 2**(n+1)*y[-1] + 2**n*x[0] + 2**(n-1)*x[1] + ... + x[n]

It may not be immediately obvious, but you can build the above sequence with numpy doing something like:

n = len(x)
y_1 = 50
pot = 2**np.arange(n-1, -1, -1)
y = np.cumsum(pot * x) / pot + y_1 * 2**np.arange(1, n+1)
>>> y
array([  101,   204,   411,   826,  1657,  3320,  6647, 13302, 26613, 53236])

The down side to this type of solutions is that they are not very general: a small change in your problem may render the whole approach useless. But whenever you can solve a problem with a little algebra, it is almost certainly going to beat any algorithmic approach by a far margin.

`numpy.cumsum()` is much faster than using a `for` loop, ~8x faster for my test application. — Steven C. Howell, Jul 23 '20 at 16:13

score 5 · Answer 2 · edited May 23 '17 at 12:06

5

If you need a recursive computation, if your y[i] should depend on the computed y[i-1] from the same run, then there seems to be no built-in solution in numpy, and you will need to compute it using a simple for loop:

y = np.empty(x.size)
last = 50
for i in range(x.size):
    y[i] = last = last * 2 + x[i]

See this question: Is a "for" loop necessary if elements of the a numpy vector are dependant upon the previous element?

Otherwise, you can implement your formula in one line using numpy:

y = np.concatenate(([50], y[:-1])) * 2 + x

Explanation:

y[:-1]

Creates a N-1-sized array: y_0, y_1, ... y_N-1.

np.concatenate(([50], y[:-1]))

Creates a N-sized array with the first element your starting value 50. So this expression basically is your y[i-1].

Then you can do the math element-wise using numpy array arithmetics.

edited May 23 '17 at 12:06

Community

1
1

answered May 06 '15 at 05:29

Ferdinand Beyer

64,979
15
154
145

Nice! And then change the first element of `y`. – Divakar May 06 '15 at 05:30
1

Shouldn't the output be `[101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236]`? Mismatch with second element I guess? Should be 204 and not 102 – Zero May 06 '15 at 05:32
3

This doesn't seem to handle the recurrence properly. It looks like it's designed to construct a whole new `y` array from an old `y` and `x`, while the desired behavior seems to be for the given relationship to hold within a single `y` array; the desired result would be a fixed point of this solution. – user2357112 May 06 '15 at 05:34
OK, I think I misunderstood the requirements then. If the OP wants a recursive formula, this answer won't work for him. (There was no example result when I wrote this answer) – Ferdinand Beyer May 06 '15 at 05:37
Indeed what I wanted is recursive, I didn't realize that and I expect numpy to be capable of this through some magic slicing. Thank you for the explanation. I'll stay with the for loop. – He Shiming May 06 '15 at 05:52

score 3 · Answer 3 · answered Dec 10 '18 at 12:13

Perhaps the fastest and most concise way is to use scipy.signal.lfilter, which implements exactly the kind of recursive relationship you described:

from scipy.signal import lfilter
import numpy as np

x = np.array([1,2,3,4,5,6,7,8,9,10])

b = [1., 0.]
a = [1., -2.]
zi = np.array([2*50])  # initial condition
y, _ = lfilter(b, a, x, zi=zi)

Result will be np.float64, but you can cast to e.g. np.int32 if that's what you need:

>>> y.astype(np.int32)
array([  101,   204,   411,   826,  1657,  3320,  6647, 13302, 26613, 53236])

score -3 · Answer 4 · answered Nov 08 '16 at 14:56

-3

This is how to do it with numpy:

import numpy as np
x = np.array([ 1, 2, 3, 4, 5, 6, 7, 8 ,9, 10 ])
y = np.array([ 50 ])
for i in np.arange(len(x)):
    y = np.append(
                  y,
                  ( y[-1] * 2 + x[i] )
                  )
y = y[1:]

print(y)

answered Nov 08 '16 at 14:56

KirkLab

9
4

never use for loops, when it comes to performance in python – Helmut Jun 20 '23 at 06:23

Numpy: calculate based on previous element?

4 Answers4

Linked