Reputation: 8660
Good evening,
I am trying to speed up the loop in this code. I have read through the numpy docs but to no avail. np.accumulate looks like it is almost what I need, but not quite.
What could I do to speed up the loop?
import numpy as np
N = 1000
AR_part = np.random.randn(N+1)
s2 = np.ndarray(N+1)
s2[0] = 1.0
beta = 1.3
old_s2 = s2[0]
for t in range( 1, N+1 ):
s2_t = AR_part[ t-1 ] + beta * old_s2
s2[t] = s2_t
old_s2 = s2_t
In response to Warren, I updated my code:
import numpy as np
from scipy.signal import lfilter, lfiltic
N = 1000
AR_part = np.random.randn(N+1)
beta = 1.3
def method1( AR_part):
s2 = np.empty_like(AR_part)
s2[0] = 1.0
old_s2 = s2[0]
for t in range( 1, N+1 ):
s2_t = AR_part[ t-1 ] + beta * old_s2
s2[t] = s2_t
old_s2 = s2_t
return s2
def method2( AR_part):
y = np.empty_like(AR_part)
b = np.array([0, 1])
a = np.array([1, -beta])
# Initial condition for the linear filter.
zi = lfiltic(b, a, [1.0], AR_part[:1])
y[:1] = 1.0
y[1:], zo = lfilter(b, a, AR_part[1:], zi=zi)
return y
s2 = method1( AR_part )
y = method2( AR_part )
np.alltrue( s2==y )
Timing the code:
%timeit method1( AR_part )
100 loops, best of 3: 1.63 ms per loop
%timeit method2( AR_part )
10000 loops, best of 3: 129 us per loop
That shows that Warren's method is over 10 times faster! Very impressive!
Upvotes: 2
Views: 425
Reputation: 114976
Your recurrence relation is linear, so it can be viewed as a linear filter. You can use scipy.signal.lfilter
to compute s2
. I recently answered a similar question here: python recursive vectorization with timeseries
Here's a script that shows how to use lfilter
to compute your series:
import numpy as np
from scipy.signal import lfilter, lfiltic
np.random.seed(123)
N = 4
AR_part = np.random.randn(N+1)
s2 = np.ndarray(N+1)
s2[0] = 1.0
beta = 1.3
old_s2 = s2[0]
for t in range( 1, N+1 ):
s2_t = AR_part[ t-1 ] + beta * old_s2
s2[t] = s2_t
old_s2 = s2_t
# Compute the result using scipy.signal.lfilter.
# Transfer function coefficients.
# `b` is the numerator, `a` is the denominator.
b = np.array([0, 1])
a = np.array([1, -beta])
# Initial condition for the linear filter.
zi = lfiltic(b, a, s2[:1], AR_part[:1])
# Apply lfilter to AR_part.
y = np.empty_like(AR_part)
y[:1] = s2[:1]
y[1:], zo = lfilter(b, a, AR_part[1:], zi=zi)
# Compare the results
print "s2 =", s2
print "y =", y
Output:
s2 = [ 1. 0.2143694 1.27602566 1.94181186 1.0180607 ]
y = [ 1. 0.2143694 1.27602566 1.94181186 1.0180607 ]
Upvotes: 5
Reputation: 3682
I agree with GHL that you won't get much more performance (though if N
was really big, and you were only computing some parts of the vector s2
, definitely use his method), but here's a different way to do what you're looking at:
import numpy as np
N = 1000
AR_part = np.random.randn(N+1)
beta = 1.3
def seq_gen(beta, constants, first_element = 1.0):
next_element = first_element
yield next_element
for j in constants:
next_element = j + beta * next_element
yield next_element
s2 = np.array([j for j in seq_gen(beta, AR_part, 1.0)])
Upvotes: 0
Reputation: 572
I am not sure there is much to do to speed up the loop... The only way I see would be to avoid the recursion, ie compute s2[t] "directly" for each t. But this is expensive as well...
You have
s2[t] = AR_part[t-1] + beta * s2[t-1]
= AR_part[t-1] + beta * (AR_part[t-2] + beta * s2[t-2])
= AR_part[t-1] + beta * AR_part[t-2] + beta^2 * s2[t-2]
= np.dot( AR[:t-1], beta_powers[-(t-1):] )
Where beta_powers contains [beta^1000, beta^999, ... 1.0]. You can create beta_powers this way:
np.power(beta, np,arange(1000))[::-1].
But I can't see a way to compute that stuff faster than what your loop does...
However you can rewrite it:
for t in range(N):
s2[t+1] = AR_part[t] + beta * s2[t]
Upvotes: 3