RcppEigen中的有效加权协方差

Question

我正在尝试生成一个可以计算一系列加权产品的函数

其中W是对角矩阵.有许多W矩阵但只有一个X矩阵.

为了有效,我可以将W表示为包含对角线部分的数组(w).然后在R中这将是 crossprod(X, w*X)

要不就 crossprod(X * sqrt(w))

我可以循环W系列,但这似乎效率低下.整个产品可以作为只有w改变才能回收第i列和第j列的产品X_i*X_j.我想要制作的功能看起来像这样

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 k = 0; k < K; k++) {
    Eigen::SparseMatrix<double> matprod(p, p);
    for (int i = 0; i < p; i++) {
      Eigen::SparseVector<double> prod = X.col(i).cwiseProduct(W.col(k));
      for (int j = i; j < p; j++) {
        double out = prod.dot(X.col(j));
        matprod.coeffRef(i,j) = out;
        matprod.coeffRef(j,i) = out;
      }
    }
    matprod.makeCompressed();
    crossprods[k] = matprod;
  }

  return crossprods;
}

它返回正确的产品,并且因为对中间prod变量进行操作而应该是高效的.然而,crossprod尽管没有利用回收利用,但使用R的循环似乎仍然要快得多.如何更好地优化此功能？

Answer 1

一般来说，如果产品中有对角矩阵，则应该仅传递对角系数w并将其用作w.asDiagonal()：

Eigen::MatrixXd foo(Eigen::SparseMatrix<double> const & X, Eigen::VectorXd const & w)
{
    return X.transpose() * w.asDiagonal() * X;
}

如果您想预先计算除与 相乘之外的所有内容w，您可以尝试存储每行的外积X并按需累加：

class ProductHelper
{
    std::vector<Eigen::SparseMatrix<double> > matrices;
public:
    ProductHelper(Eigen::SparseMatrix<double> const& X_)
    {
        // The loop below is much more efficient with row-major X
        Eigen::SparseMatrix<double, Eigen::RowMajor> const &X = X_;
        matrices.reserve(X.rows());
        for(int i=0; i<X.rows(); ++i)
        {
            matrices.push_back(X.row(i).transpose()*X.row(i));
        }
    }

    Eigen::MatrixXd multiply(Eigen::VectorXd const& w) const
    {
        assert(w.size()==matrices.size());
        assert(w.size()>0);
        Eigen::MatrixXd A = w[0]*matrices[0]; 
        for(int i=1; i<w.size(); ++i)
        {
            A+=w[i]*matrices[i];
        }
        return A;
    }
};