犰狳中的有效加权叉积

Question

我正在尝试找到一种计算 X^T * W * X 的有效方法，其中 X 是一个密集的 mat 大小，例如10,000 x 10 并且 W 是对角矩阵（我只将对角线存储在 vec 中）。

目前，我使用这个功能

arma::mat& getXtW(const arma::mat& covar,
                  const arma::vec& w,
                  arma::mat& tcovar,
                  size_t n, size_t K) {
  size_t i, k;

  for (i = 0; i < n; i++) {
    for (k = 0; k < K; k++) {
      tcovar(k, i) = covar(i, k) * w(i);
    }
  }

  return tcovar;
}

并计算

tcovar = getXtW(covar, w, tcovar, n, K);
cprod = tcovar * covar;

然而，这似乎不是最优的。

PS：你可以看到整个代码there。

Edit1：似乎我可以使用covar.t() * (covar.each_col() % w)，但这似乎并没有快多少。

Edit2：如果我自己在 Rcpp 中使用循环实现它：

arma::mat testProdW2(const arma::mat& x, const arma::vec& w) {

  int n = x.n_rows;
  int K = x.n_cols;
  arma::mat res(K, K);
  double tmp;
  for (int k = 0; k < K; k++) {
    for (int j = k; j < K; j++) {
      tmp = 0;
      for (int i = 0; i < n; i++) {
        tmp += x(i, j) * w[i] * x(i, k);
      }
      res(j, k) = tmp;
    }
  }

  for (int k = 0; k < K; k++) {
    for (int j = 0; j < k; j++) {
      res(j, k) = res(k, j);
    }
  }

  return res;
}

这比第一个实现要慢。

Answer 1

根据BLAS matrix by matrix transpose multiply，没有直接执行此操作的 BLAS 例程。相反，建议遍历X 的行并使用dsyr。我发现这是一个有趣的问题，因为我知道如何link BLAS in Rcpp，但还没有使用 RcppArmadillo 这样做。 Stack Overflow 也知道答案：Rcpparmadillo: can't call Fortran routine "dgebal"?。注意：我尚未检查，但我希望 dsyr 不是 R 附带的 BLAS 子集的一部分。因此，只有当您的 R 链接到完整的 BLAS 实现时，这才有效。

结合我们得到：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp/Benchmark/Timer.h>

#ifdef ARMA_USE_LAPACK
#if !defined(ARMA_BLAS_CAPITALS)
 #define arma_dsyr dsyr
#else
 #define arma_dsyr DSYR
#endif

extern "C"
void arma_fortran(arma_dsyr)(const char* uplo, const int* n, const double* alpha, const double* X, const int* incX, const double* A, const int* ldA);

#endif


// [[Rcpp::export]]
Rcpp::NumericVector getXtWX(const arma::mat& X, const arma::vec& w) {
  Rcpp::Timer timer;
  timer.step("start"); 

  arma::mat result1 = X.t() * (X.each_col() % w);
  timer.step("Armadillo result");

  const int n = X.n_rows;
  const int k = X.n_cols;
  arma::mat result(k, k, arma::fill::zeros);  
  for (size_t i = 0; i < n; ++i) {
    F77_CALL(dsyr)("U", &k, &w(i), &X(i,0), &n, result.memptr(), &k);
  }
  result = arma::symmatu(result);
  timer.step("BLAS result");
  Rcpp::NumericVector res(timer);
  return res;
}

/*** R
n <- 10000
k <- 10
X <- matrix(runif(n*k), n, k)
w <- runif(n)
Reduce(rbind, lapply(1:6, function(x) diff(getXtWX(X, w))/1e6))
*/

但是，对我来说，BLAS 解决方案要慢一些：

> Reduce(rbind, lapply(1:6, function(x) diff(getXtWX(X, w))/1e6))
     Armadillo result BLAS result
init         1.291243    6.666026
             1.176143    6.623282
             1.102111    6.644165
             1.094917    6.612596
             1.098619    6.588431
             1.069286    6.615529

我试图通过首先转置矩阵来改进这一点，希望在循环列矩阵时内存访问会更快，但这对我的（低功耗）系统没有影响。