时间戳

首页 > 解决方案 > 如何加快我的 Rcpp 代码，它只执行简单的操作？

r - 如何加快我的 Rcpp 代码，它只执行简单的操作？

问题描述

我正在尝试编写一个函数，该函数接受一个矩阵并为每对列计算一个值。矩阵总是有 2000 行，但可能有非常多的列（最多 100,000 左右）。我开始使用的 R 代码如下：

x_dist <- data.frame(array(0,dim=c(ncol(x),ncol(x))))
cs <- colSums(x)
for (i in 1:ncol(x)) {
  p_i <- x[,i]
  for (j in 1:ncol(x)) {
    p_j <- x[,j]
    s <- p_i+p_j
    fac <- cs[i]/(cs[i]+cs[j])
    N1 <- fac*s
    N2 <- (1-fac)*s
    d1 <- (p_i+1)/(N1+1)
    d2 <- (p_j+1)/(N2+1)
    x_dist[i,j] <- sum(N1+N2-N1*d1-N2*d2+p_i*log(d1)+p_j*log(d2))
  }
}

这个功能很慢。当矩阵只有 400 列时x，大约需要 32 秒，并且列数明显呈二次增长。

因为我听说 Rcpp 对加速循环和矩阵运算很有好处，所以我决定试一试。我对它完全陌生，但最终将以下功能放在一起：

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix wdist(NumericMatrix x) {
    int nrow = x.nrow(),ncol=x.ncol();
    NumericMatrix m = no_init_matrix(ncol,ncol);
    NumericVector v1 = no_init_vector(nrow);
    NumericVector v2 = no_init_vector(nrow);
    NumericVector s = no_init_vector(nrow);
    NumericVector N1 = no_init_vector(nrow);
    NumericVector N2 = no_init_vector(nrow);
    NumericVector d1 = no_init_vector(nrow);
    NumericVector d2 = no_init_vector(nrow);
    for(int i=0; i<ncol; ++i){
      v1 = x(_,i);
      for(int j=0; j<i; ++j){
        v2 = x(_,j);
        s = v1+v2;
        N1 = sum(v1)*s/(sum(v1)+sum(v2));
        N2 = s-N1;
        d1 = (v1+1)/(N1+1);
        d2 = (v2+1)/(N2+1);
        m(i,j) = sum(N1+N2-N1*d1-N2*d2+v1*log(d1)+v2*log(d2));
      }
    }
    return m;
}

这肯定会产生很大的不同。现在有 400 列，这大约需要 8 秒。我对改进感到高兴，但是对于我当前感兴趣的测试用例（32,000 列）来说，这仍然非常缓慢。我觉得我在做一些相对简单的操作，所以我很困惑为什么我的代码仍然这么慢。我尝试阅读有关编写高效 Rcpp 代码的内容，但没有找到任何有助于解决我的问题的内容。请让我知道我是否做错了什么，或者我可以进行任何改进以使我的代码更快（甚至是 R 代码本身，如果它可以比 Rcpp 代码更快的话！）

一些示例数据可能是：

set.seed(121220) 
x <- array(rpois(2000*400,3),dim=c(2000,400))

标签： rdataframercpp

解决方案

我重构了您的基本 R 代码，希望它可以加快速度

f <- function(...) {
  p <- x[, t(...)]
  N <- matrix(rowSums(p), ncol = 1) %*% colSums(p) / sum(p)
  d <- (p + 1) / (N + 1)
  sum(N - N * d + p * log(d))
}
x_dist <- diag(0, ncol(x))
x_dist[lower.tri(x_dist)] <- combn(ncol(x), 2, FUN = f)
x_dist <- pmax(x_dist, t(x_dist))

为了加快您的 Rcpp 代码，您可以在将矩阵初始化为全零矩阵for后尝试以下嵌套循环：m

    for(int i=0; i<ncol-1; ++i){
      v1 = x(_,i);
      for(int j=i+1; j<ncol; ++j){
        v2 = x(_,j);
        s = v1+v2;
        N1 = sum(v1)*s/sum(s);
        N2 = s-N1;
        d1 = (v1+1)/(N1+1);
        d2 = (v2+1)/(N2+1);
        val = sum(N1+N2-N1*d1-N2*d2+v1*log(d1)+v2*log(d2));
        m(i,j) = val;
        m(j,i) = val;
      }
    }

它应用了矩阵是对称的属性，从而将计算复杂度降低了一半。

推荐阅读