Large Matrices in Rcpparmadillo via The Arma_64Bit_Word Define

Large Matrices in RcppArmadillo via the ARMA_64BIT_WORD define

1. RcppArmadillo is an Rcpp-only package due to the contents in the /src directory. To enable C++11 use // [[Rcpp::plugins(cpp11)]]
2. ARMA_64BIT_WORD is not defined. To define it add #define ARMA_64BIT_WORD 1 prior to the #include <RcppArmadillo.h>.

Sample implementation using sourceCpp()

#define ARMA_64BIT_WORD 1
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(cpp11)]] 

// [[Rcpp::export]] 
arma::mat cosCpp(const arma::mat& X) {

    arma::mat Y = arma::trans(X) * X; // matrix product
    arma::mat res = Y / (arma::sqrt(arma::diagvec(Y)) * arma::trans(arma::sqrt(arma::diagvec(Y))));

    return res;
}

To define it within the /src/Makevars{.win} for a package use:

PKG_CPPFLAGS = -DARMA_64BIT_WORD=1

Error message: Cube::init(): requested size is too large; suggest to enable ARMA_64BIT_WORD using sommer package for GWAS

My first take Ana is that you're trying to fit an unstructured multivariate model with 26 traits when you use cbind(), that means that if you have 1000 records, this will be a model of 309 x 26 = 8,034 records which would be a bit too big for the direct inversion algorithm that sommer uses, plus the number of parameters to estimate are a lot (think all the covariance parameters (26*25)/2 = 325. I would suggest fitting a GWAS per trait in a for loop to solve your issue. Unless you have a good justification to run a multivariate GWAS this is the issue with your analysis more than the C++ code behind. For example:

var_cov <- A.mat(m_matrix) ## aditive relationship matrix
traits <- c(DW20, PLA07, PLA08, PLA09, PLA10, PLA11, PLA12, PLA13, PLA14, PLA15, PLA16, PLA17, PLA18, RGR07_09, RGR08_10, RGR09_11, RGR10_12, RGR11_13, RGR12_14, RGR13_15, RGR14_16, RGR15_17, RGR16_18, SA, SL, SW)

for(itrait in traits){

model <- GWAS(as.formula(paste(itrait,"~1")), random = ~ vs(accession, Gu = var_cov), data = pheno2, M = m_matrix, gTerm = "u:accession", n.PC = 5)

}

If it turns out that even with a single trait the arma::cube function presents memory issues then definitely we need to look at why the armadillo library cannot deal with those dimensions.

Cheers,
Eduardo

Handling extremely large and sparse matrices in RcppArmadillo

The originally posted code was actually correct, but somehow the machine it was running on was not---see the discussion above for details.

As an add-on, here is a slightly edited version of the code without global namespace, actual access to the matrix (passing was enough, this is simply more explicit) and returning void. We also added the usual 'call it from R for me' trick one can do with Rcpp Attributes.

Code

// lightly edited version of question, working fine for me
#define ARMA_64BIT_WORD 1
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
void spmatTest(arma::sp_mat M) {
  // show a simple access of matrix,
  // here we just sum elements on diagonal
  Rcpp::Rcout << "Sum of diag is "
              << arma::as_scalar(arma::sum(M.diag()))
              << std::endl;
}

/*** R
spmatTest( Matrix::Diagonal( 1e6 ))
spmatTest( Matrix::Diagonal( 1e7 ))
spmatTest( Matrix::Diagonal( 1e8 ))
*/

Output

R> Rcpp::sourceCpp("~/git/stackoverflow/64428776/answer.cpp")

R> spmatTest( Matrix::Diagonal( 1e6 ))
Sum of diag is 1e+06

R> spmatTest( Matrix::Diagonal( 1e7 ))
Sum of diag is 1e+07

R> spmatTest( Matrix::Diagonal( 1e8 ))
Sum of diag is 1e+08
R> 

Large SpMat object with RcppArmadillo

Your matrix seems odd. By saying

 Matrix(rnorm(p),sparse=TRUE)

you end up with a p x 1 matrix, albeit sparse. If I just assign 10 rows or colums things
just work.

R> p <- 100000
R> X <- Matrix(rnorm(p),nrow=10,sparse=TRUE)
R> dim(X)
[1]    10 10000
R> norm2(X)
SpMat Xsp:
100832
NULL
R> 

So I think you just need a better sparse matrix to work with -- the conversion code, and Armadillo's sparse Matrix type, are fine.

Update on 2013-04-30: This was actually an Armadillo bug, which was just fixed upstream. A new RcppArmadillo verion 0.3.810.2 is now in SVN, and should migrate soon to CRAN shortly. You no longer need to define ARMA_64BIT_WORD.

passing Matrix::sparseMatrix to Rcpp

Ok, so I think I foud the answer. Basically armadillo throws this error if the total number of elements is larger than the size of the variable that stores the number of elements as seen here:
https://gitlab.com/conradsnicta/armadillo-code/-/blob/9.900.x/include/armadillo_bits/SpMat_meat.hpp

Someone realized this before and provided a solution here:
Large Matrices in RcppArmadillo via the ARMA_64BIT_WORD define

Why does Rcpp function segfault only when called from R package but not when sourced directly via sourceCpp?

I hinted in the comments that should try to decompose your problem into smaller ones. I just tried a little bit. Given the persistent

 error: arma::memory::acquire(): out of memory
 Error in Znew_gen2(U, Z, group, cols, n, q, d, Znew@address, J_SpMat) : 
   std::bad_alloc
 Execution halted

I added code to check your parameters, and indeed for

 Rcpp::Rcout << "n*nMC is " << n*nMC << ", J.ncols is " << J.n_cols << std::endl;

I see

 n*nMC is 50000, J.ncols is 25769803830

so currently your problem is not with a bigmemory object but rather the sparse matrix.

Edit a little while later: Turns out that the problem is likely your use
of #define ARMA_64BIT_WORD which makes the matrix dimensions be outside of int and hence the appearance of these YUGE values. If I removed that, your code runs.

So lesson re-learned: make the problem smaller and smaller, here it meant decompose your use of a sparse matrix from the use of bigmemory matrix.

Execution time of sparse matrix multiplication

No. The issue being documented here is two fold:

1. the conversion time required to transform a sparse matrix in R to an RcppArmadillo object is larger than Matrix's use of C routines.
2. the cost of copying vs. reference creation.

Regarding 2., the construction of the arma::sp_mat uses a copy since it does not have a reference (e.g. &) as a suffix. In particular, note:

#define ARMA_64BIT_WORD
#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::vec sp_multiply_copy(arma::sp_mat A, arma::vec nodes_status) {
  return A * nodes_status;
}

// [[Rcpp::export]]
arma::vec sp_multiply_ref(const arma::sp_mat& A, const arma::vec& nodes_status) {
  return A * nodes_status;
}

From here, the slight a slight performance difference between the two is noted:

benchmark(sp_multiply_copy(A, nodes.status),
          sp_multiply_ref(A, nodes.status),
          multiplyR(A, nodes.status), order = "relative")[, 1:4]

#                                test replications elapsed relative
# 3        multiplyR(A, nodes.status)          100   1.240    1.000
# 2  sp_multiply_ref(A, nodes.status)          100   2.766    2.231
# 1 sp_multiply_copy(A, nodes.status)          100   3.141    2.533

With this being said, we return to the first point: the Matrix functions for sparse matrices are highly optimized and directly use C. Examples of said routines can be viewed at R/products.R, src/Csparse.c, src/dtrMatrix.c. As a result, Matrix operations will be far more efficient.

Now, that is not to say speed cannot be obtained in C++. In particular, if the matrix object is repetitively used in C++ for multiplication with an instantiation of pass-by-reference (e.g. &), then it should be faster than a call to Matrix's multiplication routine.

---

Related Topics