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

您可以尝试计算权重矩阵的 Cholesky 分解，将矩阵乘以该分解，然后计算 RcppEigen 文档中列出的叉积。一些使用 RcppEigen 的示例代码可以是

#include <RcppEigen.h>

using Eigen::MatrixXd;
using Eigen::VectorXd;

//[[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
MatrixXd weightedCovariance(MatrixXd & X, MatrixXd & W) {
  int p = X.cols(); //assuming each row is a unique observation
  MatrixXd L = W.llt().matrixL();
  MatrixXd XtWX = MatrixXd(p, p).setZero().selfadjointView<Eigen::Lower>().rankUpdate(X.transpose() * L);
  return(XtWX);
}

// [[Rcpp::export]]
MatrixXd diag_weightedCovariance(MatrixXd & X, VectorXd & W) {
  int p = X.cols(); //assuming each row is a unique observation
  VectorXd w = W.cwiseSqrt();
  MatrixXd XtWX = MatrixXd(p, p).setZero().selfadjointView<Eigen::Lower>().rankUpdate(X.transpose() * w.asDiagonal());
  return(XtWX);
}

Eigen 在后台做了很多优化，所以告诉它结果是对称的应该会加快速度。使用微基准检查 R 中的时序：

set.seed(23847) #for reproducibility
require(microbenchmark)

#Create R version of Cpp function
Rcpp::sourceCpp('weighted_covar.cpp')

#generate data
p <- 100
n <- 1000
X <- matrix(rnorm(p*n), nrow=n, ncol=p)
W <- diag(1, n, n)
w <- diag(W)

R_res   <- crossprod(chol(W) %*% X ) #general weighted covariance
R_res_diag <- crossprod(sqrt(w) * X ) #utilizing your optimization, if we know it's diagonal
Cpp_res <- weightedCovariance(X, W)
Cpp_res_diag <- diag_weightedCovariance(X, w)

#make sure all equal
all.equal(R_res, Cpp_res)
#[1] TRUE
all.equal(R_res, R_res_diag)
#[1] TRUE
all.equal(Cpp_res_diag, R_res_diag)
#[1] TRUE

#check timings
microbenchmark(crossprod(chol(W) %*% X ))
# Unit: milliseconds
#                     expr      min      lq     mean  median       uq      max neval
# crossprod(chol(W) %*% X) 251.6066 262.739 275.1719 268.615 276.4994 479.9318   100

microbenchmark(crossprod(sqrt(w) * X ))
# Unit: milliseconds
#                   expr      min       lq     mean   median       uq     max neval
# crossprod(sqrt(w) * X) 5.264319 5.394289 5.499552 5.430885 5.496387 6.42099   100

microbenchmark(weightedCovariance(X, W))
# Unit: milliseconds
#                     expr      min       lq     mean   median       uq      max neval
# weightedCovariance(X, W) 26.64534 27.84632 31.99341 29.44447 34.59631 51.39726   100

microbenchmark(diag_weightedCovariance(X, w), unit = "ms")
# Unit: milliseconds
#                          expr     min       lq      mean   median        uq      max neval
# diag_weightedCovariance(X, w) 0.67571 0.702567 0.7469946 0.713579 0.7405515 1.321888   100

我也没有在这个实现中使用你的稀疏结构，所以在考虑到这一点之后你可能会得到更快的速度。

Answer 2

一般来说，如果你的产品中有一个对角矩阵，你应该只传递对角系数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;
    }
};