Fast correlation in R using C and parallelization

Question

My project for today was to write a fast correlation routine in R using the basic skillset I have. I have to find the correlation between almost 400 variables each having almost a million observations (i.e. a matrix of size p=1MM rows & n=400 cols).

R's native correlation function takes almost 2 mins for 1MM rows and 200 observations per variable. I have not run for 400 observations per column, but my guess is it will take almost 8 mins. I have less than 30 secs to finish it.

Hence, I want to do do things.

1 - write a simple correlation function in C and apply it in blocks parallely (see below).

2 - The blocks - split the correlation matrix in three blocks (top left square of size K*K, bottom right square of size (p-K)(p-K), and top right rectangular matrix of size K(p-K)). This covers all cells in the correlation matrix corr since I only need the upper triangle.

3 - run the C function via a .C call parallely using snowfall.

n = 100
p = 10
X = matrix(rnorm(n*p), nrow=n, ncol=p)
corr = matrix(0, nrow=p, ncol=p)

# calculation of column-wise mean and sd to pass to corr function
mu = colMeans(X)
sd = sapply(1:dim(X)[2], function(x) sd(X[,x]))

# setting up submatrix row and column ranges
K = as.integer(p/2)

RowRange = list()
ColRange = list()
RowRange[[1]] = c(0, K)
ColRange[[1]] = c(0, K)

RowRange[[2]] = c(0, K)
ColRange[[2]] = c(K, p+1)

RowRange[[3]] = c(K, p+1)
ColRange[[3]] = c(K, p+1)

# METHOD 1. NOT PARALLEL
########################
# function to calculate correlation on submatrices
BigCorr <- function(x){
  Rows = RowRange[[x]]
  Cols = ColRange[[x]]    
  return(.C("rCorrelationWrapper2", as.matrix(X), as.integer(dim(X)), 
            as.double(mu), as.double(sd), 
            as.integer(Rows), as.integer(Cols), 
            as.matrix(corr)))
}

res = list()
for(i in 1:3){
  res[[i]] = BigCorr(i)
}

# METHOD 2
########################
BigCorr <- function(x){
    Rows = RowRange[[x]]
    Cols = ColRange[[x]]    
    dyn.load("./rCorrelation.so")
    return(.C("rCorrelationWrapper2", as.matrix(X), as.integer(dim(X)), 
          as.double(mu), as.double(sd), 
          as.integer(Rows), as.integer(Cols), 
          as.matrix(corr)))
}

# parallelization setup
NUM_CPU = 4
library('snowfall')
sfSetMaxCPUs() # maximum cpu processing
sfInit(parallel=TRUE,cpus=NUM_CPU) # init parallel procs
sfExport("X", "RowRange", "ColRange", "sd", "mu", "corr")  
res = sfLapply(1:3, BigCorr)
sfStop()  

Here is my problem:

for method 1, it works, but not the way I want it to. I believed, that when I pass the corr matrix, I am passing an address and C would be making changes at source.

# Output of METHOD 1
> res[[1]][[7]]
      [,1]      [,2]        [,3]       [,4]         [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1 0.1040506 -0.01003125 0.23716384 -0.088246793    0    0    0    0     0
 [2,]    0 1.0000000 -0.09795989 0.11274508  0.025754150    0    0    0    0     0
 [3,]    0 0.0000000  1.00000000 0.09221441  0.052923520    0    0    0    0     0
 [4,]    0 0.0000000  0.00000000 1.00000000 -0.000449975    0    0    0    0     0
 [5,]    0 0.0000000  0.00000000 0.00000000  1.000000000    0    0    0    0     0
 [6,]    0 0.0000000  0.00000000 0.00000000  0.000000000    0    0    0    0     0
 [7,]    0 0.0000000  0.00000000 0.00000000  0.000000000    0    0    0    0     0
 [8,]    0 0.0000000  0.00000000 0.00000000  0.000000000    0    0    0    0     0
 [9,]    0 0.0000000  0.00000000 0.00000000  0.000000000    0    0    0    0     0
[10,]    0 0.0000000  0.00000000 0.00000000  0.000000000    0    0    0    0     0
> res[[2]][[7]]
      [,1] [,2] [,3] [,4] [,5]        [,6]        [,7]        [,8]       [,9]       [,10]
 [1,]    0    0    0    0    0 -0.02261175 -0.23398448 -0.02382690 -0.1447913 -0.09668318
 [2,]    0    0    0    0    0 -0.03439707  0.04580888  0.13229376  0.1354754 -0.03376527
 [3,]    0    0    0    0    0  0.10360907 -0.05490361 -0.01237932 -0.1657041  0.08123683
 [4,]    0    0    0    0    0  0.18259522 -0.23849323 -0.15928474  0.1648969 -0.05005328
 [5,]    0    0    0    0    0 -0.01012952 -0.03482429  0.14680301 -0.1112500  0.02801333
 [6,]    0    0    0    0    0  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
 [7,]    0    0    0    0    0  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
 [8,]    0    0    0    0    0  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
 [9,]    0    0    0    0    0  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
[10,]    0    0    0    0    0  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
> res[[3]][[7]]
      [,1] [,2] [,3] [,4] [,5] [,6]       [,7]        [,8]        [,9]       [,10]
 [1,]    0    0    0    0    0    0 0.00000000  0.00000000  0.00000000  0.00000000
 [2,]    0    0    0    0    0    0 0.00000000  0.00000000  0.00000000  0.00000000
 [3,]    0    0    0    0    0    0 0.00000000  0.00000000  0.00000000  0.00000000
 [4,]    0    0    0    0    0    0 0.00000000  0.00000000  0.00000000  0.00000000
 [5,]    0    0    0    0    0    0 0.00000000  0.00000000  0.00000000  0.00000000
 [6,]    0    0    0    0    0    1 0.03234195 -0.03488812 -0.18570151  0.14064640
 [7,]    0    0    0    0    0    0 1.00000000  0.03449697 -0.06765511 -0.15057244
 [8,]    0    0    0    0    0    0 0.00000000  1.00000000 -0.03426464  0.10030619
 [9,]    0    0    0    0    0    0 0.00000000  0.00000000  1.00000000 -0.08720512
[10,]    0    0    0    0    0    0 0.00000000  0.00000000  0.00000000  1.00000000

But the original corr matrix remains unchanged:

> corr
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    0    0    0    0    0    0     0
 [4,]    0    0    0    0    0    0    0    0    0     0
 [5,]    0    0    0    0    0    0    0    0    0     0
 [6,]    0    0    0    0    0    0    0    0    0     0
 [7,]    0    0    0    0    0    0    0    0    0     0
 [8,]    0    0    0    0    0    0    0    0    0     0
 [9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

Question #1: Is there any way to ensure that the C function changes values of corr at source? I can still merge these three to create an upper triangular correlation matrix, but I wanted to know if change at source is possible. Note: this does not help me accomplish fast correlation since I am merely running a loop.

Question #2: For METHOD 2, how do I load the shared object to each core for parallel jobs on each core at the init step (and not how I have done it)?

Question #3: What does this error mean? I need some pointers, and I would love to debug it myself.

Question #4: Is there a fast way of calculating correlation over matrices 1MM by 400, in less then 30 secs?

When I run METHOD 2, I get the following error:

R(6107) malloc: *** error for object 0x100664df8: incorrect checksum for freed object - object was probably modified after being freed.
*** set a breakpoint in malloc_error_break to debug
Error in unserialize(node$con) : error reading from connection

Attached below is my plain vanilla C code for correlation:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stddef.h>
#include <R.h> // to show errors in R


double calcMean (double *x, int n);
double calcStdev (double *x, double mu, int n);
double calcCov(double *x, double *y, int n, double xmu, double ymu);        

void rCorrelationWrapper2 ( double *X, int *dim, double *mu, double *sd, int *RowRange, int *ColRange, double *corr) {

    int i, j, n = dim[0], p = dim[1];
    int RowStart = RowRange[0], RowEnd = RowRange[1], ColStart = ColRange[0], ColEnd = ColRange[1];
    double xyCov;

    Rprintf("\n p: %d, %d <= row < %d, %d <= col < %d", p, RowStart, RowEnd, ColStart, ColEnd);

    if(RowStart==ColStart && RowEnd==ColEnd){
        for(i=RowStart; i<RowEnd; i++){
            for(j=i; j<ColEnd; j++){
                Rprintf("\n i: %d, j: %d, p: %d", i, j, p);
                xyCov = calcCov(X + i*n, X + j*n, n, mu[i], mu[j]);
                *(corr + j*p + i) = xyCov/(sd[i]*sd[j]);
            }
        }
    } else {
        for(i=RowStart; i<RowEnd; i++){
            for (j=ColStart; j<ColEnd; j++){
                xyCov = calcCov(X + i*n, X + j*n, n, mu[i], mu[j]);
                *(corr + j*p + i) = xyCov/(sd[i]*sd[j]);
            }
        }
    }
}


// function to calculate mean

double calcMean (double *x, int n){
    double s = 0;
    int i;
    for(i=0; i<n; i++){     
        s = s + *(x+i);
    }
    return(s/n);
}

// function to calculate standard devation

double calcStdev (double *x, double mu, int n){
    double t, sd = 0;
    int i;

    for (i=0; i<n; i++){
        t = *(x + i) - mu;
        sd = sd + t*t;
    }    
    return(sqrt(sd/(n-1)));
}


// function to calculate covariance

double calcCov(double *x, double *y, int n, double xmu, double ymu){
    double s = 0;
    int i;

    for(i=0; i<n; i++){
        s = s + (*(x+i)-xmu)*(*(y+i)-ymu);
    }
    return(s/(n-1));
}

Panto

@MartinMorgan - R's native cor function (based on the build I have) takes more time as I mentioned above. I am using Andrey's suggestion below and it is taking about 2 mins for 1MM by 400 vars. Will update. – Panto

Answer 1

Using a fast BLAS (via Revolution R or Goto BLAS) you can calculate all these correlations fast in R without writing any C code. On my first generation Intel i7 PC it takes 16 seconds:

n = 400;
m = 1e6;

# Generate data
mat = matrix(runif(m*n),n,m);
# Start timer
tic = proc.time();
# Center each variable
mat = mat - rowMeans(mat);
# Standardize each variable
mat = mat / sqrt(rowSums(mat^2));   
# Calculate correlations
cr = tcrossprod(mat);
# Stop timer
toc = proc.time();

# Show the results and the time
show(cr[1:4,1:4]);
show(toc-tic)

The R code above reports the following timing:

 user  system elapsed 
31.82    1.98   15.74 

I use this approach in my MatrixEQTL package.
http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/

More information about various BLAS options for R is available here:
http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/runit.html#large

Tears

Without building R using any of the optimized BLAS, it is taking about 2 mins on my machine (2.9Ghz i7). I will install R with optimized BLAS and let you know. – Panto

Yes, @user1971988, I'd be curious of the performance of this code for you with BLAS. – Tears

Also, it is a custom on this site to accept an answer if you like it. – Tears

I am trying to replicate your times after re0installing R from source using an optimized BLAS. Give me a couple of days and I will update my results and accept your answer. – Panto

What method does it use? – Brigitta

Shouldn't you have cr = tcrossprod(mat) / nrow(mat) (alternatively cr = tcrossprod(mat) / (nrow(mat) - 1) depending on the estimator)? – Mun

@user2280549, no, I should not. The code is correct as written and you can test it to confirm. – Tears

@AndreyShabalin I did. n = 1000; m = 10; set.seed(12345); mat = matrix(runif(m*n),n,m); corr1 <- cor(mat)[1:4, 1:4]; mat = mat - rowMeans(mat); mat = mat / sqrt(rowSums(mat^2)); corr2 = tcrossprod(mat); all(corr1 == corr2) You can test on your own that it's not equal to standard implementation. Even dimensions of your matrix are wrong. – Mun

I correlate rows, not columns. – Tears

cr = tcrossprod(mat) calculates pearson correlations or spearman? – Tessitura

@Apex, it calculates matrix product, which gives correlations if matrices are normalized to have rows with zero mean and unit sum of squares. – Tears

It will produce Pearson correlations, unless you replaces the values in each row with ranks first. Then it would be spearman correlations. – Tears

Answer 2

A few things.

First, if you are using the .C interface for external calls, then by default it makes copies of all of the arguments. That's why the object corr isn't being modified. If you want to avoid this then you have to set DUP=false in the .C call. However, in general using .C to modify existing R objects is not the preferred way to do things. Instead, you probably want to create a new array and allow the external call to fill it in, like this.

corr<-.C("rCorrelationWrapper2", as.double(X), as.integer(dim(X)), 
        as.double(mu), as.double(sd), 
        as.integer(Rows), as.integer(Cols), 
        result=double(p*q))$result
corr<-array(corr,c(p,q))

Second, as far as writing a fast correlation function, the first thing you should try is compiling R with an efficient BLAS implementation. This will not just make your correlation function faster, it will make all of your linear algebra faster. Good free candidates are ACML from AMD, or ATLAS. Either of those will be able to compute correlation matrices very quickly. The speedup is more than just parallelization -- these libraries also are smart about things cache usage and are optimized at the assembly level so even with just one core you'll see a big improvement. http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acml/ http://math-atlas.sourceforge.net/

Finally, if you really want to write your own C code, I would suggest using openMP to automatically split up the computation among different threads, rather than doing it by hand. But, for something as basic as matrix multiplication, it's probably better to go with an available optimized library.

Amenity