[R] Summary: non-negative least squares

[Robert Abugov]

Thank you Brian Ripley, Gardar Johannesson, and Marcel Wolbers for your
prompt
and friendly help! I will share any further learnings as I move through
these suggestions.    -Bob Abugov


Brian Ripley wrote:

I just use optim() on the sum of squares with non-negativity constraints.
That did not exist in 1999.


Gardar Johannesson wrote:

You can always just use the quadratic programing library in R (quadprog).
That is, if you have

y as your data vector
X as your design matrix

Then do,

D <- t(X) %*% X
d <- t(t(y) %*% X)
A <- diag(ncol(X))
b <- rep(0,ncol(X))
fit <- solve.QP(D=D,d=d,A=t(A),b=b,meq=0)

and the solution is in fit$solution

Good luck,
Gardar


Marcel Wolbers wrote:

Hi Bob

I use the routine below which uses library quadprog for
non-negative least squares.
The scaling (xscale, yscale) is in principle unnecessary; it's
only there because there's a minor bug in the library quadprog
and the scaling is often a workaround.
I don't know how this compares to Brian Ripleys suggestion (use
optim) in times of speed and accuracy, but it should be competitive.

Yours,
Marcel


#---------------------------------------------------------------------------
--
nnls.fit <- function(x,y,wsqrt=1,eps=0,rank.tol=1e-07) {
  ## Purpose: Nonnegative Least Squares (similar to the S-Plus function
  ## with the same name) with the help of the R-library quadprog

 ## ------------------------------------------------------------------------
  ## Attention:
  ## - weights are square roots of usual weights
  ## - the constraint is coefficient>=eps

 ## ------------------------------------------------------------------------
  ## Author: Marcel Wolbers, July 99
  ##
========================================================================
  require ("quadprog")
  m <- NCOL(x)
  if (length(eps)==1) eps <- rep(eps,m)
  x <- x * wsqrt
  y <- y * wsqrt
  #sometimes a rescaling of x and y helps (if solve.QP.compact fails
otherwise)
  xscale <- apply(abs(x),2,mean)
  yscale <- mean(abs(y))
  x <- t(t(x)/xscale)
  y <- y/yscale
  Rinv <- backsolve(qr.R(qr(x)),diag(m))
  cf <- solve.QP.compact(Dmat=Rinv,dvec=t(x)%*%y,Amat=rbind(rep(1,m)),
                   Aind=rbind(rep(1,m),1:m),bvec=eps*xscale/yscale,
                         factorized=TRUE)$sol
  cf <- cf*yscale/xscale  #scale back
  cf
}

##- #PS: Why the scaling in nnls.fit?
##- # Consider the following
##- y <-  c(1659799.4,2741065.6,652875.7)
##- x <-  matrix(c(24950.79,34825.50,14387.05,102962.14,15525.89,24074.87),
##- ncol=2)
##- nnls.fit(x,y) #would unfortunately give an error without the scaling











On Fri, 16 Nov 2001, Robert Abugov wrote:

>
> In July of 1999 Douglas Bates invited R users to implement an algorithm
for
> non-negative least squares based on Bates and Wolf, 1984: Communications
in
> Statistics, Part B 13:841-850.
> <http://www.ens.gu.edu.au/robertk/R/help/99b/0058.html>
>
> I'm wondering if anybody has implemented this or something similar so I
> won't have to reinvent the wheel.
>
> Thanks!
>
> Bob Abugov

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._