Rcpp versus C - Outlier

Viewed 93

I have a general question with a specific example. I wrote a function to calculate columnwise variance of a matrix in C (using .Call interface) and C++ (using Rcpp interface). Looking at the following benchmarks I wonder:

> microbenchmark(times = 1000,
+                colVar(AB), # .Call Interface
+                colV(AB, ncol(AB), nrow(AB)), #Rcpp
+                apply(AB, 2, var)) #R
Unit: milliseconds
                         expr       min        lq      mean    median        uq        max neval
                   colVar(AB)  3.245000  3.350793  3.474891  3.433126  3.543796   5.110652  1000
 colV(AB, ncol(AB), nrow(AB))  4.064942  4.408336 10.215952  5.934169  6.383477  99.651530  1000
            apply(AB, 2, var) 28.260730 30.740058 46.674155 31.464449 33.586160 129.343892  1000
> 

In distribution and mean the C and the C++ function perform pretty similar, however when it comes to a maximum value there is a huge difference. Can anybody explain to me why? This is especially interesting since I am trying to learn C/C++ but also because I want to write more complicated functions in C/C++, where this could actually matter. AB is a matrix with dimension 1000 x 1000, created with 1 000 000 rnorm() values. Below you find the Codes for my C and Rcpp functions:

C (R-Level):

colVar <- function(x){
  .Call("colV", x, ncol(x), nrow(x))
}

C (C-Level):

#include <R.h>
#include <Rinternals.h>
#include <math.h>


SEXP colV(SEXP y, SEXP n, SEXP r){
    int *nc = INTEGER(n);
    double *x = REAL(y);
    int d = length(y);
    int *nr = INTEGER(r);
    int i, j, z;
    //int d = nr * nc;

    double xSq[(d)];
    SEXP result;
    PROTECT(result = allocVector(REALSXP, (*nc)));
    memset(REAL(result), 0, (*nc) * sizeof(double));
    double *colVar = REAL(result);
    int fr = ((*nr) - 1);


    for(z = 0; z < (d); z++){
        xSq[z] = pow(x[z], 2);
    }

    for(i = 0; i < (*nc); i++){
        double colMean = 0;
        double xSm = 0;
        double colMsq = 0;
        for(j = 0; j < (*nr); j++){
            colMean += ((x[(j + ((*nr) * i)) ]) / (*nr));
            xSm += (xSq[(j + (*nr * i))]);
        }
        colMsq = (*nr) * (pow(colMean, 2));
        colVar[i] = ((xSm - colMsq) / fr);
    }
    UNPROTECT(1);
    return(result);
}

And the Rcpp-Function:

cppFunction(plugins = "unwindProtect",'NumericVector colV(NumericVector y, int n, int r){
            int nc = n;
            NumericVector x = y;
            int nr = r;
            int d = n * r;
            int i, j, z;

            // NumericVector colMean (nc);
            NumericVector xSq (d);
            // NumericVector colMsq (nc);
            // NumericVector xSm (nc);

            NumericVector colVar (nc);

            int fr = ((nr) - 1);


            for(z = 0; z < (d); z++){
               xSq[z] = x[z] * x[z];
            }

            for(i = 0; i < (nc); i++){
                double colMean = 0;
                double xSm = 0;
                double colMsq = 0;
                for(j = 0; j < (nr); j++){
                    colMean += ((x[(j + ((nr) * i)) ]) / (nr));
                    xSm += (xSq[(j + (nr * i))]);
                }
                colMsq = (nr) * (colMean * colMean);
                colVar[i] = ((xSm - colMsq) / fr);
            }
            return colVar;
            }')

I have commented out stuff in the C++ function to make it as similar as possible to the C function. If anybody of you can help me with my question I would be very thankful.

0 Answers
Related