Using the function kronecker() crashes my R session

[geyh96]

Adding the line contains the kronecker() function crashes my R session. How can I use the functions in the unsupported module of the Eigen in Rcpp?

#include "math.h"
#include <Rcpp.h>
#include <RcppEigen.h>

#include <omp.h>
using namespace Rcpp;
using namespace std;
using namespace Eigen;
//[[Rcpp::depends(RcppEigen)]]
//[[Rcpp::plugins(openmp)]]

//[[Rcpp::export]]
Eigen::MatrixXd myFun(Eigen::Map<Eigen::MatrixXd> alpha, Eigen::Map<Eigen::MatrixXd> Xdata, Eigen::Map<Eigen::MatrixXd> Kmatrix, double sigma, int nparal)
{
  Eigen::MatrixXd Alpha, X,K;
  Alpha=alpha;
  X=Xdata;
  K=Kmatrix;
  int n = X.rows();
  int p = X.cols();
  int k =0;
  int i;
  Eigen::MatrixXd tempA;
  Eigen::MatrixXd In;

  tempA = Eigen::kroneckerProduct(X,X);

 return X;
}

[eddelbuettel]

See https://eigen.tuxfamily.org/dox/unsupported/group__KroneckerProduct__Module.html and in particular try this:

Warning If you want to replace a matrix by its Kronecker product with some matrix, do NOT do this: A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect instead, use eval() to work around this: A = kroneckerProduct(A,B).eval();

[geyh96]

Thank you for your help. It is not the cause of my problem. I tried the eval() before. I rewrite the function and it works and I don't know why. If it happens again, I will update the issue. Thanks for your help and your contribution to the R community! @eddelbuettel

[eddelbuettel]

I recommend you try a self-contained five-liner in C++ to multiply two 2x2 or 3x3 matrices. Then do the same from RcppEigen. Only then do it with an argument from R. You should be able to drill down to the actual issue, and I suspect it will have to do more with Eigen than with the interface to it we provide.

[eddelbuettel]

There is a test file unsupported/test/kronecker_product.cpp in the Eigen upstream sources with a simple use case calling kroneckerProduct() at the end. Adapting that in a minimal and verifiable example (as opposed to your example including OpenMP and what not):

Code

#include <RcppEigen.h>

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

//[[Rcpp::export]]
Eigen::MatrixXd myFun(Eigen::Map<Eigen::MatrixXd> a, Eigen::Map<Eigen::MatrixXd> b) {
  Eigen::MatrixXd res = kroneckerProduct(a,b);
  return res;
}

/*** R
a <- matrix(as.numeric(1:4), 2, 2)
b <- matrix(as.numeric(4:1), 2, 2)
myFun(a,b)
*/

Output

> Rcpp::sourceCpp("/tmp/kron.cpp")
[...bazillion line of compiler noise under -Wignored-attributes and other flags omitted ...]
> a <- matrix(as.numeric(1:4), 2, 2)

> b <- matrix(as.numeric(4:1), 2, 2)

> myFun(a,b)
     [,1] [,2] [,3] [,4]
[1,]    4    2   12    6
[2,]    3    1    9    3
[3,]    8    4   16    8
[4,]    6    2   12    4
> 

So I am closing this. Please do feel free to reopen if you find an issue in RcppEigen related to this.