python – 迭代使用自己的输出的数组的最佳方法

for n in range(1,len(X)):
    Y[n] = X[n] + X[n - 1] + Y[n-1]

访问单个NumPy数组元素或(元素 – )迭代NumPy数组很慢(就像真的很慢).如果你想对NumPy数组进行手动迭代：就是不要这样做！

但你有一些选择.最简单的方法是将数组转换为Python列表并迭代列表(听起来很愚蠢,但请留在我身边 – 我将在答案1的末尾提供一些基准测试)：

X = X.tolist()
Y = Y.tolist()
for n in range(1,len(X)):
    Y[n] = X[n] + X[n - 1] + Y[n-1]

如果您还对列表使用直接迭代,则可能更快：

X = X.tolist()
Y = Y.tolist()
for idx,(Y_n_m1,X_n,X_n_m1) in enumerate(zip(Y,X[1:],X),1):
    Y[idx] = X_n + X_n_m1 + Y_n_m1

然后有更复杂的选项需要额外的包.最值得注意的是Cython和Numba,它们被设计为直接处理数组元素并尽可能避免Python开销.例如,使用Numba,您可以在jitted(即时编译)函数中使用您的方法：

import numba as nb

@nb.njit
def func(X,Y):
    for n in range(1,len(X)):
        Y[n] = X[n] + X[n - 1] + Y[n-1]

X和Y可以是NumPy数组,但numba将直接在缓冲区上工作,超过其他方法(可能是数量级).

Numba比Cython更“重”,但它可以更快更容易使用.但没有conda,很难安装numba … YMMV

然而,这里也是一个Cython版本的代码(使用IPython魔法编译,如果你不使用IPython,它会有点不同)：

In [1]: %load_ext cython

In [2]: %%cython
   ...:
   ...: cimport cython
   ...:
   ...: @cython.boundscheck(False)
   ...: @cython.wraparound(False)
   ...: cpdef cython_indexing(double[:] X,double[:] Y):
   ...:     cdef Py_ssize_t n
   ...:     for n in range(1,len(X)):
   ...:         Y[n] = X[n] + X[n - 1] + Y[n-1]
   ...:     return Y

仅举一个例子(基于the timing framework from my answer to another question),关于时间：

import numpy as np
import numba as nb
import scipy.signal

def numpy_indexing(X,len(X)):
        Y[n] = X[n] + X[n - 1] + Y[n-1]
    return Y

def list_indexing(X,Y):
    X = X.tolist()
    Y = Y.tolist()
    for n in range(1,len(X)):
        Y[n] = X[n] + X[n - 1] + Y[n-1]
    return Y

def list_direct(X,Y):
    X = X.tolist()
    Y = Y.tolist()
    for idx,1):
        Y[idx] = X_n + X_n_m1 + Y_n_m1
    return Y

@nb.njit
def numba_indexing(X,len(X)):
        Y[n] = X[n] + X[n - 1] + Y[n-1]
    return Y


def numpy_cumsum(X,Y):
    Y[1:] = X[1:] + X[:-1]
    np.cumsum(Y,out=Y)
    return Y

def scipy_lfilter(X,Y):
    a = [1,-1]
    b = [1,1]
    return Y[0] - X[0] + scipy.signal.lfilter(b,a,X)

# Make sure the approaches give the same result
X = np.random.random(10000)
Y = np.zeros(10000)
Y[0] = np.random.random()

np.testing.assert_array_equal(numba_indexing(X,Y),numpy_indexing(X,Y))
np.testing.assert_array_equal(numba_indexing(X,numpy_cumsum(X,Y))
np.testing.assert_almost_equal(numba_indexing(X,scipy_lfilter(X,cython_indexing(X,Y))

# Timing setup
timings = {numpy_indexing: [],list_indexing: [],list_direct: [],numba_indexing: [],numpy_cumsum: [],scipy_lfilter: [],cython_indexing: []}
sizes = [2**i for i in range(1,20,2)]

# Timing
for size in sizes:
    X = np.random.random(size=size)
    Y = np.zeros(size)
    Y[0] = np.random.random()
    for func in timings:
        res = %timeit -o func(X,Y)
        timings[func].append(res)

# Plottig absolute times

%matplotlib notebook
import matplotlib.pyplot as plt

fig = plt.figure(1)
ax = plt.subplot(111)

for func in timings:
    ax.plot(sizes,[time.best for time in timings[func]],label=str(func.__name__))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('size')
ax.set_ylabel('time [seconds]')
ax.grid(which='both')
ax.legend()
plt.tight_layout()

# Plotting relative times

fig = plt.figure(1)
ax = plt.subplot(111)

baseline = numba_indexing # choose one function as baseline
for func in timings:
    ax.plot(sizes,[time.best / ref.best for time,ref in zip(timings[func],timings[baseline])],label=str(func.__name__))
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_xlabel('size')
ax.set_ylabel('time relative to "{}"'.format(baseline.__name__))
ax.grid(which='both')
ax.legend()

plt.tight_layout()

结果如下：

绝对运行时

相对运行时间(与numba函数相比)

因此,只需将其转换为列表,您将大约快3倍！通过直接在这些列表上进行迭代,您可以获得另一个(更小的)加速,在此基准测试中只有20％,但与原始解决方案相比,我们现在快了近4倍.使用numba,与列表操作相比,您可以将速度提高100倍以上！ Cython只比numba慢一点(约40-50％),可能是因为我没有用Cython来完成所有可能的优化(通常它的速度不会超过10-20％).但是对于大型阵列,差异会变小.

1我在another answer确实进入了更多的细节.那个Q A是关于转换到一个集合但是因为set使用(隐藏)“手动迭代”它也适用于此处.

我包括了NumPy cumsum和Scipy lfilter方法的时间安排.与numba函数相比,对于小阵列,这些大约慢20倍,对于大阵列大约慢4倍.但是,如果我正确地解释了这个问题,那么你不仅要寻找一些方法,而不仅仅是在例子中应并非每个自引用循环都可以使用NumPy或SciPys过滤器的cum *函数来实现.但即便如此,似乎他们无法与Cython和/或numba竞争.