Description Usage Arguments Value Source References See Also Examples

An implementation of an algorithm for linear least squares problems
with non-negative and non-positive
constraints based on the Lawson-Hanson
NNLS algorithm. Solves *\min{\parallel A x - b \parallel_2}*
with the constraint *x_i ≥ 0*
if *con_i ≥ 0* and *x_i ≤ 0* otherwise, where
*x, con \in R^n, b \in R^m*, and *A* is an *m \times n* matrix.

1 | ```
nnnpls(A, b, con)
``` |

`A` |
numeric matrix with |

`b` |
numeric vector of length |

`con` |
numeric vector of length |

`nnnpls`

returns
an object of class `"nnnpls"`

.

The generic accessor functions `coefficients`

,
`fitted.values`

, `deviance`

and `residuals`

extract
various useful features of the value returned by `nnnpls`

.

An object of class `"nnnpls"`

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 |

`bound` |
vector of the indices of |

`nsetp` |
the number of elements of |

This is an R interface to Fortran77 code distributed with the book referenced below by Lawson CL, Hanson RJ (1995), obtained from Netlib (file ‘lawson-hanson/all’), with some trivial modifications to allow for the combination of constraints to non-negative and non-positive values, and to return the variable NSETP.

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.

nnls, 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 | ```
## 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
## 2 of the Gaussian shapes are non-negative and 1 is non-positive
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])
X[,2] <- -X[,2]
## 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 NNNPLS criteria
nnnpls_sol <- function(matdat, A) {
X <- matrix(0, nrow = 51, ncol = 3)
for(i in 1:ncol(matdat))
X[i,] <- coef(nnnpls(A,matdat[,i],con=c(1,-1,1)))
X
}
X_nnnpls <- nnnpls_sol(matdat,A)
## Not run:
## can solve the same problem with L-BFGS-B algorithm
## but need starting values for x and
## impose a very low/high bound where none is desired
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, -1000, 0), upper=c(1000,0,1000),
method="L-BFGS-B")$par
X
}
X_bfgs <- bfgs_sol(matdat,A)
## the RMS deviation under NNNPLS is less than under L-BFGS-B
sqrt(sum((X - X_nnnpls)^2)) < sqrt(sum((X - X_bfgs)^2))
## and L-BFGS-B is much slower
system.time(nnnpls_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(c(1,-1,1))
# 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 NNNPLS is about the same as under quadprog
# sqrt(sum((X - X_nnnpls)^2))
# sqrt(sum((X - X_quadprog)^2))
## and quadprog requires about the same amount of time
# system.time(nnnpls_sol(matdat,A))
# system.time(quadprog_sol(matdat,A))
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.