高效的轧制窗产品的总和

Question

给定i = 0到N-1的数字序列a [i]，我试图计算以下总和：

a[0] * a[1] * a[2] +
a[1] * a[2] * a[3] +
a[2] * a[3] * a[4] +
...
a[N-4] * a[N-3] * a[N-2] +
a[N-3] * a[N-2] * a[N-1]

我想使乘法组的大小G（在上面的例子中为3）成为一个可变参数。然后，可以使用简单的O（N * G）算法天真地获得结果，该算法可以用伪代码编写，如下所示：

sum = 0
for i from 0 to (N-G-1):
  group_contribution = 1
  for j from 0 to (G-1):
    group_contribution *= a[i+j]
  sum += group_contribution

然而，对于大G，显然该算法非常低效，特别是假设序列a [i]的数量事先不知道并且必须在运行时昂贵地计算。

出于这个原因，我考虑使用以下复杂度O（N + G）算法，它通过计算轧制产品来回收序列a [i]的值：

sum = 0
rolling_product = 1
for i from 0 to (G-1):
  rolling_product *= a[i]
sum += rolling_product
for i from G to (N-1):
  rolling_product /= a[i-G]
  rolling_product *= a[i]
  sum += rolling_product

然而，我关注标准浮点表示中除法的数值稳定性。

我有兴趣知道是否有一种稳定，快速的方法来计算这个总和。这对我来说感觉像是一项基本的数字任务，但目前我不知道它是如何高效的。

谢谢你的任何想法！

Answer 1

是的，如果您仔细计算反向部分产品，则无需进行划分。

def window_products(seq, g):
    lst = list(seq)
    reverse_products = lst[:]
    for i in range(len(lst) - 2, -1, -1):
        if i % g != len(lst) % g:
            reverse_products[i] *= reverse_products[i + 1]
    product = 1
    for i in range(len(lst) - g + 1):
        yield reverse_products[i] * product
        if i % g == len(lst) % g:
            product = 1
        else:
            product *= lst[i + g]


print(list(window_products(range(10), 1)))
print(list(window_products(range(10), 2)))
print(list(window_products(range(10), 3)))

Answer 2

作为序言，您可以考虑在两种算法上运行一些测试用例并比较结果（例如，作为相对误差）。

接下来，如果你有额外的内存和时间，这里是O（N log2 G）时间和内存的稳定方法。它类似于range minimum query问题的恒定时间，线性空间的方法。

预计算两种功率范围的产品

设B [i] [j]是从位置i开始的2j个元素的乘积，所以

B [i] [j] = a [i] x a [i + 1] x ... x a [i + 2j-1]

我们感兴趣的是B中的N log2 G值，即0≤j≤log2G的值。我们可以在O（1）中计算这些值中的每一个，因为

B [i] [j] = B [i] [j-1] x B [i + 2j-1] [j-1]

计算总和

为了计算总和中的一个项，我们将G分解为两个幂大小的块。例如，如果G = 13，那么第一项是

a [0]×...×a [12] =（a [0]×...×a [7]）×（a [8]×...×a [11]）×a [12] = B [0] [3]×B [8] [2]×B [12] [0]

每个O（N）项可以在O（log 2 G）时间内计算，因此找到和的总复杂度是O（N log 2 G）。

Answer 3

你的滚动产品是个好主意，但正如你所说，它存在稳定性问题。我会这样修复：

• 单独使用类似的系统来跟踪零和负数的数量。那些是整数和，所以没有稳定性问题。
• 而不是计算所有a[i]的滚动乘积，计算log(abs(a[i]))的滚动总和，不包括零。然后，当你需要产品时，它是(num_zeros > 0 ? 0.0 : exp(log_sum)) * sign。这将解决主要的不稳定因素。
• 当您从聪明的滚动log_sum算法生成输出时，您应该同时构建一个新的log_sum，您没有从中减去任何内容。当新总和中的元素数量达到G时，用该数字覆盖滚动的llog_sum并将其重置为零。这将消除长期累积的任何舍入误差。

Answer 4

创建一个新的序列b，其中b [0] = a [0] * a [1] * a [2] * ... a [G-1]等。

现在你有一个更简单的问题，在哪里计算你可以保持总和的b值的总和，每次你添加一个值你减去b [0]并添加新值并将它们全部向下滑动一个（使用循环缓冲区，所以没有移动）。 [典型滑动窗移动平均型代码]

保持最后一个G a []值的缓存并计算要添加到结尾的新值只是O（G）操作，并且您只计算一次[i]。