英文:
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 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?
答案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;
}
/*** R
v <- 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;
}
/*** 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.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论