bvls {bvls} | R Documentation |

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))
```

`A` |
numeric matrix with |

`b` |
numeric vector of length |

`bl` |
numeric vector of length |

`bu` |
numeric vector of length |

`key` |
If |

`istate` |
vector of length |

`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:

`x` |
the parameter estimates. |

`deviance` |
the residual sum-of-squares. |

`residuals` |
the residuals, that is response minus fitted values. |

`fitted` |
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 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.

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

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
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)
```

[Package *bvls* version 1.4 Index]