Moving average function in C++ for use with R

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.

• Not a Rcpp expert so I can't help with a formal review... But if you were interested in a base R implementation, you could use the filter function. roll_mean <- function(x, w) as.vector(filter(x, w)) / sum(w) will get you almost there. You'll just have to deal with the end parts. Mar 24, 2018 at 0:26
• filter, rollapply and most R built-in function are about 10000 slower than this for my 10^7 sized vectors Mar 24, 2018 at 0:29
• filter is implemented in C so that's a bit surprising. I was expecting a bit of a slowdown compared to Rcpp but not that much... Rcpp it is, then! I hope you get help with your review. Mar 24, 2018 at 0:31
• @flodel: stats::filter() converts its input to ts, has multiple sanity checks for the various types of filters it supports, and loops over columns of the mts object before it gets to the C code. That results in a lot of intermediate memory allocations, which (I'm sure you know) can slow things down considerably for larger objects. The speed gains aren't from Rcpp, per-se. They're from omitting all the conversions and checks. Apr 1, 2018 at 15:04
• Thanks Joshua. Nonetheless, filter will be orders of magnitude faster than rollapply so the OP's "10,000x slower" comment can't apply to both (and I seriously doubt it applies to filter).. I wish the accepted answer's timings had included a filter based solution to see where it stands. As one has to weigh-in speed improvements versus the added complexity of writing and maintaining Rcpp code. Apr 2, 2018 at 23:01

2 Answers

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."

• You swapped the ordering in your benchmark, which requires a second look ;) Mar 26, 2018 at 12:06
• @miscco, good point, I will update now Mar 26, 2018 at 12:07
• @MaxFt, Sorry for responding so late. I will say that you should try to use iterators when you can, however indexing still has its value. If you find yourself writing extra code just to use iterators in place of a very simple indexing solution, I think you've found your answer. Mar 26, 2018 at 23:23
• Note that 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. Apr 10, 2018 at 23:30
• @HongOoi, I'm glad you addressed this as I was unsure of the difference. I will add another example with []. Feel free to edit my answer with any guidance you may have regarding this topic. Apr 11, 2018 at 17:05

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.

• Although you don't directly address the code, I think this still counts as an answer. I would recommend putting more information about FFT instead of just the link; links can rot. Apr 10, 2018 at 15:49