bvls {bvls}R Documentation

The Stark-Parker implementation of bounded-variable least squares


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.


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



numeric matrix with m rows and n columns


numeric vector of length m


numeric vector of length n specifying the lower bound on each element of x


numeric vector of length n specifying the upper bound on each element of x


If key > 0 the routine initializes using the user's guess about which components of x are active, i.e. are strictly within their bounds, which are at their lower bounds, and which are at their upper bounds. This information is supplied through the array istate.


vector of length ncol(A)+1. If key > 0, istate is as follows: the last contains the total number of components at their bounds (the bound variables). The absolute values of the first nbound <- tail(istate,1) entries of istate are the indices of these bound components of x. The sign of istate[1:nbound] indicates whether x(abs(istate[1:nbound])) is at its upper or lower bound. istate[1:nbound] is positive if the component is at its upper bound, negative if the component is at its lower bound. istate[(nbound+1):ncol(A)] contain the indices of the components of x that are active (i.e. are expected to lie strictly within their bounds). When key > 0, the routine initially sets the active components to the averages of their upper and lower bounds.


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:


the parameter estimates.


the residual sum-of-squares.


the residuals, that is response minus fitted values.


the fitted values.


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 The code was modified slightly to allow compatibility with the gfortran compiler. The authors have agreed to distribution under GPL version 2 or newer.


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


## 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

## 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_bvls <- bvls_sol(matdat,A) 

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,
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 

## End(Not run)

[Package bvls version 1.4 Index]