Эффективные взвешенные ковариации в RcppEigen

Я пытаюсь создать функцию, которая может вычислять серию взвешенных продуктов

где W - диагональная матрица. Существует много W-матриц, но только одна X-матрица.

Чтобы быть эффективным, я могу представить W как массив (w), содержащий диагональную часть. Тогда в R это будет crossprod(X, w*X)

или просто crossprod(X * sqrt(w))

Я мог бы перебрать цикл W, но это кажется неэффективным. Весь продукт может быть как Изменяется только w, поэтому продукты X_i * X_j для столбцов я и 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. Однако для циклирования в R с использованием crossprod, похоже, все еще намного быстрее, несмотря на то, что не использует переработку. Как я могу оптимизировать эту функцию больше?

Ответ 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;
    }
};

Ответ 2

Вы можете попробовать вычислить декомпозицию 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);
}

Эйген проводит большую оптимизацию под капотом, поэтому сообщение о том, что результат симметричен, должно ускорить процесс. Проверка времени в 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 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

Я также не использовал вашу разреженную структуру в этой реализации, поэтому вы можете получить большую скорость после учета этого.