Moving average function in C++ for use with R

Question

I am trying to improve the code and speed up in C++ (Rcpp) a (centered) weighted moving average function I coded.

An example of what the roll_mean function does. Note that the function works no matter what the size of x is and adapts to both tails of my data

w=c(1/2,1,1/2)
x=c(4,2,6,12)
res=c(2,5,7,3) 
res=c((x[1:2]*w[2:3])/sum(w[2:3]),x[1:3]*w[1:3]/sum(w[1:3]),x[2:4]*w[1:3]/sum(w[1:3]),x[3:4]*w[1:2]/sum(w[1:2]))

The file PartialMA.cpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector roll_mean(const NumericVector& x,
                        const NumericVector& w) {

  int n = x.size();
  int w_size = w.size();
  int size = (w_size - 1) / 2;

  NumericVector res(n);
  int i, ind_x, ind_w;

  double tmp_wsum, tmp_xwsum, tmp_w;

  for (i = 0; i < n; i++) {
    tmp_xwsum = 0;
    tmp_wsum = 0;
    for (ind_x = i - size, ind_w = 0; ind_x < i + size; ind_x++, ind_w++) {
      if((ind_x >= 0) & (ind_x < n)){
      tmp_w = w(ind_w);
      tmp_xwsum += x(ind_x) *  tmp_w;
      tmp_wsum += tmp_w;
      }
    }
    res[i] = tmp_xwsum / tmp_wsum;
  }

  return res;
}

I tried to replace the loop + if statement with this to minimize the number of iteration:

for (ind_x = std::max(0, i - size), ind_w = std::max(0, size-1); ind_x < std::min(n, i + size); ind_x++, ind_w++) {

I feel like I am not rigorous enough and I would be very grateful if someone could help me improving my code and eventually speed up the function as much as possible.

Answer 1

Overall your code is not bad, but we can do better. Let's start with the coding style.

Your variables are clearly named and your code is easy to follow, however the use of & is not necessary.

const NumericVector &x

R objects are always passed by reference in Rcpp even without the &. Dirk has some slides (see 29,30) with more information about this topic here. It is also in the Rcpp FAQ. The point is, is that we can safely drop these ampersands.

Now, let's talk about how we can improve the performance of your code. We first note that there are many unnecessary checks in the inner for loop.

if((ind_x >= 0) & (ind_x < n))

This can be avoided outside of the loop like so :

// [[Rcpp::export]]
NumericVector roll_meanReduceChecks(const NumericVector x,
                                    const NumericVector w) {
    
    int n = x.size();
    int w_size = w.size();
    int size = (w_size - 1) / 2;
    
    NumericVector res(n);
    int i, ind_x, ind_w, strt, endx;
    
    double tmp_wsum, tmp_xwsum, tmp_w;
    
    for (i = 0; i < n; i++) {
        tmp_xwsum = 0;
        tmp_wsum = 0;
        
        if ((i - size) <= 0) {
            strt = 0;
            ind_w = size - i;
        } else {
            strt = i - size;
            ind_w = 0;
        }
    
        endx = ((i + size) >= n) ? n : (i + size);
     
        for (ind_x = strt; ind_x < endx; ind_x++, ind_w++) {
            // This check is no longer necessary
            // if((ind_x >= 0) & (ind_x < n)){ 
            tmp_w = w(ind_w);
            tmp_xwsum += x(ind_x) * tmp_w;
            tmp_wsum += tmp_w;
        }
        
        res[i] = tmp_xwsum / tmp_wsum;
    }
    
    return res;
}

With this modification, we get about 10% faster timings:

set.seed(42)
x <- sample(10^3, 10^4, TRUE)
w <- sample(100, 10^3 , TRUE) / 100

## Gives the same result
all.equal(roll_mean(x, w), roll_meanReduceChecks(x, w))
[1] TRUE

library(microbenchmark)
microbenchmark(roll_mean(x, w), roll_meanReduceChecks(x, w), 
               times = 50, unit = "relative")
    Unit: relative
                       expr      min       lq     mean   median       uq      max neval
            roll_mean(x, w) 1.138508 1.135057 1.119695 1.132882 1.123308 1.149893    50
roll_meanReduceChecks(x, w) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    50

We can still do better. The real savings comes into play with the use of iterators. There is a lot of good information out there regarding this topic and here is just one resource written by Hadley Wickham that is a good starting place.

Putting it all together we have:

// [[Rcpp::export]]
NumericVector roll_meanIterator(const NumericVector x,
                                const NumericVector w) {
    
    int n = x.size();
    int w_size = w.size();
    int size = (w_size - 1) / 2;
    
    NumericVector res(n);
    NumericVector::const_iterator itw, itx, itEnd;
    double tmp_wsum, tmp_xwsum;
    
    for (int i = 0; i < n; i++) {
        tmp_xwsum = 0;
        tmp_wsum = 0;
        
        if ((i - size) <= 0) {
            itx = x.begin();
            itw = w.begin() + size - i;
        } else {
            itx = x.begin() + i - size;
            itw = w.begin();
        }
        
        itEnd = ((i + size) >= n) ? x.end() : x.begin() + i + size;
        
        for (; itx < itEnd; itx++, itw++) {
            tmp_xwsum += (*itx) * (*itw);
            tmp_wsum += (*itw);
        }
        
        res[i] = tmp_xwsum / tmp_wsum;
    }
    
    return res;
}

Here is a sanity check:

all.equal(roll_mean(x, w), roll_meanIterator(x, w))
[1] TRUE

And here are the benchmarks:

microbenchmark(roll_mean(x, w), roll_meanIterator(x, w), 
                          times = 50, unit = "relative")
Unit: relative
                   expr      min       lq     mean   median       uq      max neval
        roll_mean(x, w) 8.480686 8.458674 8.493975 8.894612 8.548443 8.453865    50
roll_meanIterator(x, w) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    50

That is nearly 9x faster. Not bad, considering we didn't have to change the code that much.

Update

Below is an example of over optimizing and a lesson on "when to stop" as the code appears to be more efficient but is much messier than the code above with little to no trade-off in efficiency.

If we examine tmp_wsum more closely, we will find that on every iteration, it is only changing at the endpoints. In roll_meanReduceChecks above, if we insert std::cout << ind_w << ' '; just before the inner for loop and std::cout << ind_w << std::endl; just after it, here is the output for the first 5 iterations:

499 998
498 998
497 998
496 998
495 998

And the last 5 iterations:

0 504
0 503
0 502
0 501
0 500

This seems like an opportunity for even further optimizing our code as we can alter our code such that we won't have to reconstruct tmp_wsum from scratch every time. However, in order to implement this, we must add two extra variable, and two extra checks to every iteration. Observe:

NumericVector roll_meanOverOptimized(const NumericVector x,
                                     const NumericVector w) {
    
    int n = x.size();
    int w_size = w.size();
    int size = (w_size - 1) / 2;
    
    NumericVector res(n);    // must add variables for checks below
    NumericVector::const_iterator itw, itwBeg, itwEnd, itx, itxEnd;
    double tmp_wsum = 0, tmp_xwsum = 0;
    
    // We must first initialize tmp_wsum, itwBeg,
    // itwEnd, as well as populate res[0]
    itx = x.begin();
    itwBeg = itw = w.begin() + size;
    itxEnd = (size >= n) ? x.end() : x.begin() + size;
    
    for (; itx < itxEnd; itx++, itw++) {
        tmp_xwsum += (*itx) * (*itw);
        tmp_wsum += (*itw);
    }
    
    res[0] = tmp_xwsum / tmp_wsum;
    itwEnd = itw;
    
    // Start i @ 1 instead of 0 as the first
    // iteration was taken care of above
    for (int i = 1; i < n; i++) {
        tmp_xwsum = 0;
        
        if ((i - size) < 0) {
            itx = x.begin();
            itw = w.begin() + size - i;
        } else {
            itx = x.begin() + i - size;
            itw = w.begin();
        }
        
        // first check
        if (itw != itwBeg)
            tmp_wsum += (*itw);
        itwBeg = itw;

        itxEnd = ((i + size) > n) ? x.end() : x.begin() + i + size;
        
        // N.B. only one variable is being updated now
        for (; itx < itxEnd; itx++, itw++)
            tmp_xwsum += (*itx) * (*itw);
        
        // second check
        if (itw != itwEnd)
            tmp_wsum -= (*itw);
        itwEnd = itw;
        
        res[i] = tmp_xwsum / tmp_wsum;
    }
    
    return res;
}

Let's see if this extra optimization paid off:

## Gives the same results
all.equal(roll_meanIterator(x, w), roll_meanOverOptimized(x, w))
[1] TRUE

microbenchmark(roll_mean(x, w),
               roll_meanIterator(x, w),
               roll_meanOverOptimized(x, w), 
               times = 50, unit = "relative")
               Unit: relative
                        expr      min       lq     mean   median       uq      max neval
             roll_mean(x, w) 8.679318 8.751559 8.605352 8.796965 8.509649 7.683454    50
     roll_meanIterator(x, w) 1.015420 1.028507 1.041876 1.019177 1.017208 1.133470    50
roll_meanOverOptimized(x, w) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    50

Virtually identical results as roll_meanIterator. IMO, the ~2% efficiency gain is not worth it as we have made our code more difficult to read.

Update 2 (Many thanks to @HongOoi)

In the comments @HongOoi states, "the speedup isn't because of iterators per se, but because OP is using bounds-checked array access via the operator (). Switching to the regular [] should give you the same result as iterators. Let us test!!!

// [[Rcpp::export]]
NumericVector roll_meanBrackets(const NumericVector x,
                                const NumericVector w) {
    
    int n = x.size();
    int w_size = w.size();
    int size = (w_size - 1) / 2;
    
    NumericVector res(n);
    double tmp_wsum, tmp_xwsum;
    unsigned long int ind_x, ind_w, xEnd;
    
    for (int i = 0; i < n; i++) {
        tmp_xwsum = 0;
        tmp_wsum = 0;
        
        if ((i - size) < 0) {
            ind_x = 0;
            ind_w = size - i;
        } else {
            ind_x = i - size;
            ind_w = 0;
        }
        
        xEnd = ((i + size) > n) ? x.size() : i + size;
        
        for (; ind_x < xEnd; ind_x++, ind_w++) {
            tmp_xwsum += x[ind_x] * w[ind_w];
            tmp_wsum += w[ind_w];
        }
        
        res[i] = tmp_xwsum / tmp_wsum;
    }
    
    return res;
}

And here are the benchmarks:

microbenchmark(roll_mean(x, w),
               roll_meanIterator(x, w),
               roll_meanBrackets(x, w),
               times = 50, unit = "relative")
Unit: relative
                   expr      min       lq     mean   median        uq      max neval
        roll_mean(x, w) 8.596342 8.425527 8.374452 8.414136 8.3626927 7.479161    50
roll_meanIterator(x, w) 1.012213 1.007219 1.007606 1.000564 0.9973372 1.177148    50
roll_meanBrackets(x, w) 1.000000 1.000000 1.000000 1.000000 1.0000000 1.000000    50

As predicted by @HongOoi, the efficiency is nearly identical. This begs the question

Which method should you use? Iterators or C-style indexing?

Fortunately for us, there is some great information out there addressing this very topic. In fact, I found the answer to the stackoverflow question Iterator Loop vs Index Loop extremely helpful. Here is the summary by @TemplateRex:

"if you really need the index (e.g. access the previous or next element, printing/logging the index inside the loop etc.) or you need a stride different than 1, then I would go for the explicitly indexed-loop, otherwise I'd go for the range-for loop.

For generic algorithms on generic containers I'd go for the explicit iterator loop unless the code contained no flow control inside the loop and needed stride 1, in which case I'd go for the STL for_each + a lambda."

Answer 2

Normally, the weighted average weights should add to 1, and in your sample it adds to 2

w=c(1/2,1,1/2)

sum(w)=2

maybe it should be?

w=c(1/4,1/2,1/4) 

and you can get a moving average with the filter function

F <- filter(x, filter = w, method = c("convolution"), sides = 2)

If your convolution kernel w is too large, and you want to check it for speed, I would try Fast Fourier convolution. FFT convolution should already be implemented somewhere in all languages.

convolve(x, w, conj = TRUE, type = c( "open"))

The FFT has the property that

#pseudocode
FFT(F) = FFT(x) * FFT(w)

so, to get F you do

#pseudocode
F <- inverseFFT( FFT(x) * FFT(w) )

Some disadvantages with the FFT are that

• commonly it needs length(x) == 2^n (n integer)
• it is periodic.

Those problems tend to be solved by padding the data x with 0

Also, in your Rcpp code, frequently a weighted average kernel (w) is symmetrical and has repeated weights, so you can take advantage of it by saving the multiplications to avoid repeating them. I do not do c++ code, but somebody else may use it to improve your code.