与 NumPy “累积”类似的功能，但在多个数组上答案

【问题标题】：Similar functionality to NumPy "accumulate" but on multiple arrays与 NumPy “累积”类似的功能，但在多个数组上
【发布时间】：2020-07-27 13:50:19
【问题描述】：

我正在尝试找到一种模仿以下代码的快速方法（最好不使用for 循环）：

# a and b are arrays of the same length
# f is function of 3 scalar variables

res = np.zeros(len(a))
temp = 1

for i in range(len(a)):
    temp = f(temp, a[i], b[i])
    res[i] = temp

请注意，temp 变量会随着每次迭代而更新。这段代码与 NumPy 的累加机制可以实现的功能非常相似。但是，在我的例子中，函数f 中使用了几个数组。

函数f 可以是这样的

def f(temp, a_scalar, b_scalar):
    return temp * a_scalar + b_scalar

【问题讨论】：

显示你的 f 函数
对于任意 python 函数 f，没有真正的方法可以让它快速，因为 python 必须在每次迭代时运行解释函数 f。 numpy.ufuncs 可以快速累积，因为它们是用 C 实现的，因此 numpy 不必在 python 中进行循环。如果您想加快某些特定功能f，则可能有一种方法可以使用numpy 功能的组合来实现。如果没有，并且您仍然确实需要更快的代码，请考虑用 C 语言编写它并使用 ctypes 调用它。
@mrip。使用 C 或 numba 的另一个原因是，即使对于这个微不足道的例子，数值稳定性也可能是一个问题。几个大数的累积乘积确实会导致一些扩展问题。

标签： python arrays numpy accumulate

【解决方案1】：

与大多数矢量化问题一样，您必须根据每个功能制定解决方案。对于您提供的示例，请注意 res[i] = res[i - 1] * a[i] + b[i] 可以重写为以下序列：

res[0] = a[0] + b[0]
res[1] = (a[0] + b[0]) * a[1] + b[1] = a[0] * a[1] + b[0] * a[1] + b[1]
res[2] = a[0] * a[1] * a[2] + b[0] * a[1] * a[2] + b[1] * a[2] + b[2]
...
res[n] = a[0] * ... * a[n] + b[0] * (a[1] * ... an]) + b[1] * (a[2] * ... * a[n]) + ... b[n-1] * a[n] + b[n]

这些项是cumprod(a[:n:-1]) 和b 的乘积的累积和。您也可以将a 的累积积的尾部视为cumprod(a) / cumprod(a[:n]。所以我们有

res[n] = prod(a[:n]) + sum(prod(a[:n]) * b[:n] / cumprod(a[:n]))

由于这个关系由乘积和商组成，我们可以把它分成

res[n] = cumprod(a)[n] * (1 + sum(b[:n] / cumprod(a[:n])))

剩余的递归项也可以重写为累积积的累积和：

res[n] = cumprod(a)[n] * (1 + cumsum(b / cumprod(a))[n])

递归现在已经被分成了numpy支持的一系列累加。

警告

请记住，与原始公式相比，分离在数值上的稳定性要差得多。如果这很重要，请使用原始递归并使用 numba 编译，或者直接用 C 编写代码。

示例

np.random.seed(0xBEEF)
a = np.random.uniform(-2, 2, size=10)
b = np.random.uniform(-3, 3, size=10)
res = np.cumprod(a) * (1 + np.cumsum(b / np.cumprod(a)))

问题中的原始实现和此处显示的代码产生几乎相同的结果：

array([-0.04044657,  2.5090314 ,  0.66147071, -1.40795856,  0.93920848,
       -1.75817347, -4.45795915,  2.2258663 ,  2.38540097,  1.43741821])

差异微乎其微，最多为 1e-16 量级，这是由于累积与增量更新导致精度损失所致。

【讨论】：