Rcpp Eigen inplace matrix multiplication error

Question

I'm new to Rcpp and Eigen and I'm trying to run a inplace matrix multiplication while I run into this error. R didn't tell me what the error is about.

SEXP cpp_hom_crit(
    const Eigen::Map<Eigen::MatrixXd> Q, Eigen::Map<Eigen::VectorXd> d, Eigen::Map<Eigen::MatrixXd> Qinv) {

    Q *= d.asDiagonal() * Qinv; //error happens here

    return Rcpp::wrap(Q);
}

I tried this one also, but it does not work either

SEXP cpp_hom_crit(
    const Eigen::Map<Eigen::MatrixXd> Q, Eigen::Map<Eigen::VectorXd> d, Eigen::Map<Eigen::MatrixXd> Qinv) {

    Q = Q * d.asDiagonal() * Qinv; //error happens here

    return Rcpp::wrap(Q);
}

How can I resolve this error and possibly make this faster?

Remove the `const` from `const Eigen::Map Q`. (The other arguments could have a `const` added or even be changed to `const Eigen::Map` or `const Eigen::Ref&` (note the extra `&` in `const Ref&`). — chtz, Jul 06 '23 at 07:16

score 2 · Accepted Answer · answered Jul 06 '23 at 13:42

Here is one way to fix your problem. As @chtz alluded to in the comment, you cannot declare a variable const and assign to it. So here we assign to a new temporary variable and return that result.

Code

#include <RcppEigen.h>

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

// [[Rcpp::export]]
Eigen::MatrixXd cpp_hom_crit(const Eigen::Map<Eigen::MatrixXd> Q,
                             const Eigen::Map<Eigen::VectorXd> d,
                             const Eigen::Map<Eigen::MatrixXd> Qinv) {

    auto res = Q * d.asDiagonal() * Qinv;
    return res;
}

/*** R
v <- as.numeric(1:3)
M <- v %*% t(v)
Mi <- chol2inv(M)
cpp_hom_crit(M, v, Mi)
*/

Output

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

> v <- as.numeric(1:3)

> M <- v %*% t(v)

> Mi <- chol2inv(M)

> cpp_hom_crit(M, v, Mi)
     [,1]      [,2]      [,3]
[1,] 0.75 0.0694444 0.0370370
[2,] 1.50 0.1388889 0.0740741
[3,] 2.25 0.2083333 0.1111111
>

The values are of course nonsensical but they illustrate at least the correct dimension. Note that libraries like LAPACK often have dedicated functions for common vector-matrix product expressions, do the looping internally and may also avoid the temporary so not it is not clear that Eigen will automatically be (much) faster.

Rcpp Eigen inplace matrix multiplication error

1 Answers1

Code

Output