recursion - python recursive vectorization with timeseries

Question

I have a Timeseries (s) which need to be processed recursively to get a timeseries result (res). Here is my sample code:

res=s.copy()*0  
res[1]=k # k is a constant  
for i in range(2,len(s)):  
    res[i]=c1*(s[i]+s[i-1])/2 +c2*res[i-1]+c3*res[i-2]

where c1,c2,c3 are constants. It works properly but I'd like to use vectorization and I tried with:

res[2:]=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]

but I get "ValueError: operands could not be broadcast together with shapes (1016) (1018) "
if I try with

res=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]

doesn't give any error, but I don't get a correct result, because res[0] and res[1] have to be initialized before the calculation will take place. Is there a way to process it with vectorization?
Any help will be appreciated, thanks!

Update this with the values you are using for the constants. If you make it copy-paste runnable, you'll probably get more help. Also - do you know about numpy? — mrKelley
Is this because you want to speed things up? (memoization could speed up stuff like this if you rewrite it a bit i guess?) — usethedeathstar
For such problems I like using numexpr (code.google.com/p/numexpr) often improves speed significantly with little effort. — Dietrich

Warren Weckesser Warren Weckesser · Accepted Answer · 2014-01-24T17:22:58

This expression

    res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]

says that res is the output of a linear filter (or ARMA process) with input s. Several libraries have functions for computing this. Here's how you can use the scipy function scipy.signal.lfilter.

From inspection of the recurrence relation, we get the coefficients of the numerator (b) and denominator (a) of the filter's transfer function:

b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

We'll also need an appropriate initial condition for lfilter to handle res[:2] == [0, k]. For this, we use scipy.signal.lfiltic:

zi = lfiltic(b, a, [k, 0], x=s[1::-1])

In the simplest case, one would call lfilter like this:

y = lfilter(b, a, s)

With an initial condition zi, we use:

y, zo = lfilter(b, a, s, zi=zi)

However, to exactly match the calculation provided in the question, we need the output y to start with [0, k]. So we'll allocate an array y, initialize the first two elements with [0, k], and assign the output of lfilter to y[2:]:

y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)

Here's a complete script with the original loop and with lfilter:

import numpy as np
from scipy.signal import lfilter, lfiltic


c1 = 0.125
c2 = 0.5
c3 = 0.25

np.random.seed(123)
s = np.random.rand(8)
k = 3.0

# Original version (edited lightly)

res = np.zeros_like(s)
res[1] = k  # k is a constant  
for i in range(2, len(s)):  
    res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]


# Using scipy.signal.lfilter

# Coefficients of the filter's transfer function.
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

# Create the initial condition of the filter such that
#     y[:2] == [0, k]
zi = lfiltic(b, a, [k, 0], x=s[1::-1])

y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)

np.set_printoptions(precision=5)
print "res:", res
print "y:  ", y

The output is:

res: [ 0.       3.       1.53206  1.56467  1.24477  1.08496  0.94142  0.84605]
y:   [ 0.       3.       1.53206  1.56467  1.24477  1.08496  0.94142  0.84605]

lfilter accepts an axis argument, so you can filter an array of signals with a single call. lfiltic does not have an axis argument, so setting up the initial conditions requires a loop. The following script shows an example.

import numpy as np
from scipy.signal import lfilter, lfiltic
import matplotlib.pyplot as plt


# Parameters
c1 = 0.2
c2 = 1.1
c3 = -0.5
k = 1

# Create an array of signals for the demonstration.
np.random.seed(123)
nsamples = 50
nsignals = 4
s = np.random.randn(nsamples, nsignals)

# Coefficients of the filter's transfer function.
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

# Create the initial condition of the filter for each signal
# such that
#     y[:2] == [0, k]
# We need a loop here, because lfiltic is not vectorized.
zi = np.empty((2, nsignals))
for i in range(nsignals):
    zi[:, i] = lfiltic(b, a, [k, 0], x=s[1::-1, i])

# Create the filtered signals.
y = np.empty_like(s)
y[:2, :] = np.array([0, k]).reshape(-1, 1)
y[2:, :], zo = lfilter(b, a, s[2:], zi=zi, axis=0)

# Plot the filtered signals.
plt.plot(y, linewidth=2, alpha=0.6)
ptp = y.ptp()
plt.ylim(y.min() - 0.05*ptp, y.max() + 0.05*ptp)
plt.grid(True)
plt.show()

Plot:

filtered signals

recursion - python recursive vectorization with timeseries

1 Answers