I'm implementing a package in R that uses Eigen (v3.3.5) to do various calculations related to Markov chains with very large sparse matrices. Across the board, Eigen has provided a reasonable performance boost, but I'm seeing an issue with several calculations that involve matrix powers. An example of one these calculations is as simple as M^t, where M is a sparse matrix and t is a number of time steps. Since calculating this directly would produce an impractically large dense matrix, I'm using elementary vectors to get specific rows and columns via vM^t and (M^t)v. This then let's me do something like this: vM^t = vMMMMMM.... = ((((vM)M)M)M)..., which gets around having to calculate a dense matrix. Overall, it works as expected, but for some reason it suddenly slows down somewhere between 45000 and 50000 time steps, like so:

https://imgur.com/a/LDcg41U

Regardless of the matrix dimensions, this inflection is occurring at the same number of time steps (overall time does still increase with larger dimensions), and I'm not sure why. It's happening across all my functions that use a similar approach for dealing with matrix powers. Here's the code:

Code: Select all

Rcpp::NumericVector qpow_row(Eigen::Map<Eigen::SparseMatrix<double> > &M, const int row, const int steps)
{
  Eigen::RowVectorXd res = M.row(row-1); // Subtract 1 because R indexes start at 1

  for(int i = 1; i < steps; i++) {
    res = res * M;
  }

  return Rcpp::wrap(res);
}


So, trying various things to rule out things like unintentionally modifying my matrix, cache issues, etc, I narrowed the issue down to my result vector. Despite the values from my calculations seeming perfectly reasonable and normal in terms of hardware limits, I decided to try this handy-dandy bit of code to monitor the result at varying time steps:

Code: Select all

  int inf_count = 0;
  int nan_count = 0;
  int norm_count = 0;
  int subnorm_count = 0;
  int zero_count = 0;
  int unk_count = 0;

  for (int i = 0; i < res.cols(); i++) {
    switch(std::fpclassify(res(i))) {
    case FP_INFINITE:  inf_count++; break;
    case FP_NAN:       nan_count++; break;
    case FP_NORMAL:    norm_count++; break;
    case FP_SUBNORMAL: subnorm_count++; break;
    case FP_ZERO:      zero_count++; break;
    default:           unk_count++; break;
    }
  }

  Rcpp::Rcout << inf_count << " inf numbers\n";
  Rcpp::Rcout << nan_count << " nan numbers\n";
  Rcpp::Rcout << norm_count << " normal numbers\n";
  Rcpp::Rcout << subnorm_count << " subnormal numbers\n";
  Rcpp::Rcout << zero_count << " zero numbers\n";
  Rcpp::Rcout << unk_count << " unknown numbers\n\n";


credit: https://en.cppreference.com/w/cpp/numeric/math/fpclassify

Sure enough, right around the 47000th step, the values in the vector go from all normal to denormal/subnormal, resulting in a 50x slowdown. What I have to figure out is why these numbers are considered subnormal, given that I've see numbers as large as 1e-11 be considered subnormal, which is well within what the exponent for doubles supports. Ultimately, it's just an issue of me needing to gain a better understanding of subnormal numbers and whether simd hardware has any special issues I need to be aware of.

Are there any recommendations for how to handle this issue with Eigen? There appears to be a few general options listed on Intel's site here: https://software.intel.com/en-us/cpp-compiler-developer-guide-and-reference-denormal-numbers