快速累加和幂运算符
Fast cumulative sum and power operator
我有一个预测算法,它使用以下代码计算时间序列在给定范围内的趋势:
import numpy as np
horizon = 91
phi = 0.2
trend = -0.004
trend_up_to_horizon = np.cumsum(phi ** np.arange(horizon) + 1) * self.trend
在此示例中,前两个 trend_up_horizon
值是:
array([-0.008 , -0.0128])
是否有计算速度更快的方法来实现此目的?目前这需要很长时间,因为我猜使用 np.cumsum
方法和 **
运算符很昂贵。
感谢您的帮助
您可以使用 Cython 使其更快一点,但速度并不快
运行 %timeit
在你的基础上 np.cumsum(phi ** np.arange(horizon) + 1) * trend
说我的笔记本电脑需要 17.5µs,这并不多
具有等效功能的 Cython 版本是:
import numpy as np
cimport numpy as np
cimport cython
@cython.boundscheck(False)
def do_cumsum(size_t horizon, double phi, double trend):
cdef np.ndarray[double, ndim=1] out = np.empty(horizon, dtype=np.float)
cdef double csum = 0
cdef int i
for i in range(horizon):
csum += phi ** i + 1
out[i] = csum * trend
return out
将 do_cumsum(horizon, phi, trend)
的时间减少到 6.9µs,而如果我切换到单个 precision/32bit 浮动,这将减少到 4.5µs
就是说,微秒并不多,您最好将精力集中在其他地方
您可以更快地执行此操作。正如您已经假设的(不必要的)电力运营商是这里的主要问题。
除此之外,Numpy 没有 power(float64,int64) 的特殊实现,其中指数是一个小的正整数。相反,Numpy 总是计算 power(float64,float64),这是一项复杂得多的任务。
Numba 和 Numexpr 有一个针对简单大小写 power(float64,int64) 的特殊实现,所以让我们先尝试一下。
第一种方法
import numpy as np
import numba as nb
horizon = 91
phi = 0.2
trend = -0.004
@nb.njit()
def pow_cumsum(horizon,phi,trend):
out=np.empty(horizon)
csum=0.
for i in range(horizon):
csum+=phi**i+1
out[i]=csum*trend
return out
之前说过直接计算幂是不必要的,可以重写算法来完全避免这种情况。
第二种方法
@nb.njit()
def pow_cumsum_2(horizon,phi,trend):
out=np.empty(horizon)
out[0]=2.*trend
TMP=2.
val=phi
for i in range(horizon-1):
TMP=(val+1+TMP)
out[i+1]=TMP*trend
val*=phi
return out
计时
%timeit np.cumsum(phi ** np.arange(horizon) + 1) * trend
7.44 µs ± 89.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit pow_cumsum(horizon,phi,trend)
882 ns ± 4.91 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit pow_cumsum_2(horizon,phi,trend)
559 ns ± 3.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
我有一个预测算法,它使用以下代码计算时间序列在给定范围内的趋势:
import numpy as np
horizon = 91
phi = 0.2
trend = -0.004
trend_up_to_horizon = np.cumsum(phi ** np.arange(horizon) + 1) * self.trend
在此示例中,前两个 trend_up_horizon
值是:
array([-0.008 , -0.0128])
是否有计算速度更快的方法来实现此目的?目前这需要很长时间,因为我猜使用 np.cumsum
方法和 **
运算符很昂贵。
感谢您的帮助
您可以使用 Cython 使其更快一点,但速度并不快
运行 %timeit
在你的基础上 np.cumsum(phi ** np.arange(horizon) + 1) * trend
说我的笔记本电脑需要 17.5µs,这并不多
具有等效功能的 Cython 版本是:
import numpy as np
cimport numpy as np
cimport cython
@cython.boundscheck(False)
def do_cumsum(size_t horizon, double phi, double trend):
cdef np.ndarray[double, ndim=1] out = np.empty(horizon, dtype=np.float)
cdef double csum = 0
cdef int i
for i in range(horizon):
csum += phi ** i + 1
out[i] = csum * trend
return out
将 do_cumsum(horizon, phi, trend)
的时间减少到 6.9µs,而如果我切换到单个 precision/32bit 浮动,这将减少到 4.5µs
就是说,微秒并不多,您最好将精力集中在其他地方
您可以更快地执行此操作。正如您已经假设的(不必要的)电力运营商是这里的主要问题。
除此之外,Numpy 没有 power(float64,int64) 的特殊实现,其中指数是一个小的正整数。相反,Numpy 总是计算 power(float64,float64),这是一项复杂得多的任务。
Numba 和 Numexpr 有一个针对简单大小写 power(float64,int64) 的特殊实现,所以让我们先尝试一下。
第一种方法
import numpy as np
import numba as nb
horizon = 91
phi = 0.2
trend = -0.004
@nb.njit()
def pow_cumsum(horizon,phi,trend):
out=np.empty(horizon)
csum=0.
for i in range(horizon):
csum+=phi**i+1
out[i]=csum*trend
return out
之前说过直接计算幂是不必要的,可以重写算法来完全避免这种情况。
第二种方法
@nb.njit()
def pow_cumsum_2(horizon,phi,trend):
out=np.empty(horizon)
out[0]=2.*trend
TMP=2.
val=phi
for i in range(horizon-1):
TMP=(val+1+TMP)
out[i+1]=TMP*trend
val*=phi
return out
计时
%timeit np.cumsum(phi ** np.arange(horizon) + 1) * trend
7.44 µs ± 89.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit pow_cumsum(horizon,phi,trend)
882 ns ± 4.91 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit pow_cumsum_2(horizon,phi,trend)
559 ns ± 3.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)