sparse matrix - умножение разреженной матрицы.улучшение производительности в R - PullRequest
Новая
       31

sparse matrix - умножение разреженной матрицы.улучшение производительности в R

0 голосов
/ 08 октября 2018

R моего mac связан с openblas.Когда я смотрю на использование «% CPU» при выполнении разреженного умножения в R или в Rcpp с использованием Armadillo , не похоже, что многопоточность используется в отличие от плотно-плотного умножения.С точки зрения скорости, умножение одного потока на R или Armadillo в одном потоке также медленнее, чем в Matlab.

Чтобы решить эту проблему, я реализовал алгоритм Ф. Г. Густавсона (https://dl.acm.org/citation.cfm?id=355796) для выполнения умножения разреженных матриц в Rcpp с использованием контейнера spMat от Armadillo.

Я вижу улучшение (см. ниже), если я игнорирую упорядочение строк, которое является прямой реализацией алгоритма, однако стандартное упорядочение делает его более медленным, чем R ( отредактировано согласно комментарию mtall ). Я не эксперт в Rcpp /RcppArmadillo / C ++ и я ищу помощи в двух конкретных вещах:

  • Программно, как я могу сделать функцию sp_sp_gc_ord более эффективной и быстрой на основе однопоточного приложения?

  • Моя неудачная попытка многопоточности sp_sp_gc_ord с openmp вызывает сбой R. Я прокомментировал omp Команды ниже.Я посмотрел обсуждения галереи Rcpp по OpenMP http://gallery.rcpp.org/tags/openmp/, но не смог выяснить проблему

Я woulБуду признателен за любую помощь.Ниже приведен воспроизводимый пример кода и соответствующего микробенчмарка: 

#### Rcpp functions

#include <RcppArmadillo.h>
#include<omp.h>
#include<Rcpp.h>

using namespace Rcpp;
using namespace arma;

// [[Rcpp::plugins(openmp)]]
// [[Rcpp::depends(RcppArmadillo)]]


// [[Rcpp::export]]
sp_mat sp_sp_gc_ord(const arma::sp_mat &A, const arma::sp_mat &B, double p){

  // This function evaluates A * B where both A & B are sparse and the resultant
  // product is also sparse

  // define matrix sizes
  const int mA= A.n_rows;
  const int nB= B.n_cols;

  // number of non-zeros in the resultant matrix
  const int nnzC = ceil(mA * nB * p);

  // initialize colptr, row_index and value vectors for the resultant sparse matrix
  urowvec colptrC(nB+1);
  colptrC.zeros();
  uvec rowvalC(nnzC);
  rowvalC.zeros();
  colvec nzvalC(nnzC);

  //setenv("OMP_STACKSIZE","500M",1);

  // counters and other variables
  unsigned int i, jp, j, kp, k, vp; 
  unsigned int ip = 0;
  double nzB, nzA; 
  ivec xb(mA);
  xb.fill(-1);
  vec x(mA);

  // loop logic: outer loop over columns of B and inner loop over columns of A and then aggregate

  //  #pragma omp parallel for shared(colptrC,rowvalC,nzvalC,x,xb,ip,A,B) private(j,nzA,nzB,kp,i,jp,kp,k,vp) default(none) schedule(auto) 
  for(i=0; i< nB; i++) { 

    colptrC.at(i) = ip;

    for ( jp = B.col_ptrs[i]; jp < B.col_ptrs[i+1]; jp++) {

      j = B.row_indices[jp];
      nzB = B.values[jp];

      for ( kp = A.col_ptrs[j]; kp < A.col_ptrs[j+1]; kp++ ){

        k = A.row_indices[kp];
        nzA = A.values[kp];

        if (xb.at(k) != i){
          rowvalC.at(ip) = k;
          ip +=1;
          // Rcpp::print(wrap(ip));
          xb.at(k) = i;
          x.at(k) = nzA * nzB;
        } else {
          x.at(k) += nzA * nzB;
        }
      }
    }

    // put in the value vector of resultant matrix

    if(ip>0){        
      for ( vp= colptrC.at(i); vp <= (ip-1); vp++ ) {
        nzvalC.at(vp) = x(rowvalC.at(vp));
      }

    }
  }

  // resize and put in the spMat container
  colptrC.at(nB) = ip;
  sp_mat C(rowvalC.subvec(0,(ip-1)),colptrC,nzvalC.subvec(0,(ip-1)),mA,nB);

  // Gustavson's algorithm produces unordered rows for each column: a standard way to address this is: (X.t()).t()

  return (C.t()).t();
}  


 // [[Rcpp::export]]
sp_mat sp_sp_arma(const sp_mat &A, const sp_mat &B){

  return A * B; 

} 


// [[Rcpp::export]]
mat dense_dense_arma(const mat &A, const mat &B){

  return A * B; 

} 

#### End

Соответствующая часть микробенчмарка в R: 

#### Microbenchmark 

library(Matrix)
library(microbenchmark) 

## define two matrices
m<- 1000
n<- 6000
p<- 2000

A<-  matrix(runif(m*n),m,n)
B<-  matrix(runif(n*p),n,p)
A[abs(A)> .01] = B[abs(B)> .01] = 0
A <- as(A,'dgCMatrix')
B<- as(B,'dgCMatrix')
Adense<- as.matrix(A)
Bdense<- as.matrix(B)

## sp_sp_gc is the function without ordering 

microbenchmark(sp_sp_gc(A,B,.5),sp_sp_arma(A,B),A%*%B,
dense_dense_arma(Adense,Bdense),Adense %*% Bdense,Adense %*% B, times=100)

Unit: milliseconds
                             expr       min        lq      mean    median        uq       max neval
              sp_sp_gc(A, B, 0.5)  16.09809  21.75001  25.76436  24.44657  26.96300  99.30778   100
          sp_sp_gc_ord(A, B, 0.5)  36.78781  44.64558  49.82102  47.64348  51.87361 116.85013   100
                 sp_sp_arma(A, B)  47.45203  52.77132  59.37077  59.24010  62.41710  86.15647   100
                          A %*% B  23.64307  28.99649  32.88566  32.10017  35.21816  59.16251   100
 dense_dense_arma(Adense, Bdense) 286.22358 302.95170 345.66766 317.75786 340.50143 862.15116   100
                Adense %*% Bdense 292.32099 317.10795 342.48345 329.80950 342.21333 697.56468   100
                     Adense %*% B 167.87248 186.63499 219.11872 195.19197 212.50286 843.17172   100   
####

sessionInfo (): 

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /usr/local/Cellar/openblas/0.3.3/lib/libopenblas_haswellp-r0.3.3.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Matrix_1.2-14           RcppArmadillo_0.8.500.0 Rcpp_0.12.18           

loaded via a namespace (and not attached):
[1] compiler_3.5.1  grid_3.5.1      lattice_0.20-35

Rcpp и RcppArmadillo устанавливаются из исходного кода после установки clang4 для Mac по ссылке Coatless https://github.com/coatless/r-macos-rtools

...