Efficient sum of expensive rolling window products

Question

Given a sequence of numbers a[i] for i = 0 to N-1, I'm trying to calculate the following sum:

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]

I'd like to make the size G of a multiplied group (in the example above 3) a variable parameter. Then, the result can be naively obtained using a simple O(N*G) algorithm, which can be written in pseudocode like so:

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

However, for large G's it is obvious that this algorithm is terribly inefficient, especially assuming that numbers of the sequence a[i] aren't known in advance and have to be expensively calculated at runtime.

For this reason, I've considered using the following algorithm of complexity O(N+G) which recycles the values of the sequence a[i] by calculating a rolling product:

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

I'm however concerned about numerical stability of division in standard floating-point representation.

I'd be interested in knowing if there is a stable, faster way to calculate this sum. It feels like a basic numerical task to me, yet currently I'm at a loss how do it efficiently.

Thank you for any ideas!

Answer 1

Yes, if you compute the reverse partial products carefully, you don't need to divide.

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

As a preamble, you could consider running some test cases on both algorithms and compare the results (as a relative error, for example).

Next, if you have the additional memory and time, here is a stable approach in O(N log₂G) time and memory. It is similar to the approach in constant time, linearithmic space to the range minimum query problem.

Precomputing products of power-of-two ranges

Let B[i][j] be the product of 2^j elements of a starting at position i, so

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

We are interested in N log₂G values in B, namely those for 0 ≤ j ≤ log₂G. We can compute each of these values in O(1), since

B[i][j] = B[i][j - 1] × B[i + 2^{j - 1}][j - 1]

Computing the sum

To compute one term in the sum, we decompose G into power-of-two sized chunks. For example, if G = 13, then the first term is

a[0] × ... × a[12] = (a[0] × ... × a[7]) × (a[8] × ... × a[11]) × a[12] = B[0][3] × B[8][2] × B[12][0]

Each of the O(N) terms can be computed in O(log₂G) time, thus the total complexity of finding the sum is O(N log₂G).