我正在 RcppEigen 中编写一个用于加权协方差的函数。在其中一个步骤中,我想获取矩阵 X 的第 i 列和第 j 列,并计算应该返回某种向量的 cwiseProduct。cwiseProduct 的输出将进入一个可以多次重复使用的中间变量。从文档看来,cwiseProduct返回一个CwiseBinaryOp,它本身有两种类型。我的 cwiseProduct 对两个列向量进行操作,所以我认为正确的返回类型应该是Eigen::CwiseBinaryOp<Eigen::ColXpr, Eigen::ColXpr>,但是我收到错误no member named ColXpr in namespace Eigen
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
Rcpp::List Crossprod_sparse(Eigen::MappedSparseMatrix<double> X, Eigen::Map<Eigen::MatrixXd> W) {
int K = W.cols();
int p = X.cols();
Rcpp::List crossprods(W.cols());
for (int i = 0; i < p; i++) {
for (int j = i; j < p; j++) {
Eigen::CwiseBinaryOp<Eigen::ColXpr, Eigen::ColXpr> prod = X.col(i).cwiseProduct(X.col(j));
for (int k = 0; k < K; k++) {
//double out = prod.dot(W.col(k));
}
}
}
return crossprods;
}
Run Code Online (Sandbox Code Playgroud)
我也试过保存到 SparseVector
Eigen::SparseVector<double> prod = X.col(i).cwiseProduct(X.col(j));
Run Code Online (Sandbox Code Playgroud)
以及计算,但根本不节省
X.col(i).cwiseProduct(X.col(j));
Run Code Online (Sandbox Code Playgroud)
如果我根本不保存产品,函数会很快返回,暗示 cwiseProduct 不是一个昂贵的函数。当我将它保存到 SparseVector 中时,该函数非常慢,让我认为 SparseVector 不是正确的返回类型,而 Eigen 正在做额外的工作以将其转换为该类型。
回想一下,Eigen 依赖于表达式模板,因此如果您不分配表达式,则该表达式本质上是无操作的。在您的情况下,将其分配给 aSparseVector是正确的做法。关于速度,请确保在编译器优化开启的情况下进行编译(如-O3)。
尽管如此,我相信有一种更快的方法来编写您的整体计算。例如,你确定所有的X.col(i).cwiseProduct(X.col(j))都是非空的吗?如果不是,则应重写第二个循环以仅迭代重叠列的稀疏集。循环也可以互换以利用有效的矩阵产品。