4

I have an array:

a = np.array([2,3,5,8,3,5])

What is the most efficient (vectorized) way to calculate an array where each resulting element is (Pseudocode):

result[0] = a[0]
for i > 0:
    result[i] = result[i-1] + (a[i] - result[i-1]) * factor

I could do this with the following inefficient code (factor = 0.5):

a = np.array([2,3,5,8,3,5])
result = np.array([a[0]])
for k in a[1:]:
   result = np.append(result, result[-1]+(k-result[-1])*0.5)

The result of this damping function would be:

array([ 2.,  2.5,  3.75,  5.875,  4.4375,  4.71875])
Lafexlos
  • 7,618
  • 5
  • 38
  • 53
leofields
  • 659
  • 1
  • 7
  • 12
  • 2
    Not every recurrence relation can be vectorized. Whether yours can is really more of a math question than a programming question (i.e., can you convert the recurrence relation into an explicit formula for `f(n)`). – BrenBarn Dec 30 '15 at 07:32
  • Theoretically, this is O(n) and you can't do better than that, if the requirement is to calculate all of the n values. – Selçuk Cihan Dec 30 '15 at 08:07

2 Answers2

6

You a looking for Haskell's scanl1 alternative in Python (Haskell example):

Prelude> scanl1 (\a  b -> a + (b - a) * 0.5) [2, 3, 5, 8, 3, 5]
[2.0,2.5,3.75,5.875,4.4375,4.71875]

There is an accumulate function in itertools module:

In [1]: import itertools

In [2]: itertools.accumulate([2, 3, 5, 8, 3, 5], lambda a, b: a + (b - a) * 0.5)
Out[2]: <itertools.accumulate at 0x7f1fc1fc1608>

In [3]: list(itertools.accumulate([2, 3, 5, 8, 3, 5], lambda a, b: a + (b - a) * 0.5))
Out[3]: [2, 2.5, 3.75, 5.875, 4.4375, 4.71875]

With NumPy you may use numpy.ufunc.accumulate function, however, according to this answer, there is a bug in the implementation, that is why we should use a cast. Unfortunately, I'm not very familiar with NumPy, and, probably, there is a better way:

In [9]: import numpy as np

In [10]: uf = np.frompyfunc(lambda a, b: a + (b - a) * 0.5, 2, 1)

In [11]: uf.accumulate([2,3,5,8,3,5], dtype=np.object).astype(np.float)
Out[11]: array([ 2.     ,  2.5    ,  3.75   ,  5.875  ,  4.4375 ,  4.71875])
Community
  • 1
  • 1
awesoon
  • 32,469
  • 11
  • 74
  • 99
1

I would like to post how the @soon code works, or how to implement it with reduce:

def scanl(f, v):
    return reduce(lambda (l, v1), v:(l+[f(v1, v)], f(v1, v)), v[1:],  ([v[0]], v[0]))[0]

>>> scanl(lambda a, b: a + (b - a) * 0.5,[2, 3, 5, 8, 3, 5])
[2, 2.5, 3.75, 5.875, 4.4375, 4.71875]

Its not the best performance nor the cutier one but gives you an idea of how to do it.

Netwave
  • 40,134
  • 6
  • 50
  • 93