The Stark-Parker algorithm for bounded-variable least squares

Description

An R interface to the Stark-Parker implementation of an algorithm for bounded-variable least squares that solves \min{\parallel A x - b \parallel_2} with the constraint l \le x \le u, where l,x,u \in R^n, b \in R^m and A is an m \times n matrix.

References

Stark PB, Parker RL (1995). Bounded-variable least-squares: an algorithm and applications, Computational Statistics, 10, 129-141.

See Also

bvls, the method "L-BFGS-B" for optim, solve.QP, nnls

The Stark-Parker implementation of bounded-variable least squares

Description

An R interface to the Stark-Parker implementation of bounded-variable least squares that solves the least squares problem \min{\parallel A x - b \parallel_2} with the constraint l \le x \le u, where l,x,u \in R^n, b \in R^m and A is an m \times n matrix.

Usage

bvls(A, b, bl, bu, key=0, istate=rep(0,ncol(A)+1))

Arguments

Value

bvls returns an object of class "bvls".

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

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

Source

This is an R interface to the Fortran77 code accompanying the article referenced below by Stark PB, Parker RL (1995), and distributed via the statlib on-line software repository at Carnegie Mellon University (URL http://lib.stat.cmu.edu/general/bvls). The code was modified slightly to allow compatibility with the gfortran compiler. The authors have agreed to distribution under GPL version 2 or newer.

References

Stark PB, Parker RL (1995). Bounded-variable least-squares: an algorithm and applications, Computational Statistics, 10, 129-141.

See Also

the method "L-BFGS-B" for optim, solve.QP, nnls

Examples

## 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
X <- matrix(nrow = 50, ncol = 3) 
wavenum <- seq(18000,28000, length=nrow(X))
location <- c(25000, 22000) 
delta <- c(1000,1000)
Xcolfunc <- function(wavenum, location, delta)
  exp( - log(2) * (2 * (wavenum - location)/delta)^2)
for(i in 1:2) X[,i] <- Xcolfunc(wavenum, location[i], delta[i])

X[1:40,3] <- Xcolfunc(wavenum, 23000, 1000)[11:nrow(X)]
X[41:nrow(X),3]<- - Xcolfunc(wavenum, 23000, 1000)[21:30]

## 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 BVLS criteria 
bvls_sol <- function(matdat, A) {
  X <- matrix(0, nrow = ncol(matdat), ncol = ncol(A) )
  bu <- c(Inf,Inf,.75)
  bl <- c(0,0,-.75)
  for(i in 1:ncol(matdat)) 
     X[i,] <- coef(bvls(A,matdat[,i], bl, bu))
  X
}
X_bvls <- bvls_sol(matdat,A) 

matplot(X,type="p",pch=20)
matplot(X_bvls,type="l",pch=20,add=TRUE)
legend(10, -.5,
c("bound <= zero", "bound <= zero", "bound <= -.75 <= .75"),
col = c(1,2,3), lty=c(1,2,3),
text.col = "blue")

## 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 = ncol(matdat), ncol = 3)
  bu <- c(1000,1000,.75)
  bl <- c(0,0,-.75)
  for(i in 1:ncol(matdat))  
    X[i,] <-  optim(startval, fn = fn1, b=matdat[,i], A=A,
                  upper = bu, lower = bl,
                  method="L-BFGS-B")$par
    X
}
X_bfgs <- bfgs_sol(matdat,A) 

## the RMS deviation under BVLS is less than under L-BFGS-B 
sqrt(sum((X - X_bvls)^2)) < sqrt(sum((X - X_bfgs)^2))

## and L-BFGS-B is much slower 
system.time(bvls_sol(matdat,A))
system.time(bfgs_sol(matdat,A))

## End(Not run)