**2022年1月编辑：**使用CRAN上RcppML R包中的nnls函数。这是下面给出的函数的更快的基于特征的实现，随后是坐标下降最小二乘改进。

为了研究这个问题，我花了将近一个星期的时间。
我还花了将近两天的时间试图解析multiway::fnnls实现，并且将避免在R礼仪、可解释性和内存使用方面使用选择词。
我不明白为什么multiway::fnnls的作者认为他们的实现应该很快，考虑到fortranslawson/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代码并将其运行到Armadillo中，但失败了。
我正在努力实现这个方法作为一个Rcpp包，并将张贴时，它使它成为一个稳定的Github版本。
附注：使用Eigen可以使这更快。Armadillo求解器使用Cholesky分解和直接替换。Eigen的Cholesky求解器更快，因为它在原地做了更多的操作。

我要求OpenAI的代码-davinci-002 https://beta.openai.com/playground?model=code-davinci-002将multiway包中的fnnls纯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;  
}

型