我可以使用 numpy 来加速这个循环吗? [英] Can I use numpy to speed this loop?

问题描述

晚上好，

我正在尝试加速此代码中的循环.我已经通读了 numpy 文档，但无济于事.np.accumulate 看起来几乎是我需要的，但不完全是.

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

为了回应 Warren，我更新了我的代码:

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 )

定时代码:

%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

这表明沃伦的方法快了 10 倍以上！非常令人印象深刻！

That shows that Warren's method is over 10 times faster! Very impressive!

推荐答案

您的递推关系是线性的，因此可以将其视为线性过滤器.您可以使用 scipy.signal.lfilter 来计算 s2.我最近在这里回答了一个类似的问题:python recursive vectorization with timeseries

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

这是一个脚本，展示了如何使用 lfilter 来计算您的系列:

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

输出:

s2 = [ 1.          0.2143694   1.27602566  1.94181186  1.0180607 ]
y  = [ 1.          0.2143694   1.27602566  1.94181186  1.0180607 ]

这篇关于我可以使用 numpy 来加速这个循环吗?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！