8

I'm writing an R package using Rcpp and RcppEigen and I'm having a problem with matrix slicing and sub-setting. I need to get an off-diagonal square matrix from a larger square matrix

The Eigen::MatrixXd slicing methods seem to be the problem. Using Eigen::seq() method doesn't work at all from R because "no member named 'seq' in namespace 'Eigen'" and the MatrixXd.block(i, j, n, n) method is causing a crash with certain large(ish) values of n. Sometimes it works perfectly fine, but if I increase the size, it causes a fatal crash.

Here's an example of the C++ code:

// [[Rcpp::depends(RcppEigen)]]
#include <iostream>
#include <RcppEigen.h>

using namespace Rcpp;

using Eigen::Map;
using Eigen::MatrixXd;

typedef Map<MatrixXd> MapMatd;

// [[Rcpp::export]]
List crosspart_worker_cpp(const MapMatd& Vij, ...){

  int n = Vij.cols()/2;

  /* ... some arbitrary code ... */
 

  // extract sub-block of varcovar matrix (only unique pairs)
  MatrixXd Vsub = Vij.block(1, n + 1, n, n);

  /* ... more code ... */

  List out_lst = List::create(Named("Vsub") = Vsub, ...);

  return out_lst;
}

and R code that uses it:

# test_crosspart-worker-cpp.R

# setup ----
## source relevant file
  # normally build and load package instead of sourceCpp()
if(!exists("crosspart_worker_cpp")){
  Rcpp::sourceCpp("crosspart-worker.cpp")
}

## make reproducible
set.seed(75)

## set parameters
n = 100 # rows in original model matrix X

## ... additional setup code ...

# Generate dummy data with arbitrary contents, but correct structure ----
## varcovar matrix
Vij <- matrix(abs(rnorm(2*n * 2*n)), nrow = 2*n)

## ... other code ...

# use the function ----
result <- crosspart_worker_cpp(Vij = Vij, ...)

Why doesn't this work consistently? Are there other options for sub-setting a matrix in RcppEigen?

My original post contained a larger C++ file and I didn't know where the error was. I've been able to identify the problem line of code and the function works perfectly if I remove it. Currently, I am obtaining the matrix subset in R and then passing that as an argument to the C++ function. However, it would be be highly beneficial for my package to do this in all in a C++ function.

I used Rcout to find the problem line of code:

MatrixXd Vsub = VDiag.block(1, n + 1, n, n);

Which is supposed to take the block matrix starting at row 1 and column np and contain np rows and columns according to Eigen: Slicing and Indexing

this is meant to replicate the R code Vsub[1:np, (n + 1):(2*n)].

Why would this cause a crash? Does the Eigen slice method not work in R? is there an RcppEigen specific way? I haven't found one online.

System specs:

OS: Microsoft Windows 10 Education Edition
Processor: AMD Ryzen 7 3700X 8-core
Installed Physical Memory (RAM): 16.0 GB

R version: 4.0.2 (2020-06-22) -- "Taking Off Again"
Rstudio version: 1.3.1056
morrowcj
  • 169
  • 6
  • 6
    Your first step should be to try and find out which line of code causes the issue. I'd add `Rcout` commands to do that. – Roland Aug 25 '20 at 06:20
  • 3
    In addition to the great advice from @Roland, it may be helpful to comment out most of your code, run it with the large data, if it doesn't crash, uncomment out a bit, run it, etc. At this point you can quickly find the line of code that is bombing. Now, you can add `Rcpp::print(Rcpp::wrap(objectOfInterest))` to get an idea of what your data structures look like. I have found this method to be very effective. – Joseph Wood Aug 25 '20 at 15:42
  • Re Eigen::seq ; maybe due to versions; is Eigen v3.4 required https://eigen.tuxfamily.org/dox-devel/group__TutorialSlicingIndexing.html and RcppEigen on cran is 3.3.4 – – user2957945 Sep 02 '20 at 09:23
  • That is a really good observation! I should have caught that. In that case, what are my options for slicing? It appears that the github version of `RcppEigen` is also not to 3.4 yet. Do I need to use RcppArmadillo for that task? Can I update the/use `Eigen` package that Rtools C++ compiler uses? – morrowcj Sep 02 '20 at 17:05

0 Answers0