这个循环处理非常缓慢，我希望能加速它。但是，由于一个值的结果依赖于前一个值的结果，我不知道如何将其向量化。有什么建议吗？
import numpy as np
sig = np.random.randn(44100)
alpha = .9887
beta = .999
out = np.zeros_like(sig)
for n in range(1, len(sig)):
if np.abs(sig[n]) >= out[n-1]:
out[n] = alpha * out[n-1] + (1 - alpha) * np.abs( sig[n] )
else:
out[n] = beta * out[n-1]
Numba的即时编译器应该很好地处理你面临的索引开销，通过在第一次执行期间将函数编译为本机代码：
In [1]: %cpaste
Pasting code; enter '--' alone on the line to stop or use Ctrl-D.
:import numpy as np
:
:sig = np.random.randn(44100)
:alpha = .9887
:beta = .999
:
:def nonvectorized(sig):
: out = np.zeros_like(sig)
:
: for n in range(1, len(sig)):
: if np.abs(sig[n]) >= out[n-1]:
: out[n] = alpha * out[n-1] + (1 - alpha) * np.abs( sig[n] )
: else:
: out[n] = beta * out[n-1]
: return out
:--
In [2]: nonvectorized(sig)
Out[2]:
array([ 0. , 0.01862503, 0.04124917, ..., 1.2979579 ,
1.304247 , 1.30294275])
In [3]: %timeit nonvectorized(sig)
10 loops, best of 3: 80.2 ms per loop
In [4]: from numba import jit
In [5]: vectorized = jit(nonvectorized)
In [6]: np.allclose(vectorized(sig), nonvectorized(sig))
Out[6]: True
In [7]: %timeit vectorized(sig)
1000 loops, best of 3: 249 µs per loop
编辑：根据评论建议，添加jit基准测试。 jit（非矢量化） 创建了一个轻量级的包装器，因此是一项廉价操作。
In [8]: %timeit jit(nonvectorized)
10000 loops, best of 3: 45.3 µs per loop
该函数本身在第一次执行时进行编译（因此是即时编译），这需要一些时间，但可能不会太长：
In [9]: %timeit jit(nonvectorized)(sig)
10 loops, best of 3: 169 ms per loop
一旦分析了依赖关系，大部分“向量化”并行性就被排除在外了。（即使JIT编译器也无法“反对”这种依赖障碍进行向量化）
您可以以向量化方式预先计算一些重复使用的值，但是没有直接的Python语法方式（没有外部JIT编译器的解决方法），将前向移位依赖循环计算安排到CPU向量寄存器对齐的共同并行计算中：
from zmq import Stopwatch # ok to use pyzmq 2.11 for [usec] .Stopwatch()
aStopWATCH = Stopwatch() # a performance measurement .Stopwatch() instance
sig = np.abs(sig) # self-destructive calc/assign avoids memalloc-OPs
aConst = ( 1 - alpha ) # avoids many repetitive SUB(s) in the loop
for thisPtr in range( 1, len( sig ) ): # FORWARD-SHIFTING-DEPENDENCE LOOP:
prevPtr = thisPtr - 1 # prevPtr->"previous" TimeSlice in out[] ( re-used 2 x len(sig) times )
if sig[thisPtr] < out[prevPtr]: # 1st re-use
out[thisPtr] = out[prevPtr] * beta # 2nd
else:
out[thisPtr] = out[prevPtr] * alpha + ( aConst * sig[thisPtr] ) # 2nd
一个很好的向量化加速示例可以在这种情况下看到，即计算策略可以沿着本地numpy数组的1D、2D甚至3D结构并行/广播。要获得大约100倍的加速效果，请参见RGBA-2D矩阵加速处理的向量化代码PNG图片处理（OpenGL着色器管道）。
即使这个简单的python代码修订已经将速度提高了大约2.8倍（现在，即没有进行安装以允许使用特定的JIT优化编译器）：
>>> def aForwardShiftingDependenceLOOP(): # proposed code-revision
... aStopWATCH.start() # ||||||||||||||||||.start
... for thisPtr in range( 1, len( sig ) ):
... # |vvvvvvv|------------# FORWARD-SHIFTING-LOOP DEPENDENCE
... prevPtr = thisPtr - 1 #|vvvvvvv|--STEP-SHIFTING avoids Numpy syntax
... if ( sig[ thisPtr] < out[prevPtr] ):
... out[ thisPtr] = out[prevPtr] * beta
... else:
... out[ thisPtr] = out[prevPtr] * alpha + ( aConst * sig[thisPtr] )
... usec = aStopWATCH.stop() # ||||||||||||||||||.stop
... print usec, " [usec]"
>>> aForwardShiftingDependenceLOOP()
57593 [usec]
57879 [usec]
58085 [usec]
>>> def anOriginalForLOOP():
... aStopWATCH.start()
... for n in range( 1, len( sig ) ):
... if ( np.abs( sig[n] ) >= out[n-1] ):
... out[n] = out[n-1] * alpha + ( 1 - alpha ) * np.abs( sig[n] )
... else:
... out[n] = out[n-1] * beta
... usec = aStopWATCH.stop()
... print usec, " [usec]"
>>> anOriginalForLOOP()
164907 [usec]
165674 [usec]
165154 [usec]
%timeit vectorized = jit(nonvectorized)的结果吗？谢谢。 - user3666197