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?

Answer 1

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.