Corrélation rapide dans R avec C et parallélisation

Mon projet d’aujourd’hui était d’écrire une routine de corrélation rapide en R en utilisant les compétences de base que j’ai. Je dois trouver la corrélation entre près de 400 variables ayant chacune près d’un million d’observations (c’est-à-dire une masortingce de taille p = 1 MM lignes & n = 400 colonnes).

La fonction de corrélation native de R prend presque 2 minutes pour des lignes de 1MM et 200 observations par variable. Je n’ai pas couru pour 400 observations par colonne, mais je suppose que cela prendra presque 8 minutes. J’ai moins de 30 secondes pour le finir.

Par conséquent, je veux faire des choses.

1 – écrivez une fonction de corrélation simple en C et appliquez-la en blocs parallèlement (voir ci-dessous).

2 – Les blocs – divisez la masortingce de corrélation en trois blocs (carré haut gauche de taille K * K, carré bas droite de taille (pK) (pK) et masortingce rectangular haut droite de taille K (pK)). Cela couvre toutes les cellules de la masortingce de corrélation car je n’ai besoin que du sortingangle supérieur.

3 – exécutez la fonction C via un appel .C en parallèle en utilisant des chutes de neige.

 n = 100 p = 10 X = masortingx(rnorm(n*p), nrow=n, ncol=p) corr = masortingx(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 submasortingx 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 submasortingces 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() 

Voici mon problème:

pour la méthode 1, cela fonctionne, mais pas comme je le veux. Je croyais que lorsque je passais la masortingce corr, je passais une adresse et que C apporterait des changements à la 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 

Mais la masortingce corr origine rest inchangée:

 > 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 n ° 1: Y a-t-il un moyen de s’assurer que la fonction C change les valeurs de corr à la source? Je peux encore fusionner ces trois pour créer une masortingce de corrélation sortingangular supérieure, mais je voulais savoir si un changement à la source est possible. Remarque: cela ne m’aide pas à réaliser une corrélation rapide car je ne fais que lancer une boucle.

Question n ° 2: pour la MÉTHODE 2, comment puis-je charger l’object partagé sur chaque cœur pour des tâches parallèles sur chaque cœur à l’étape d’init (et non comment je l’ai fait)?

Question n ° 3: que signifie cette erreur? J’ai besoin de quelques pointeurs et j’aimerais pouvoir le déboguer moi-même.

Question n ° 4: Existe-t-il un moyen rapide de calculer la corrélation sur des masortingces de 1 MM par 400, en moins de 30 secondes?

Lorsque j’exécute METHOD 2, le message d’erreur suivant s’affiche:

 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 

Ci-dessous, mon code C simple et vanillé pour la corrélation:

 #include  #include  #include  #include  #include  // 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("\np: %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("\ni: %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)); } 

En utilisant un BLAS rapide (via Revolution R ou Goto BLAS), vous pouvez calculer rapidement toutes ces corrélations en R sans écrire de code C. Sur mon PC Intel i7 de première génération, cela prend 16 secondes:

 n = 400; m = 1e6; # Generate data mat = masortingx(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) 

Le code R ci-dessus indique le minutage suivant:

  user system elapsed 31.82 1.98 15.74 

J’utilise cette approche dans mon paquet MasortingxEQTL .
http://www.bios.unc.edu/research/genomic_software/Masortingx_eQTL/

Plus d’informations sur les différentes options BLAS pour R sont disponibles ici:
http://www.bios.unc.edu/research/genomic_software/Masortingx_eQTL/runit.html#large

Quelques choses.

Tout d’abord, si vous utilisez l’interface .C pour les appels externes, il copie par défaut tous les arguments. C’est pourquoi l’object corr n’est pas en cours de modification. Si vous voulez éviter cela, vous devez définir DUP = false dans l’appel .C. Cependant, en général, utiliser .C pour modifier des objects R existants n’est pas la méthode préférée pour faire les choses. Au lieu de cela, vous souhaiterez probablement créer un nouveau tableau et permettre à l’appel externe de le remplir, comme ceci.

 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)) 

Deuxièmement, en ce qui concerne l’écriture d’une fonction de corrélation rapide, la première chose à essayer est de comstackr R avec une implémentation BLAS efficace. Cela ne fera pas que rendre votre fonction de corrélation plus rapide, cela rendra toute votre algèbre linéaire plus rapide. Les bons candidats libres sont ACML d'AMD ou ATLAS. Chacun de ceux-ci sera capable de calculer des masortingces de corrélation très rapidement. L'accélération est plus qu'une simple parallélisation - ces bibliothèques sont également intelligentes en matière d'utilisation du cache et sont optimisées au niveau de l'assemblage. Ainsi, même avec un seul cœur, vous constaterez une grande amélioration. http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acml/ http://math-atlas.sourceforge.net/

Enfin, si vous voulez vraiment écrire votre propre code C, je suggérerais d'utiliser openMP pour scinder automatiquement le calcul entre différents threads, plutôt que de le faire à la main. Mais pour quelque chose d'aussi fondamental que la multiplication masortingcielle, il est probablement préférable d'utiliser une bibliothèque optimisée disponible.