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

Question

我正在寻找 R 中的快速（基于活动集的）非负最小二乘算法的快速实现 Bro, R., & de Jong, S. (1997) A fast non-negativity-constrained least squares algorithm. Journal of Chemometrics, 11, 393-401. 在 multiway package I found this pure R implementation:

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

  }

但在我的测试中它比 plain nnls function in the nnls package（用 fortran 编码）慢得多，尽管算法 fnnls 应该更快。我想知道是否有人碰巧有 fnnls 的 Rcpp 端口可用，理想情况下使用犰狳类并允许 X 稀疏并且可能还支持 Y 有多个列？

Answer 1

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

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

我也花了将近两天的时间来尝试解析 multiway::fnnls 实现，并将避免在 R 礼节、可解释性和内存使用方面使用选择词。

我不明白为什么 multiway::fnnls 的作者认为他们的实现应该很快。考虑到 fortran Lawson/Hanson 实现，R-only 实现似乎毫无用处。

这是我写的一个 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 版本时，我将 post 启动。

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

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

Rcpp implementation of fast nonnegative least squares?

r

constraints

rcpp

least-squares