Реализация RCPP быстрых неотрицательных наименьших квадратов? - PullRequest
Новая
       50

Реализация RCPP быстрых неотрицательных наименьших квадратов?

0 голосов
/ 19 сентября 2019

Я искал быструю реализацию в R алгоритма быстрого (на основе активного набора) неотрицательного наименьших квадратов Bro, R., & de Jong, S. (1997) Быстрый неотрицательность.алгоритм наименьших квадратов с ограничениями.Journal of Chemometrics, 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

  }

, но в моих тестах она намного медленнее, чем *Функция 1008 * plain nnls в пакете nnls (закодировано на фортране), хотя алгоритмически fnnls должна быть быстрее.Мне было интересно, есть ли у кого-нибудь доступный Rcpp порт fnnls, в идеале с использованием классов броненосцев и разрешением разреженности X, а также, возможно, поддержка Y для нескольких столбцов?

...