Microbenchmarks: Matrix R package

Conclusion: Any operation or subview in C++ using Rcpp::dgCMatrix is always faster than the equivalent operation using R in the Matrix package (except for retrieving large non-contiguous subviews).

For convenience, the microbenchmarks below use a single 1000 x 1000 random square matrix that is 10% dense. Results are comparable for other densities and sizes.

Element access

Retrieve a single element at random, a random subset of 100 non-contiguous elements within a single row or column, a random non-contiguous subview of 100 rows and columns, or a single row or a single column.

C++

#include <RcppSparse.h>

//[[Rcpp::export]]
double Rcpp_SingleElement(Rcpp::dgCMatrix& mat, int row, int col){
  return mat(row, col);
}
  
//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_RowSubview(Rcpp::dgCMatrix& mat, Rcpp::IntegerVector& rows, int col){
  return mat(rows, col);
}

//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_ColSubview(Rcpp::dgCMatrix& mat, int row, Rcpp::IntegerVector& cols){
  return mat(row, cols);
}

//[[Rcpp::export]]
Rcpp::NumericMatrix Rcpp_RowColSubview(Rcpp::dgCMatrix& mat, Rcpp::IntegerVector& rows, Rcpp::IntegerVector& cols){
  return mat(rows, cols);
}

//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_Row(Rcpp::dgCMatrix& mat, int i){
  return mat.row(i);
}


//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_Col(Rcpp::dgCMatrix& mat, int i){
  return mat.col(i);
}

R

library(microbenchmark)
library(Matrix)
mat <- rsparsematrix(1000, 1000, 0.1)

microbenchmark(
  "Rcpp_SingleElement" = Rcpp_SingleElement(mat, sample(1:1000, 1), sample(1:1000, 1)),
  "Matrix_SingleElement" = mat[sample(1:1000, 1), sample(1:1000, 1)],
  "Rcpp_RowSubview" = Rcpp_RowSubview(mat, sample(1:1000, 100), sample(1:1000, 1)),
  "Matrix_RowSubview" = mat[sample(1:1000, 100), sample(1:1000, 1)],
  "Rcpp_ColSubview" = Rcpp_ColSubview(mat, sample(1:1000, 1), sample(1:1000, 100)),
  "Matrix_ColSubview" = mat[sample(1:1000, 1), sample(1:1000, 100)],
  "Rcpp_RowColSubview" = Rcpp_RowColSubview(mat, sample(1:1000, 100), sample(1:1000, 100)),
  "Matrix_RowColSubview" = mat[sample(1:1000, 100), sample(1:1000, 100)],
  "Rcpp_Row" = Rcpp_Row(mat, sample(1:1000, 1)),
  "Matrix_Row" = mat[sample(1:1000, 1), ],
  "Rcpp_Col" = Rcpp_Col(mat, sample(1:1000, 1)),
  "Matrix_Col" = mat[, sample(1:1000, 1)],
  times = 100
)

Results

Unit: microseconds
                 expr   min     lq    mean median     uq     max neval
   Rcpp_SingleElement   9.3  10.65  12.878  12.70  14.75    19.7   100
 Matrix_SingleElement 222.2 226.75 231.167 229.40 233.70   277.2   100
      Rcpp_RowSubview  20.2  22.45  24.452  25.15  26.20    29.2   100
    Matrix_RowSubview 236.3 241.05 435.863 243.50 246.70 19344.6   100
      Rcpp_ColSubview  19.1  22.55  24.569  24.40  26.65    37.1   100
    Matrix_ColSubview 253.6 257.95 265.068 261.05 264.90   353.7   100
   Rcpp_RowColSubview 445.2 458.50 469.096 466.45 474.25   525.0   100
 Matrix_RowColSubview 238.5 242.80 247.123 244.95 247.60   285.4   100
             Rcpp_Row  15.2  36.90  46.700  45.25  57.10    75.9   100
           Matrix_Row 392.7 399.50 423.887 405.05 420.30  1522.5   100
             Rcpp_Col   6.2   8.00   9.792   9.25  11.30    22.8   100
           Matrix_Col 217.9 221.65 227.109 225.50 229.60   342.2   100

Marginal totals

C++

#include <RcppSparse.h>

//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_colSums(Rcpp::dgCMatrix& mat){
  return mat.colSums();
}
  
//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_colMeans(Rcpp::dgCMatrix& mat){
  return mat.colMeans();
}

//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_rowSums(Rcpp::dgCMatrix& mat){
  return mat.rowSums();
}

//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_rowMeans(Rcpp::dgCMatrix& mat){
  return mat.rowMeans();
}

R

library(microbenchmark)
library(Matrix)
mat <- rsparsematrix(1000, 1000, 0.1)

microbenchmark(
  "Rcpp::colSums"    = Rcpp_colSums(mat),
  "Matrix::colSums"  = Matrix::colSums(mat),
  "Rcpp::rowSums"    = Rcpp_rowSums(mat),
  "Matrix::rowSums"  = Matrix::rowSums(mat),
  "Rcpp::colMeans"   = Rcpp_colMeans(mat),
  "Matrix::colMeans" = Matrix::colMeans(mat),
  "Rcpp::rowMeans"   = Rcpp_rowMeans(mat),
  "Matrix::rowMeans" = Matrix::rowMeans(mat),
  times = 100
)

Results

Unit: microseconds
             expr    min      lq     mean  median      uq    max neval
    Rcpp::colSums  275.5  276.60  280.459  277.80  281.25  316.1   100
  Matrix::colSums  360.5  365.55  371.755  370.15  374.50  402.6   100
    Rcpp::rowSums  281.2  282.95  288.330  285.25  291.15  343.0   100
  Matrix::rowSums 1199.9 1228.40 1257.487 1237.60 1271.10 1442.2   100
   Rcpp::colMeans  286.2  293.65  305.597  295.50  298.15 1168.6   100
 Matrix::colMeans  360.4  366.45  372.733  370.50  374.30  430.0   100
   Rcpp::rowMeans  307.8  309.85  313.908  311.50  314.50  361.1   100
 Matrix::rowMeans 1216.6 1234.10 1272.008 1250.35 1293.95 1545.9   100