Инверсия матрицы и транспонирование R против c ++

Question

У меня есть расчетный поток мощности по сети.Тем не менее, я обнаружил, что для больших сетей расчет идет медленно.Я попытался реализовать алгоритм с помощью RcppArmadillo.Функция Rcpp в несколько раз быстрее для небольших сетей / матриц, но становится такой же медленной для более крупных.Скорость операции влияет на полезность функций в дальнейшем.Я плохо пишу функции или это только те опции, которые у меня есть?

Приведенный ниже код R и пример времени выполнения для матрицы 10, 100, 1000.

library(igraph); library(RcppArmadillo)

#Create the basic R implementation of the equation
flowCalc <- function(A,C){

  B <- t(A) %*% C %*% A
  Imp <- solve(B)
  PTDF <- C %*% A %*% Imp
  Out <- list(Imp = Imp, PTDF= PTDF)
  return(Out)

}

#Create the c++ implementation
txt <- 'arma::mat Am = Rcpp::as< arma::mat >(A);
   arma::mat Cm = Rcpp::as< arma::mat >(C);
   arma::mat B = inv(trans(Am) * Cm * Am);
   arma::mat PTDF = Cm * Am * B;
   return Rcpp::List::create( Rcpp::Named("Imp") = B ,
                           Rcpp::Named("PTDF") = PTDF ) ; '

flowCalcCpp <- cxxfunction(signature(A="numeric",
    C="numeric"),
    body=txt,
    plugin="RcppArmadillo")


#Create a function to generate dummy data
MakeData <- function(n, edgep){#make an erdos renyi graph of size x
    g <- erdos.renyi.game(n, edgep)
    #convert to adjacency matrix
    Adjmat <- as_adjacency_matrix(g, sparse = F)
    #create random graph and mask the elements with not edge
    Cmat <- matrix(rnorm(n*n), ncol = n)*Adjmat
    ##Output a list of the two matrices
    list(A = Adjmat, C = Cmat)
}

#generate dummy data
set.seed(133)
Data10 <- MakeData(10, 1/1)
Data100 <- MakeData(100, 1/10)
Data1000 <- MakeData(1000, 1/100)

#Compare results
BenchData <- microbenchmark(
               R10 = flowCalc(Data10$A, Data10$C),
               R100 = flowCalc(Data100$A, Data100$C),
               R1000 = flowCalc(Data1000$A, Data1000$C),
               Cpp10 = flowCalcCpp(Data10$A, Data10$C),
               Cpp100 = flowCalcCpp(Data100$A, Data100$C),
               Cpp1000 = flowCalcCpp(Data1000$A, Data1000$C))

Результаты по скорости показаны ниже.

РЕДАКТИРОВАТЬ:

После прочтения ответа Ральфа и комментариев Дирка я использовал https://cran.r -project.org / web / packages / gcbd / vignettes/gcbd.pdf, чтобы получить лучший обзор различий в реализациях BLAS

Затем я использовал руководство Дирка по установке реализации Microsoft BLAS https://github.com/eddelbuettel/mkl4deb (ясно, что сейчас я живу в Eddelverse)

Как только я это сделал, я установил ArrayFire и RcppArrayFire, как предложил Ральф.Затем я запустил код и получил следующие результаты:

Unit: microseconds
    expr        min          lq        mean      median          uq         max neval
     R10     37.348     83.2770    389.6690    143.9530    181.8625   25022.315   100
    R100    464.148    587.9130   1187.3686    680.8650    849.0220   32602.678   100
   R1000 143065.901 160290.2290 185946.5887 191150.2335 201894.5840  464179.919   100
   Cpp10     11.946     30.6120    194.8566     55.6825     74.0535   13732.984   100
  Cpp100    357.880    452.7815    987.8059    496.9520    554.5410   39359.877   100
 Cpp1000 102949.722 124813.9535 136898.4688 132852.9335 142645.6450  214611.656   100
    AF10    713.543    833.6270   1135.0745    966.5920   1057.4175    8177.957   100
   AF100   2809.498   3152.5785   3662.5323   3313.0315   3569.7785   12581.535   100
  AF1000  77905.179  81429.2990 127087.2049  82579.6365  87249.0825 3834280.133   100

Скорость уменьшается для меньших матриц, в два раза быстрее для матриц порядка 100, но почти в 10 раз быстрее для больших матриц,что согласуется с результатами Ральфа.Разница с использованием C ++ также выросла.

По моим данным, при обновлении BLAS стоит использовать версию C ++, но, возможно, не версию Arrayfire, поскольку мои матрицы недостаточно велики.

Answer 1

Я попробовал ваш код на двухъядерном ноутбуке под управлением Debian Linux и R 3.5.1, связанных с OpenBLAS.Кроме того, я также использовал RcppArrayFire для выполнения этих вычислений с использованием графического процессора (ничего особенного, просто встроенная графика).Вот результаты теста:

Unit: microseconds
    expr        min         lq        mean     median         uq        max neval  cld
     R10     55.488     90.142    117.9538    133.617    138.009    161.604    10 a   
    R100    698.883    711.488   1409.8286    773.352    901.339   5647.804    10 a   
   R1000 198612.971 213531.557 225713.6705 214993.916 226526.513 306675.313    10    d
   Cpp10     16.157     31.478     38.3473     40.529     51.122     52.351    10 a   
  Cpp100    519.527    544.728    573.0099    570.956    610.985    613.400    10 a   
 Cpp1000 174121.655 184865.003 196224.3825 193142.715 207066.037 223900.830    10   c 
    AF10    700.805    790.203   1469.4980   1347.639   1410.717   3824.905    10 a   
   AF100   1376.737   1562.675   1818.8094   1898.797   1949.364   2222.889    10 a   
  AF1000 106398.431 109130.482 118817.6704 118612.679 123193.649 139474.579    10  b  

При наименьшем размере (R10 и Cpp10) ваша машина быстрее моей.Но уже на R100 и Cpp100, в частности на R1000 и Cpp1000, я получаю более быстрое выполнение.Как отметил Дирк в комментариях, вы должны изучить оптимизированные / параллельные реализации BLAS.В Debian Linux (и производных, таких как Ubuntu) это так же просто, как

  sudo apt-get install libopenblas-base

См. Также здесь .Теперь о результатах с использованием графического процессора: это довольно типично.Для более мелких матриц это хуже, чем как у основания R, так и у броненосца.Использование графического процессора приводит к значительным накладным расходам при перемещении данных между основной памятью и памятью графического процессора.Но для самого большого размера параллельное выполнение на графическом процессоре перевешивает эти издержки, и выигрыш в производительности довольно приятный.

Вот мой код для справки.Я позволил себе обновить inline::cxxfunction до Rcpp::cppFunction:

#Create the basic R implementation of the equation
flowCalc <- function(A,C){
  B <- t(A) %*% C %*% A
  Imp <- solve(B)
  PTDF <- C %*% A %*% Imp
  Out <- list(Imp = Imp, PTDF= PTDF)
  return(Out)

}

# Create Armadillo function
Rcpp::cppFunction(depends = "RcppArmadillo", code = '
Rcpp::List flowCalcCpp(const arma::mat &Am, const arma::mat &Cm) {
   arma::mat B = inv(trans(Am) * Cm * Am);
   arma::mat PTDF = Cm * Am * B;
   return Rcpp::List::create( Rcpp::Named("Imp") = B ,
                           Rcpp::Named("PTDF") = PTDF );
}')

# Create ArrayFire function
Rcpp::cppFunction(depends = "RcppArrayFire", code = '
Rcpp::List flowCalcAF(const RcppArrayFire::typed_array<f32> &A, 
                      const RcppArrayFire::typed_array<f32> &C) {
  af::array B = af::inverse(af::matmul(af::matmulTN(A, C), A));
  af::array PTDF = af::matmul(af::matmul(C, A), B);
  return Rcpp::List::create( Rcpp::Named("Imp") = B ,
                             Rcpp::Named("PTDF") = PTDF );
}')



library(igraph)
#Create a function to generate dummy data
MakeData <- function(n, edgep){#make an erdos renyi graph of size x
  g <- erdos.renyi.game(n, edgep)
  #convert to adjacency matrix
  Adjmat <- as_adjacency_matrix(g, sparse = F)
  #create random graph and mask the elements with not edge
  Cmat <- matrix(rnorm(n*n), ncol = n)*Adjmat
  ##Output a list of the two matrices
  list(A = Adjmat, C = Cmat)
}

#generate dummy data
set.seed(133)
Data10 <- MakeData(10, 1/1)
Data100 <- MakeData(100, 1/10)
Data1000 <- MakeData(1000, 1/100)

#Compare results
microbenchmark::microbenchmark(
  R10 = flowCalc(Data10$A, Data10$C),
  R100 = flowCalc(Data100$A, Data100$C),
  R1000 = flowCalc(Data1000$A, Data1000$C),
  Cpp10 = flowCalcCpp(Data10$A, Data10$C),
  Cpp100 = flowCalcCpp(Data100$A, Data100$C),
  Cpp1000 = flowCalcCpp(Data1000$A, Data1000$C),
  AF10 = flowCalcAF(Data10$A, Data10$C),
  AF100 = flowCalcAF(Data100$A, Data100$C),
  AF1000 = flowCalcAF(Data1000$A, Data1000$C),
  times = 10)