如何将矢量化函数应用于 numpy 数组的前一个元素？答案

【问题标题】：How can I apply a vectorized function to the previous element of a numpy array?如何将矢量化函数应用于 numpy 数组的前一个元素？
【发布时间】：2018-08-08 12:32:00
【问题描述】：

我想应用这样的功能：

s[i] = a*x[i] + (1 - a)*s[i-1]

其中s 和x 都是长度相同的数组。

我不想使用 for 循环，因为这些数组非常大（>50 百万）。我试过做这样的事情

def f(a,x):
    s = [0]*len(x)
    s[i] = a*x[i] + (1 - a)*s[i-1]
    return s

但当然 i 没有定义，所以这不起作用。

有没有办法使用map 或numpy.apply_along_axis 或其他一些矢量化方法来做到这一点？

我还没有遇到不使用 for 循环将函数应用于数组的当前和先前元素的方法，而这正是我想在这里了解的真正方法。

编辑

为了明确起见，这里是 for 循环实现，它有效，但我想避免

s = [0]*len(x)
a=0.45
for i in range(len(x)):
    s[i] = a*x[i] + (1-a)*s[i-1]

s[0] = x[0] # reset value of s[0]

【问题讨论】：

50k 并不是很大。顺便说一句，新的s[0] 会是什么？您的公式是否仅适用于 i>0？
是的，它只适用于 i>0
s[i-1] 应该是原始 s 的值，还是已经更新的先前条目？
您能否更改计算以使用cumsum 或其他ufunc accumulate 方法？

标签： python arrays numpy vectorization

【解决方案1】：

您可以避免循环，尽管矢量化更像是“为每个值计算所有内容”。

import numpy as np

# Loop
def fun(x, a):
    s = [0] * len(x)
    for i in range(len(x)):
        s[i] = a * x[i] + (1 - a) * s[i - 1]
    s[0] = x[0]
    return s

# Vectorized
def fun_vec(x, a):
    x = np.asarray(x, dtype=np.float32)
    n = np.arange(len(x))
    p = a * (1 - a) ** n
    # Trick from here: https://stackoverflow.com/q/49532575/1782792
    pz = np.concatenate((np.zeros(len(p) - 1, dtype=p.dtype), p))
    pp = np.lib.stride_tricks.as_strided(
        pz[len(p) - 1:], (len(p), len(p)),
        (p.strides[0], -p.strides[0]), writeable=False)
    t = x[np.newaxis] * pp
    s = np.sum(t, axis=1)
    s[0] = x[0]
    return s


x = list(range(1, 11))
a = 0.45
print(np.allclose(fun(x, a), fun_vec(x, a)))
# True

这种策略虽然占用 O(n²) 内存，而且计算量更大。根据具体情况，由于并行性，它可能会更快（我做了类似的事情来消除 TensorFlow 中的 tf.while_loop 并取得了巨大成功），但在这种情况下它实际上更慢：

x = list(range(1, 101))
a = 0.45
%timeit fun(x, a)
# 31 µs ± 85.2 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit fun_vec(x, a)
# 147 µs ± 2.42 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

所以，可以有一个非循环版本，但它更多的是一种好奇心。

【讨论】：

【解决方案2】：

正如我在answer to basically the same question 中所写，您不能：

除了显式的for 循环外，没有其他方法（通常）。这是因为没有办法在整个行（因为每一行都依赖于其他行）。

让这更难的是您可以轻松生成chaotic behavior，例如看似无辜的样子 logistic map:x_{n+1} = r * x_n * (1 - x_{n-1}).

只有当你设法找到一个封闭的形式，基本上消除了递归关系。但这必须为每个重复关系完成，我很确定你不是甚至保证存在封闭形式...

【讨论】：