英文:
Rcpp Eigen inplace matrix multiplication error
问题
我是新手学习Rcpp和Eigen，尝试运行一个原地矩阵乘法时，遇到了这个错误。R没有告诉我错误的具体信息。
SEXP cpp_hom_crit(const Eigen::Map<Eigen::MatrixXd> Q, Eigen::Map<Eigen::VectorXd> d, Eigen::Map<Eigen::MatrixXd> Qinv) {Q *= d.asDiagonal() * Qinv; //错误发生在这里return Rcpp::wrap(Q);}
我也尝试了这个，但也不起作用
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; //错误发生在这里return Rcpp::wrap(Q);}
如何解决这个错误，可能让它运行更快？
英文:
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 herereturn 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 herereturn Rcpp::wrap(Q);}
How can I resolve this error and possibly make this faster?
答案1
得分: 2
以下是已翻译的部分：
这是解决您的问题的一种方法。正如评论中@chtz所提到的，您不能声明一个变量为const然后赋值给它。因此，在这里我们赋值给一个新的临时变量，然后返回该结果。
代码
#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;}/*** Rv <- as.numeric(1:3)M <- v %*% t(v)Mi <- chol2inv(M)cpp_hom_crit(M, v, Mi)*/
输出
> 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>
这些值当然是毫无意义的，但它们至少说明了正确的维度。请注意，像LAPACK这样的库通常有专门用于常见的矢量-矩阵乘积表达式的函数，内部进行循环操作，可能还会避免临时变量，因此不清楚Eigen是否会自动更快（多）。
英文:
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;}/*** Rv <- 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.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论