I think I have finally cracked it!

Here's a vectorized version of numpy_ewma function that's claimed to be producing the correct results from @RaduS's post -

def numpy_ewma_vectorized(data, window):



    alpha = 2 /(window + 1.0)

    alpha_rev = 1-alpha



    scale = 1/alpha_rev

    n = data.shape[0]



    r = np.arange(n)

    scale_arr = scale**r

    offset = data[0]*alpha_rev**(r+1)

    pw0 = alpha*alpha_rev**(n-1)



    mult = data*pw0*scale_arr

    cumsums = mult.cumsum()

    out = offset + cumsums*scale_arr[::-1]

    return out

Further boost

We can boost it further with some code re-use, like so -

def numpy_ewma_vectorized_v2(data, window):



    alpha = 2 /(window + 1.0)

    alpha_rev = 1-alpha

    n = data.shape[0]



    pows = alpha_rev**(np.arange(n+1))



    scale_arr = 1/pows[:-1]

    offset = data[0]*pows[1:]

    pw0 = alpha*alpha_rev**(n-1)



    mult = data*pw0*scale_arr

    cumsums = mult.cumsum()

    out = offset + cumsums*scale_arr[::-1]

    return out

Runtime test

Let's time these two against the same loopy function for a big dataset.

In [97]: data = np.random.randint(2,9,(5000))

    ...: window = 20

    ...:



In [98]: np.allclose(numpy_ewma(data, window), numpy_ewma_vectorized(data, window))

Out[98]: True



In [99]: np.allclose(numpy_ewma(data, window), numpy_ewma_vectorized_v2(data, window))

Out[99]: True



In [100]: %timeit numpy_ewma(data, window)

100 loops, best of 3: 6.03 ms per loop



In [101]: %timeit numpy_ewma_vectorized(data, window)

1000 loops, best of 3: 665 µs per loop



In [102]: %timeit numpy_ewma_vectorized_v2(data, window)

1000 loops, best of 3: 357 µs per loop



In [103]: 6030/357.0

Out[103]: 16.89075630252101

There is around a 17 times speedup!

Here is an implementation using NumPy that is equivalent to using df.ewm(alpha=alpha).mean(). After reading the documentation, it is just a few matrix operations. The trick is constructing the right matrices.

It is worth noting that because we are creating float matrices, you can quickly eat through your memory if the input array is too large.

import pandas as pd

import numpy as np



def ewma(x, alpha):

    '''

    Returns the exponentially weighted moving average of x.



    Parameters:

    -----------

    x : array-like

    alpha : float {0 <= alpha <= 1}



    Returns:

    --------

    ewma: numpy array

          the exponentially weighted moving average

    '''

    # Coerce x to an array

    x = np.array(x)

    n = x.size



    # Create an initial weight matrix of (1-alpha), and a matrix of powers

    # to raise the weights by

    w0 = np.ones(shape=(n,n)) * (1-alpha)

    p = np.vstack([np.arange(i,i-n,-1) for i in range(n)])



    # Create the weight matrix

    w = np.tril(w0**p,0)



    # Calculate the ewma

    return np.dot(w, x[::np.newaxis]) / w.sum(axis=1)

Let's test its:

alpha = 0.55

x = np.random.randint(0,30,15)

df = pd.DataFrame(x, columns=['A'])

df.ewm(alpha=alpha).mean()



# returns:

#             A

# 0   13.000000

# 1   22.655172

# 2   20.443268

# 3   12.159796

# 4   14.871955

# 5   15.497575

# 6   20.743511

# 7   20.884818

# 8   24.250715

# 9   18.610901

# 10  17.174686

# 11  16.528564

# 12  17.337879

# 13   7.801912

# 14  12.310889



ewma(x=x, alpha=alpha)



# returns:

# array([ 13.        ,  22.65517241,  20.44326778,  12.1597964 ,

#        14.87195534,  15.4975749 ,  20.74351117,  20.88481763,

#        24.25071484,  18.61090129,  17.17468551,  16.52856393,

#        17.33787888,   7.80191235,  12.31088889])

Fastest EWMA 23x pandas

The question is strictly asking for a numpy solution, however, it seems that the OP was actually just after a pure numpy solution to speed up runtime.

I solved a similar problem but instead looked towards numba.jit which massively speeds the compute time

In [24]: a = np.random.random(10**7)

    ...: df = pd.Series(a)

In [25]: %timeit numpy_ewma(a, 10)               # /a/42915307/4013571

    ...: %timeit df.ewm(span=10).mean()          # pandas

    ...: %timeit numpy_ewma_vectorized_v2(a, 10) # best w/o numba: /a/42926270/4013571

    ...: %timeit _ewma(a, 10)                    # fastest accurate (below)

    ...: %timeit _ewma_infinite_hist(a, 10)      # fastest overall (below)

4.14 s ± 116 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

991 ms ± 52.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 

396 ms ± 8.39 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

181 ms ± 1.01 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)   

39.6 ms ± 979 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Scaling down to smaller arrays of a = np.random.random(100) (results in the same order)

41.6 µs ± 491 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

945 ms ± 12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

16 µs ± 93.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

1.66 µs ± 13.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

1.14 µs ± 5.57 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

It is also worth pointing out that my functions below are identically aligned to the pandas (see the examples in docstr), whereas a few of the answers here take various different approximations. For example,

In [57]: print(pd.DataFrame([1,2,3]).ewm(span=2).mean().values.ravel())

    ...: print(numpy_ewma_vectorized_v2(np.array([1,2,3]), 2))

    ...: print(numpy_ewma(np.array([1,2,3]), 2))

[1.         1.75       2.61538462]

[1.         1.66666667 2.55555556]

[1.         1.18181818 1.51239669]

The source code which I have documented for my own library

import numpy as np

from numba import jit

from numba import float64

from numba import int64





@jit((float64[:], int64), nopython=True, nogil=True)

def _ewma(arr_in, window):

    r"""Exponentialy weighted moving average specified by a decay ``window``

    to provide better adjustments for small windows via:



        y[t] = (x[t] + (1-a)*x[t-1] + (1-a)^2*x[t-2] + ... + (1-a)^n*x[t-n]) /

               (1 + (1-a) + (1-a)^2 + ... + (1-a)^n).



    Parameters

    ----------

    arr_in : np.ndarray, float64

        A single dimenisional numpy array

    window : int64

        The decay window, or 'span'



    Returns

    -------

    np.ndarray

        The EWMA vector, same length / shape as ``arr_in``



    Examples

    --------

    >>> import pandas as pd

    >>> a = np.arange(5, dtype=float)

    >>> exp = pd.DataFrame(a).ewm(span=10, adjust=True).mean()

    >>> np.array_equal(_ewma_infinite_hist(a, 10), exp.values.ravel())

    True

    """

    n = arr_in.shape[0]

    ewma = np.empty(n, dtype=float64)

    alpha = 2 / float(window + 1)

    w = 1

    ewma_old = arr_in[0]

    ewma[0] = ewma_old

    for i in range(1, n):

        w += (1-alpha)**i

        ewma_old = ewma_old*(1-alpha) + arr_in[i]

        ewma[i] = ewma_old / w

    return ewma





@jit((float64[:], int64), nopython=True, nogil=True)

def _ewma_infinite_hist(arr_in, window):

    r"""Exponentialy weighted moving average specified by a decay ``window``

    assuming infinite history via the recursive form:



        (2) (i)  y[0] = x[0]; and

            (ii) y[t] = a*x[t] + (1-a)*y[t-1] for t>0.



    This method is less accurate that ``_ewma`` but

    much faster:



        In [1]: import numpy as np, bars

           ...: arr = np.random.random(100000)

           ...: %timeit bars._ewma(arr, 10)

           ...: %timeit bars._ewma_infinite_hist(arr, 10)

        3.74 ms ± 60.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

        262 µs ± 1.54 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)



    Parameters

    ----------

    arr_in : np.ndarray, float64

        A single dimenisional numpy array

    window : int64

        The decay window, or 'span'



    Returns

    -------

    np.ndarray

        The EWMA vector, same length / shape as ``arr_in``



    Examples

    --------

    >>> import pandas as pd

    >>> a = np.arange(5, dtype=float)

    >>> exp = pd.DataFrame(a).ewm(span=10, adjust=False).mean()

    >>> np.array_equal(_ewma_infinite_hist(a, 10), exp.values.ravel())

    True

    """

    n = arr_in.shape[0]

    ewma = np.empty(n, dtype=float64)

    alpha = 2 / float(window + 1)

    ewma[0] = arr_in[0]

    for i in range(1, n):

        ewma[i] = arr_in[i] * alpha + ewma[i-1] * (1 - alpha)

    return ewma

Given alpha and windowSize, here's an approach to simulate the corresponding behavior on NumPy -

def numpy_ewm_alpha(a, alpha, windowSize):

    wghts = (1-alpha)**np.arange(windowSize)

    wghts /= wghts.sum()

    out = np.full(df.shape[0],np.nan)

    out[windowSize-1:] = np.convolve(a,wghts,'valid')

    return out

Sample runs for verification -

In [54]: alpha = 0.55

    ...: windowSize = 20

    ...: 



In [55]: df = pd.DataFrame(np.random.randint(2,9,(100)))



In [56]: out0 = df.ewm(alpha = alpha, min_periods=windowSize).mean().as_matrix().ravel()

    ...: out1 = numpy_ewm_alpha(df.values.ravel(), alpha = alpha, windowSize = windowSize)

    ...: print "Max. error : " + str(np.nanmax(np.abs(out0 - out1)))

    ...: 

Max. error : 5.10531254605e-07



In [57]: alpha = 0.75

    ...: windowSize = 30

    ...: 



In [58]: out0 = df.ewm(alpha = alpha, min_periods=windowSize).mean().as_matrix().ravel()

    ...: out1 = numpy_ewm_alpha(df.values.ravel(), alpha = alpha, windowSize = windowSize)

    ...: print "Max. error : " + str(np.nanmax(np.abs(out0 - out1)))



Max. error : 8.881784197e-16

Runtime test on bigger dataset -

In [61]: alpha = 0.55

    ...: windowSize = 20

    ...: 



In [62]: df = pd.DataFrame(np.random.randint(2,9,(10000)))



In [63]: %timeit df.ewm(alpha = alpha, min_periods=windowSize).mean()

1000 loops, best of 3: 851 µs per loop



In [64]: %timeit numpy_ewm_alpha(df.values.ravel(), alpha = alpha, windowSize = windowSize)

1000 loops, best of 3: 204 µs per loop

Further boost

For further performance boost we could avoid the initialization with NaNs and instead use the array outputted from np.convolve, like so -

def numpy_ewm_alpha_v2(a, alpha, windowSize):

    wghts = (1-alpha)**np.arange(windowSize)

    wghts /= wghts.sum()

    out = np.convolve(a,wghts)

    out[:windowSize-1] = np.nan

    return out[:a.size]  

Timings -

In [117]: alpha = 0.55

     ...: windowSize = 20

     ...: 



In [118]: df = pd.DataFrame(np.random.randint(2,9,(10000)))



In [119]: %timeit numpy_ewm_alpha(df.values.ravel(), alpha = alpha, windowSize = windowSize)

1000 loops, best of 3: 204 µs per loop



In [120]: %timeit numpy_ewm_alpha_v2(df.values.ravel(), alpha = alpha, windowSize = windowSize)

10000 loops, best of 3: 195 µs per loop

Here is another solution O came up with in the meantime. It is about four times faster than the pandas solution.

def numpy_ewma(data, window):

    returnArray = np.empty((data.shape[0]))

    returnArray.fill(np.nan)

    e = data[0]

    alpha = 2 / float(window + 1)

    for s in range(data.shape[0]):

        e =  ((data[s]-e) *alpha ) + e

        returnArray[s] = e

    return returnArray

I used this formula as a starting point. I am sure that this can be improved even more, but it is at least a starting point.

@Divakar's answer seems to cause overflow when dealing with

numpy_ewma_vectorized(np.random.random(500000), 10)

What I have been using is:

def EMA(input, time_period=10): # For time period = 10

    t_ = time_period - 1

    ema = np.zeros_like(input,dtype=float)

    multiplier = 2.0 / (time_period + 1)

    #multiplier = 1 - multiplier

    for i in range(len(input)):

        # Special Case

        if i > t_:

            ema[i] = (input[i] - ema[i-1]) * multiplier + ema[i-1]

        else:

            ema[i] = np.mean(input[:i+1])

    return ema

However, this is way slower than the panda solution:

from pandas import ewma as pd_ema

def EMA_fast(X, time_period = 10):

    out = pd_ema(X, span=time_period, min_periods=time_period)

    out[:time_period-1] = np.cumsum(X[:time_period-1]) / np.asarray(range(1,time_period))

    return out

This answer may seem irrelevant. But, for those who also need to calculate the exponentially weighted variance (and also standard deviation) with NumPy, the following solution will be useful:

import numpy as np



def ew(a, alpha, winSize):

    _alpha = 1 - alpha

    ws = _alpha ** np.arange(winSize)

    w_sum = ws.sum()

    ew_mean = np.convolve(a, ws)[winSize - 1] / w_sum

    bias = (w_sum ** 2) / ((w_sum ** 2) - (ws ** 2).sum())

    ew_var = (np.convolve((a - ew_mean) ** 2, ws)[winSize - 1] / w_sum) * bias

    ew_std = np.sqrt(ew_var)

    return (ew_mean, ew_var, ew_std)

Building on top of Divakar's great answer, here is an implementation which corresponds to the adjust=True flag of the pandas function, i.e. using weights rather than recursion.

def numpy_ewma(data, window):

    alpha = 2 /(window + 1.0)

    scale = 1/(1-alpha)

    n = data.shape[0]

    scale_arr = (1-alpha)**(-1*np.arange(n))

    weights = (1-alpha)**np.arange(n)

    pw0 = (1-alpha)**(n-1)

    mult = data*pw0*scale_arr

    cumsums = mult.cumsum()

    out = cumsums*scale_arr[::-1] / weights.cumsum()



    return out

Thanks to @Divakar's solution and that is really fast. However, it does cause overflow problem which was pointed out by @Danny. The function doesn't return correct answers when the length is greater than 13835 or so at my end.

The following is my solution based on Divakar's solution and pandas.ewm().mean()

def numpy_ema(data, com=None, span=None, halflife=None, alpha=None):

"""Summary

Calculate ema with automatically-generated alpha. Weight of past effect

decreases as the length of window increasing.



# these functions reproduce the pandas result when the flag adjust=False is set.

References:

https://stackoverflow.com/questions/42869495/numpy-version-of-exponential-weighted-moving-average-equivalent-to-pandas-ewm



Args:

    data (TYPE): Description

    com (float, optional): Specify decay in terms of center of mass, alpha=1/(1+com), for com>=0

    span (float, optional): Specify decay in terms of span, alpha=2/(span+1), for span>=1

    halflife (float, optional): Specify decay in terms of half-life, alpha=1-exp(log(0.5)/halflife), for halflife>0

    alpha (float, optional): Specify smoothing factor alpha directly, 0<alpha<=1



Returns:

    TYPE: Description



Raises:

    ValueError: Description

"""

n_input = sum(map(bool, [com, span, halflife, alpha]))

if n_input != 1:

    raise ValueError(

        'com, span, halflife, and alpha are mutually exclusive')



nrow = data.shape[0]

if np.isnan(data).any() or (nrow > 13835) or (data.ndim == 2):

    df = pd.DataFrame(data)

    df_ewm = df.ewm(com=com, span=span, halflife=halflife,

                    alpha=alpha, adjust=False)

    out = df_ewm.mean().values.squeeze()

else:

    if com:

        alpha = 1 / (1 + com)

    elif span:

        alpha = 2 / (span + 1.0)

    elif halflife:

        alpha = 1 - np.exp(np.log(0.5) / halflife)



    alpha_rev = 1 - alpha

    pows = alpha_rev**(np.arange(nrow + 1))



    scale_arr = 1 / pows[:-1]

    offset = data[0] * pows[1:]

    pw0 = alpha * alpha_rev**(nrow - 1)



    mult = data * pw0 * scale_arr



    cumsums = np.cumsum(mult)

    out = offset + cumsums * scale_arr[::-1]

return out

Here's my implementation for 1D input arrays with infinite window size. As it uses large numbers, it works only with input arrays with elements of absolute value < 1e16, when using float32, but that should normally be the case.

The idea is to reshape the input array into slices of a limited length, so that no overflow occurs, and then doing the ewm calculation with each slice separately.

def ewm(x, alpha):

    """

    Returns the exponentially weighted mean y of a numpy array x with scaling factor alpha

    y[0] = x[0]

    y[j] = (1. - alpha) * y[j-1] + alpha * x[j],  for j > 0



    x -- 1D numpy array

    alpha -- float

    """



    n = int(-100. / np.log(1.-alpha)) # Makes sure that the first and last elements in f are very big and very small (about 1e22 and 1e-22)

    f = np.exp(np.arange(1-n, n, 2) * (0.5 * np.log(1. - alpha))) # Scaling factor for each slice

    tmp = (np.resize(x, ((len(x) + n - 1) // n, n)) / f * alpha).cumsum(axis=1) * f # Get ewm for each slice of length n



    # Add the last value of each previous slice to the next slice with corresponding scaling factor f and return result

    return np.resize(tmp + np.tensordot(np.append(x[0], np.roll(tmp.T[n-1], 1)[1:]), f * ((1. - alpha) / f[0]), axes=0), len(x))