fast large matrix multiplication in R

Question

I have two matrices in R that I want to multiply:

a = matrix(rnorm(20*10000, mean=0, sd=5), 20, 10000)
b = matrix(rnorm(20*10000, mean=0, sd=5), 20, 10000)
t(a)%*%b

Given that the dimension in larger this matrix multiplication takes a lot of time, is there a specific way to make computations faster ? And are there any build in functions in R to make such multiplications faster?

Answer 1

There are many ways to approach this depending upon your code, effort, and hardware.

1. Use the 'best' function for the job

The simplest is to use crossprod which is the same as t(a)%*% b (Note - this will only be a small increase in speed)

crossprod(a,b) 

2. Use Rcpp (and likely RcppEigen/RcppArmadillo).

C++ will likely greater increase the speed of your code. Using linear algebra libraries will likely also help this further (hence Eigen and Armadillo). This however assumes you are willing to write some C++.

3. Use a better backend

After this point you are looking at BLAS backends such as OpenBLAS, Atlas, etc. Hooking these up to R varies depending upon your OS. It is quite easy if you are using a Debian system like Ubuntu. You can find a demo here. These can sometimes be leveraged further by libraries such as Armadillo and Eigen.

4. GPU Computing

If you have a GPU (e.g. AMD, NVIDIA, etc.) you can leverage the many cores within to greatly speed up your computations. There are a few that could be useful including gpuR, gputools, and gmatrix

EDIT - to address @jenesaiquoi comment on benefit of Rcpp

test.cpp

// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
    arma::mat C = A * B;

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, Eigen::MatrixXd B){
    Eigen::MatrixXd C = A * B;

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
    Eigen::MatrixXd C = A * B;

    return Rcpp::wrap(C);
}

test.R

library(Rcpp)

A <- matrix(rnorm(10000), 100, 100)
B <- matrix(rnorm(10000), 100, 100)

library(microbenchmark)
sourceCpp("test.cpp")
microbenchmark(A%*%B, armaMatMult(A, B), eigenMatMult(A, B), eigenMapMatMult(A, B))

Unit: microseconds
                  expr     min       lq     mean   median       uq      max neval
               A %*% B 885.846 892.1035 933.7457 901.1010 938.9255 1411.647   100
     armaMatMult(A, B) 846.688 857.6320 915.0717 866.2265 893.7790 1421.557   100
    eigenMatMult(A, B) 205.978 208.1295 233.1882 217.0310 229.4730  369.369   100
 eigenMapMatMult(A, B) 192.366 194.9835 207.1035 197.5405 205.2550  366.945   100

Answer 2

To add to cdeterman's answer: You can use eigen's build in parallelization for dense matrix products. In order to do so, you need to compile with open mp activated.

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

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

// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
  arma::mat C = A * B;
  
  return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, 
                  Eigen::MatrixXd B, 
                  int n_cores){
  
  Eigen::setNbThreads(n_cores);
  //qDebug()  << Eigen::nbThreads( );
  Eigen::MatrixXd C = A * B;
  
  return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMapMatMult2(const Eigen::Map<Eigen::MatrixXd> A,
                      Eigen::Map<Eigen::MatrixXd> B, 
                      int n_cores){
  
  Eigen::setNbThreads(n_cores);
  Eigen::MatrixXd C = A * B;
  return Rcpp::wrap(C);
}

Here are some benchmarks:

Note that if N = k = 100, parallelization does not necessarily improve performance. If the matrix dimensions get larger, parallelization starts to have an impact (N = k = 1000):

library(microbenchmark)

# Benchmark 1: N = k = 100

N <- 100
k <- 100

A <- matrix(rnorm(N*k), N, k)
B <- matrix(rnorm(N*k), k, N)

microbenchmark(A%*%B, 
               armaMatMult2(A, B),
               eigenMatMult2(A, B, n_cores = 1),
               eigenMatMult2(A, B, n_cores = 2),
               eigenMatMult2(A, B, n_cores = 4),
               eigenMapMatMult2(A, B, n_cores = 1),
               eigenMapMatMult2(A, B, n_cores = 2),
               eigenMapMatMult2(A, B, n_cores = 4), 
               times = 100

# Unit: microseconds
#                                 expr   min     lq    mean median     uq   max neval
#                              A %*% B 535.6 540.75 552.594 551.25 554.50 650.2   100
#                   armaMatMult2(A, B) 542.0 549.10 560.975 556.35 560.25 738.1   100
#     eigenMatMult2(A, B, n_cores = 1) 147.1 152.65 159.165 159.65 162.90 180.5   100
#     eigenMatMult2(A, B, n_cores = 2)  97.1 109.90 124.496 119.60 127.50 391.8   100
#     eigenMatMult2(A, B, n_cores = 4)  71.7  88.15 155.220 115.55 216.95 507.3   100
#  eigenMapMatMult2(A, B, n_cores = 1) 139.1 150.10 154.889 154.20 158.35 244.3   100
#  eigenMapMatMult2(A, B, n_cores = 2)  93.4 105.70 116.808 113.55 120.40 323.7   100
#  eigenMapMatMult2(A, B, n_cores = 4)  66.8  82.60 161.516 196.25 210.40 598.9   100
)

# Benchmark 2: N = k = 1000

N <- 1000
k <- 1000

A <- matrix(rnorm(N*k), N, k)
B <- matrix(rnorm(N*k), k, N)

microbenchmark(A%*%B, 
                armaMatMult2(A, B),
                eigenMatMult2(A, B, n_cores = 1),
                eigenMatMult2(A, B, n_cores = 2),
                eigenMatMult2(A, B, n_cores = 4),
               eigenMapMatMult2(A, B, n_cores = 1),
               eigenMapMatMult2(A, B, n_cores = 2),
               eigenMapMatMult2(A, B, n_cores = 4), 
               times = 100
)


Unit: milliseconds
                                expr      min        lq      mean    median        uq
                             A %*% B 597.1293 605.56840 814.52389 665.86650 1025.5896
                  armaMatMult2(A, B) 603.3894 620.25675 830.98947 693.22355 1078.4853
    eigenMatMult2(A, B, n_cores = 1) 131.4696 135.22475 186.69826 193.37870  219.8727
    eigenMatMult2(A, B, n_cores = 2)  67.8948  71.71355 114.52759  74.17380  173.3060
    eigenMatMult2(A, B, n_cores = 4)  41.8564  48.87075  79.55535  72.00705  106.8572
 eigenMapMatMult2(A, B, n_cores = 1) 125.3890 129.26125 175.09933 177.23655  213.0536
 eigenMapMatMult2(A, B, n_cores = 2)  62.2866  65.78785 115.74248  79.92470  167.0217
 eigenMapMatMult2(A, B, n_cores = 4)  35.2977  40.42480  68.21669  63.13655   97.2571
       max neval
 1217.6475   100
 1446.5127   100
  419.2043   100
  217.9513   100
  139.9629   100
  298.2859   100
  230.6307   100
  118.2553   100