用 C++ 编写的极其低效的代码

Question

以下函数的计算时间非常长。有什么改进的余地吗？我应该以不同的方式访问矩阵 X 的元素吗？我很感激任何意见或建议。

#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]


arma::mat myfunc(const int& n,
                 const int& p,
                 arma::mat& X,
                 arma::rowvec& y,
                 const arma::rowvec& types,
                 const arma::mat& rat,
                 const arma::rowvec& betas){
  
  arma::mat final(p+p,n);
  final.zeros();
  int i,j;
  
  for(i=0; i < n; ++i){
    arma::colvec finalfirst(p+p); finalfirst.zeros();
    for(j=0; j < n; ++j){
      arma::mat Xt = X * log(y(j));
      arma::mat finalX = join_rows(X,Xt);
      
      arma::rowvec Xi = finalX.row(i);
      
      if(types(i)==1 && y(j)==y(i)){
        finalfirst += (Xi.t() - rat.col(j));
      }
      if(types(i)>1 && y(j) > y(i)){
        finalfirst -= (Xi.t() - rat.col(j)) * exp(arma::as_scalar(betas*Xi.t()));
        
      }
      else if(y(j) <= y(i)){
        finalfirst -= Xi.t() * exp(arma::as_scalar(betas*Xi.t()));
      }
    }
    
    final.col(i) = finalfirst;
  }
  
  return(final);
}




/*** R
m=4000
types = runif(m,0,5)
types[types<=1] = 0 ; types[types > 1 & types < 3] = 1; types[types>=3]=2
microbenchmark(out = myfunc(n=m,p=2,X=matrix(rnorm(m*2),nrow=m,ncol=2),y=runif(m,0,3),types=types,rat=matrix(rnorm(m*4),nrow=4,ncol=m),betas=c(1,2,3,4)))
*/

Answer 1

如果您至少一些解释您正在做什么，这可能会有所帮助。

首先，我建议将代码分解成非常小的块，对每个块进行微基准测试，并找出瓶颈操作。这比一下子把整个功能都扔到墙上要好。

行与列访问

join_rows(X,Xt)

这是一个非常慢的操作，特别是因为您要将 rowvec 添加到面向列的 mat。 Armadillo 矩阵以列优先顺序存储为向量，因此在幕后 Armadillo 在 n 非连续位置调用 push_back，其中 n 是列数。

看来您可以完全避免这种情况，因为 finalX.row(i) 是此循环中唯一依赖于 finalX 的调用。只要弄清楚 .row(i) 是什么。

隐式操作

你做了很多转置，也许你应该从一开始就使用转置矩阵？

Xi.t() 在同一行被调用了两次： finalfirst -= (Xi.t() - rat.col(j)) * exp(arma::as_scalar(betas*Xi.t()));

转置行向量毫无意义，只需将其初始化为列向量并使用老式循环遍历它即可。多一点代码并不总是坏事，它通常也会使您的意图更加明确。

复制与就地操作

这是一个副本： final.col(i) = finalfirst;

为什么不在内存中操作并就地更新 final(_, i) 而不是使用临时 finalfirst 然后将深度复制到目标矩阵中？这是一个列，因此内存是连续的，编译器将能够为您优化访问模式，就像您在处理任何简单向量一样。

总而言之，我还没有完全理解你在做什么，但是@largest_prime_is_463035818 切换 for(i...) 和 [=24= 可能是正确的] 循环。您似乎可以将这两行从嵌套循环中拉出：

arma::mat Xt = X * log(y(j));
arma::mat finalX = join_rows(X,Xt);

因为他们根本不依赖 i。

用 C++ 编写的极其低效的代码

Extremely inefficient code written in C++

c++

r

rcpp

rcpparmadillo

行与列访问

隐式操作

复制与就地操作