# The Lawson-Hanson NNLS implemention of non-negative least squares

### Description

An R interface to the Lawson-Hanson NNLS implementation of an algorithm for non-negative linear least squares that solves \min{\parallel A x - b \parallel_2} with the constraint x ≥ 0, where x \in R^n, b \in R^m and A is an m \times n matrix.

### Usage

 1 nnls(A, b) 

### Arguments

 A numeric matrix with m rows and n columns b numeric vector of length m

### Value

nnls returns an object of class "nnls".

The generic accessor functions coefficients, fitted.values, deviance and residuals extract various useful features of the value returned by nnls.

An object of class "nnls" is a list containing the following components:

 x the parameter estimates. deviance the residual sum-of-squares. residuals the residuals, that is response minus fitted values. fitted the fitted values. mode a character vector containing a message regarding why termination occured. passive vector of the indices of x that are not bound at zero. bound vector of the indices of x that are bound at zero. nsetp the number of elements of x that are not bound at zero.

### Source

This is an R interface to the Fortran77 code distributed with the book referenced below by Lawson CL, Hanson RJ (1995), obtained from Netlib (file ‘lawson-hanson/all’), with a trivial modification to return the variable NSETP.

### References

Lawson CL, Hanson RJ (1974). Solving Least Squares Problems. Prentice Hall, Englewood Cliffs, NJ.

Lawson CL, Hanson RJ (1995). Solving Least Squares Problems. Classics in Applied Mathematics. SIAM, Philadelphia.

nnnpls, the method "L-BFGS-B" for optim, solve.QP, bvls
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 ## simulate a matrix A ## with 3 columns, each containing an exponential decay t <- seq(0, 2, by = .04) k <- c(.5, .6, 1) A <- matrix(nrow = 51, ncol = 3) Acolfunc <- function(k, t) exp(-k*t) for(i in 1:3) A[,i] <- Acolfunc(k[i],t) ## simulate a matrix X ## with 3 columns, each containing a Gaussian shape ## the Gaussian shapes are non-negative X <- matrix(nrow = 51, ncol = 3) wavenum <- seq(18000,28000, by=200) location <- c(25000, 22000, 20000) delta <- c(3000,3000,3000) Xcolfunc <- function(wavenum, location, delta) exp( - log(2) * (2 * (wavenum - location)/delta)^2) for(i in 1:3) X[,i] <- Xcolfunc(wavenum, location[i], delta[i]) ## set seed for reproducibility set.seed(3300) ## simulated data is the product of A and X with some ## spherical Gaussian noise added matdat <- A %*% t(X) + .005 * rnorm(nrow(A) * nrow(X)) ## estimate the rows of X using NNLS criteria nnls_sol <- function(matdat, A) { X <- matrix(0, nrow = 51, ncol = 3) for(i in 1:ncol(matdat)) X[i,] <- coef(nnls(A,matdat[,i])) X } X_nnls <- nnls_sol(matdat,A) matplot(X_nnls,type="b",pch=20) abline(0,0, col=grey(.6)) ## Not run: ## can solve the same problem with L-BFGS-B algorithm ## but need starting values for x bfgs_sol <- function(matdat, A) { startval <- rep(0, ncol(A)) fn1 <- function(par1, b, A) sum( ( b - A %*% par1)^2) X <- matrix(0, nrow = 51, ncol = 3) for(i in 1:ncol(matdat)) X[i,] <- optim(startval, fn = fn1, b=matdat[,i], A=A, lower = rep(0,3), method="L-BFGS-B")$par X } X_bfgs <- bfgs_sol(matdat,A) ## the RMS deviation under NNLS is less than under L-BFGS-B sqrt(sum((X - X_nnls)^2)) < sqrt(sum((X - X_bfgs)^2)) ## and L-BFGS-B is much slower system.time(nnls_sol(matdat,A)) system.time(bfgs_sol(matdat,A)) ## can also solve the same problem by reformulating it as a ## quadratic program (this requires the quadprog package; if you ## have quadprog installed, uncomment lines below starting with ## only 1 "#" ) # library(quadprog) # quadprog_sol <- function(matdat, A) { # X <- matrix(0, nrow = 51, ncol = 3) # bvec <- rep(0, ncol(A)) # Dmat <- crossprod(A,A) # Amat <- diag(ncol(A)) # for(i in 1:ncol(matdat)) { # dvec <- crossprod(A,matdat[,i]) # X[i,] <- solve.QP(dvec = dvec, bvec = bvec, Dmat=Dmat, # Amat=Amat)$solution # } # X # } # X_quadprog <- quadprog_sol(matdat,A) ## the RMS deviation under NNLS is about the same as under quadprog # sqrt(sum((X - X_nnls)^2)) # sqrt(sum((X - X_quadprog)^2)) ## and quadprog requires about the same amount of time # system.time(nnls_sol(matdat,A)) # system.time(quadprog_sol(matdat,A)) ## End(Not run)