Rcpp Armadillo speedup with large sparse matrix multiplications

A student of mine has noticed some strange behavior with some code we have implemented in Rcpp Armadillo. Specifically, when we use matrix multiplication with large sparse matrices, there appears to be a dimensionality at which the implementation actually runs more quickly, despite larger dimensions of the matrices being multiplied.

What is the explanation for this? Can we modify our implementation in order to see the improved speed at lower dimensions?

We wrote a simple example to illustrate this issue.

library(MASS)
library(Matrix)
library(Rcpp)
sourceCpp("MatMultFunc.cpp")
nreps <- 50
p <- 500
nSeq <- seq(80, 300, by = 20)
storeTime <- matrix(0, nreps, length(nSeq))
set.seed(1)

for (replicate in 1:nreps) {
  for (n in nSeq) {
  
    X <- matrix(rnorm(n*p), nrow = n)
    l <- n*(n-1)/2
    D <- matrix(0, nrow = l, ncol = n)

    counter <- 1
    for (j in 1:(n-1)) {
      for (k in (j+1):n) {
          D[counter, j] <- 1
          D[counter, k] <- -1
          counter <- counter + 1
      }
    }
    
    
    D <- Matrix(D, sparse=TRUE)
    
    Beta <- rep(0, p)
    Beta[sample(1:p, 50)] <- 1
    Beta <- Matrix(Beta, sparse = TRUE)
    
    ptm <- proc.time()
    RcppMatProds <- MatMultFunc(D, X, Beta)
    t0 <- proc.time() - ptm

    storeTime[replicate, which(nSeq == n)] <- t0[3]

    cat("Replicate: ", replicate, "; n = ", n,"\n")

  }
}

plot(nSeq, colMeans(storeTime), type = "l", ylab="Average runtime in seconds", xlab="n")

The C++ function we are timing is the following:

#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::depends("RcppArmadillo")]]

// [[Rcpp::export]]
List MatMultFunc(arma::sp_mat D, arma::mat X, arma::sp_mat BetaSp){

    arma::vec tXB(D * (X * BetaSp));

    for (int m = 0; m < 1000; m++) {

        tXB = D * (X * BetaSp);
    }

     return List::create(Named("tXB") = wrap(tXB));
}