快速非负最小二乘的Rcpp实现？

问题描述 投票：0回答：2

我正在寻找 R 中的快速（基于活动集）非负最小二乘算法的快速实现 Bro, R., & de Jong, S. (1997) 一种快速非负约束最小二乘算法。化学计量学杂志，11, 393-401。 在 multiway 包中我发现了这个纯 R 实现:

fnnls <- 
  function(XtX,Xty,ntol=NULL){     
    ### initialize variables
    pts <- length(Xty)
    if(is.null(ntol)){
      ntol <- 10*(.Machine$double.eps)*max(colSums(abs(XtX)))*pts
    }
    pvec <- matrix(0,1,pts)
    Zvec <- matrix(1:pts,pts,1)
    beta <- zvec <- t(pvec)
    zz <- Zvec
    wvec <- Xty - XtX%*%beta
    
    ### iterative procedure
    iter <- 0    
    itmax <- 30*pts
    
    # outer loop
    while(any(Zvec>0) && any(wvec[zz]>ntol)) {
      
      tt <- zz[which.max(wvec[zz])]
      pvec[1,tt] <- tt
      Zvec[tt] <- 0
      pp <- which(pvec>0)
      zz <- which(Zvec>0)
      nzz <- length(zz)
      zvec[pp] <- smpower(XtX[pp,pp],-1)%*%Xty[pp]
      zvec[zz] <- matrix(0,nzz,1)
      
      # inner loop
      while(any(zvec[pp]<=ntol) &&  iter<itmax) {
        
        iter <- iter + 1
        qq <- which((zvec<=ntol) & t(pvec>0))
        alpha <- min(beta[qq]/(beta[qq]-zvec[qq]))
        beta <- beta + alpha*(zvec-beta)
        indx <- which((abs(beta)<ntol) & t(pvec!=0))
        Zvec[indx] <- t(indx)
        pvec[indx] <- matrix(0,1,length(indx))
        pp <- which(pvec>0)
        zz <- which(Zvec>0)
        nzz <- length(zz)
        if(length(pp)>0){
          zvec[pp] <- smpower(XtX[pp,pp],-1)%*%Xty[pp]
        }
        zvec[zz] <- matrix(0,nzz,1)      
        
      } # end inner loop
      
      beta <- zvec
      wvec <- Xty - XtX%*%beta
      
    } # end outer loop
    
    beta
    
  }

但在我的测试中，它比 nnls

 包
（用 Fortran 编码）中的 
plain
nnls 函数慢得多，尽管算法上 
fnnls
应该更快。我想知道是否有人碰巧有 
Rcpp
的 
fnnls
端口可用，理想情况下使用犰狳类并允许 
X
稀疏，也许还支持 Y 具有多个列？或者其他一些适用于稀疏协变量矩阵和多个右侧的 nnls 实现？

r constraints rcpp least-squares
2个回答
4
投票

编辑 2022 年 1 月： 使用 CRAN 上 RcppML R 包中的 

nnls
函数。这是下面给出的函数的更快的基于特征的实现，随后是坐标下降最小二乘细化。

出于研究目的，我花了将近一周的时间来解决这个问题。

我还花了近两天的时间尝试解析 

multiway::fnnls
实现，并且不会在 R 礼仪、可解释性和内存使用方面使用选择词。

我不明白为什么

multiway::fnnls
的作者认为他们的实现应该很快。考虑到 fortran Lawson/Hanson 实现，仅 R 实现似乎毫无用处。

这是我编写的 RcppArmadillo 函数（快速近似解轨迹）NNLS，它为条件良好的系统复制 

multiway::fnnls
：

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

using namespace arma;
typedef unsigned int uint;

// [[Rcpp::export]]
vec fastnnls(mat a, vec b) {

  // initial x is the unbounded least squares solution
  vec x = arma::solve(a, b, arma::solve_opts::likely_sympd + arma::solve_opts::fast);
  
  while (any(x < 0)) {

    // define the feasible set as all values greater than 0
    arma::uvec nz = find(x > 0);
    
    // reset x
    x.zeros();
    
    // solve the least squares solution for values in the feasible set
    x.elem(nz) = solve(a.submat(nz, nz), b.elem(nz), arma::solve_opts::likely_sympd + arma::solve_opts::fast);
  }
  return x;
}

这种方法本质上是 TNT-NN 的前半部分，但在每次迭代时没有“启发式”尝试从可行集中添加或删除元素。

为了使这种方法超越简单的近似，我们可以添加顺序坐标下降，它接收上面的 

FAST
解作为初始化。一般来说，对于大多数条件良好的小问题，99% 的情况下，
FAST
都会给出精确的解决方案。

上述实现的一个独特属性是它不能给出误报，但有时（在大型或病态系统中）会给出误报。因此，可能比实际解决方案稍微稀疏。请注意，FAST 和精确解之间的损失通常在 1% 以内，因此，如果您不追求绝对精确解，那么这是您最好的选择。

上述算法的运行速度也比 Lawson/Hanson nnls 求解器快得多。这是我刚刚从 50 系数系统复制的微基准，复制了 10000 次：

Unit: microseconds
               expr   min    lq      mean median    uq     max neval
           fastnnls  53.9  56.2  59.32761   58.0  59.5   359.7 10000
 lawson/hanson nnls 112.9 116.7 125.96169  118.6 129.5 11032.4 10000

当然，性能会根据密度和负性而变化。与其他算法相比，我的算法往往会随着稀疏性的增加而变得更快，并且随着正解的减少而变得更快。

我尝试过简化 multiway::fnnls 代码并将其运行到犰狳中，但失败了。

我正在努力将此方法实现为 Rcpp 包，并将在它成为稳定的 Github 版本时发布。

p.s.：使用 Eigen 可以加快速度。犰狳求解器使用 Cholesky 分解和直接替换。 Eigen 的 Cholesky 求解器速度更快，因为它可以就地执行更多操作。

0
投票

为了它的价值，我要求 OpenAI 的 code-davinci-002 https://beta.openai.com/playground?model=code-davinci-002 将 

fnnls
包中的 
multiway
纯 R 代码翻译为
Rcpp
，这就是它的结果 - 虽然可能不是最佳的，但也许是任何对此感兴趣的人的开始......

/* Convert the following R code below to calculate a fast nonnegative least squares fit to Rcpp and use RcppArmadillo classes. */

fnnls <- 
  function(XtX,Xty,ntol=NULL){
    # Fast Non-Negative Least Squares algorithm based on 
    #   Bro, R., & de Jong, S. (1997) A fast non-negativity-constrained 
    #   least squares algorithm. Journal of Chemometrics, 11, 393-401.
    # Nathaniel E. Helwig ([email protected])
    # last updated: April 9, 2015
    
    ### initialize variables
    pts <- length(Xty)
    if(is.null(ntol)){
      ntol <- 10*(.Machine$double.eps)*max(colSums(abs(XtX)))*pts
    }
    pvec <- matrix(0,1,pts)
    Zvec <- matrix(1:pts,pts,1)
    beta <- zvec <- t(pvec)
    zz <- Zvec
    wvec <- Xty - XtX%*%beta
    
    ### iterative procedure
    iter <- 0    
    itmax <- 30*pts
    
    # outer loop
    while(any(Zvec>0) && any(wvec[zz]>ntol)) {
      
      tt <- zz[which.max(wvec[zz])]
      pvec[1,tt] <- tt
      Zvec[tt] <- 0
      pp <- which(pvec>0)
      zz <- which(Zvec>0)
      nzz <- length(zz)
      zvec[pp] <- smpower(XtX[pp,pp],-1)%*%Xty[pp]
      zvec[zz] <- matrix(0,nzz,1)
      
      # inner loop
      while(any(zvec[pp]<=ntol) &&  iter<itmax) {
        
        iter <- iter + 1
        qq <- which((zvec<=ntol) & t(pvec>0))
        alpha <- min(beta[qq]/(beta[qq]-zvec[qq]))
        beta <- beta + alpha*(zvec-beta)
        indx <- which((abs(beta)<ntol) & t(pvec!=0))
        Zvec[indx] <- t(indx)
        pvec[indx] <- matrix(0,1,length(indx))
        pp <- which(pvec>0)
        zz <- which(Zvec>0)
        nzz <- length(zz)
        if(length(pp)>0){
          zvec[pp] <- smpower(XtX[pp,pp],-1)%*%Xty[pp]
        }
        zvec[zz] <- matrix(0,nzz,1)      
        
      } # end inner loop
      
      beta <- zvec
      wvec <- Xty - XtX%*%beta
      
    } # end outer loop
    
    beta
    
  }

答案：

# [[Rcpp::export]]
arma::vec fnnls(const arma::mat& XtX, const arma::vec& Xty, double ntol = NULL){
  
  // initialize variables
  int pts = Xty.n_elem;
  if (ntol == NULL){
    ntol = 10*(std::numeric_limits<double>::epsilon())*arma::max(arma::sum(arma::abs(XtX),1))*pts;
  }
  arma::mat pvec = arma::zeros<arma::mat>(1,pts);
  arma::vec Zvec = arma::regspace<arma::vec>(1,pts);
  arma::vec beta = arma::zeros<arma::vec>(pts);
  arma::vec zvec = pvec.t();
  arma::uvec zz = find(Zvec>0);
  arma::vec wvec = Xty - XtX*beta;

  // iterative procedure
  int iter = 0;    
  int itmax = 30*pts;

  // outer loop
  while((any(Zvec>0)) && (any(wvec(zz)>ntol))){

    int tt = zz[wvec.elem(zz).index_max()];
    pvec.col(tt) += tt;
    Zvec[tt] = 0;
    arma::uvec pp = find(pvec>0);
    zz = find(Zvec>0);
    int nzz = zz.n_elem;

    zvec.elem(pp) = solve(XtX.submat(pp,pp),Xty.elem(pp));

    zvec.elem(zz).zeros();

    // inner loop
    while((any(zvec.elem(pp)<=ntol)) && (iter<itmax)){

      iter++;

      arma::uvec qq = find((zvec<=ntol) && (pvec>0));

      double alpha = min((beta.elem(qq)/((beta.elem(qq)-zvec.elem(qq)))));

      beta += alpha*((zvec-beta));

      arma::uword indx;

      for (arma::uword i=0; i < beta.n_elem; i++){

        if ((abs(beta[i]) < ntol) && (pvec[i] != 0)){

          indx = i;

          Zvec[indx] += indx;

          pvec[indx] += 0;

        }

      }

      pp = find(pvec>0);

      zz = find(Zvec>0);

      nzz = zz.n_elem;

      if (pp.n_elem > 0){

        zvec.elem(pp) = solve(XtX.submat(pp,pp),Xty.elem(pp));

      }      

      zvec.elem(zz).zeros();      

    } // end inner loop    

    beta = zvec;    

    wvec -= XtX*beta;    

  } // end outer loop  

  return beta;  
}
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.