Optimizing the exponential smoothing of a big array

Question

I have a large set of values (let's say 1M entries) where I need to apply an exponential smoothing algorithm, but only incrementing one value at a time (all others decay to zero). The trivial implementation would be (pseudo-code):

function smooth_once(index, x):
    // Decay everything
    for s in s_vec:
        s = s * (1 - alpha)
    // Only increase one value, for the other x=0
    s_vec[index] += alpha * x

function get_s(index):
    return s_vec[index]

This is rather slow, as every element of s_vec need to be scanned at every iteration, whereas only one need to be incremented. Also, I only need to check few values from s_vec from time to time (much less frequently than the frequency at which smooth_once is called).

Note: the smoothing and read sequence is unpredictable, and is usually interleaved, for example:

smooth_once(6573, 1.23)
get_s(8892)
smooth_once(3345, 2.45)
smooth_once(6874, 1.10)
get_s(3345)
get_s(1254)
...etc...

In order to speed up this, I was thinking of multiplying the x value by a factor (scale) which is multiplied every time, and dividing everything when the factor grows too large. Basically, this boils down to rescale everything by an exponentially-incrementing factor.

scale = 1
k = 1 / (1 - alpha)
function optimized_smooth(index, x):
    // Rescale everything
    scale = scale * k
    if scale > scale_max:
        // Switch back scale to 1
        for s in s_vec:
            s = s / scale
        scale = 1
    s_vec[index] = v + alpha * x * scale

function get_s(index):
    // Scale back for reading
    return s_vec[index] / scale

IMHO this would be much faster (assuming k is not too large), as you only need to scan s_vec when the scale overflow and you rescale the whole array.

• What do you think of this method?
• Is this already known and has a name?
• Apart from loosing a few bits of precision on the values (depending on beta_max), do you see any drawbacks?

---

Addendum. A small demo program shows that both algorithms end up with the same result on a controlled sequence of smoothing indexes.

Answer 1

After thinking twice about your solution, I am sure that your algorithm works well. IMHO it could be seen as a variant of lazy evaluation. Here is another variant.

Change smooth_once so it only puts the pairs (index,x) into a queue (lets call the arrays indexes and xs).

function smooth_once(index, x):
   indexes.append(index)
   xs.append(x)

Here is the method to calculcate the smoothing for the whole queue at once:

function nth_smooth(indexes, xs):
    n = size(indexes)  // equal to size(xs)
    if n==0: 
        return
    cn = (1-alpha) ** n // (1-alpha) to the power of n
    for s in s_vec:
        s = s * cn
    c=1
    for i = 0 to n-1:
       s_vec[indexes[i]] += alpha * xs[i] *c
       c = c * (1-alpha)

And here is the method to retrieve the result:

function get_s(index):
    nth_smooth(indexes,xs)
    indexes.clear()
    xs.clear()
    return s_vec[index]

I guess it is not too hard to understand why this works, it simply replaces the repeated multiplications of s_vec with 1-alpha by one multiplication with (1-alpha) ** n, and adds the "x*alpha" additions afterwards, taking the influence of 1-alpha to those pertubances into account.

However, if this will be faster or slower than your own proposal depends on the number of calls of smooth_once before the next get_s occurs, compared to the frequency of having the condition scale > scale_max fulfilled (in your variant).