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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
