nlsic: Non Linear Least Squares with Inequality Constraints

In nlsic: Non Linear Least Squares with Inequality Constraints

Non Linear Least Squares with Inequality Constraints

Description

Solve non linear least squares problem min_par ||r(par,...)$res|| with optional inequality constraints u%*%par >= co and optional equality constraints e%*%par = eco

Usage

nlsic(
  par,
  r,
  u = NULL,
  co = NULL,
  control = list(),
  e = NULL,
  eco = NULL,
  flsi = lsi,
  ...
)

Arguments

par	initial values for parameter vector. It can be in non feasible domain, i.e. for which u%*%par >= co is not always respected.
r	a function calculating residual vector by a call to r(par, cjac, ...) where par is a current parameter vector, and cjac is set to TRUE if Jacobian is required or FALSE if not. A call to r() must return a list having a field "res" containing the residual vector and optionally a field "jacobian" when cjac is set to TRUE.
u	constraint matrix in u%*%par >= co
co	constraint right hand side vector
control	a list of control parameters ('=' indicates default values): errx=1.e-7 error on l2 norm of the iteration step sqrt(pt*p). maxit=100 maximum of newton iterations btstart=1 staring value for backtracking coefficient btfrac=0.5 (0;1) by this value we diminish the step till tight up to the quadratic model of norm reduction in backtrack (bt) iterations btdesc=0.1 (0;1) how good we have to tight up to the quadratic model 0 - we are very relaxe, 1 - we are very tight (may need many bt iterations) btmaxit=15 maximum of backtrack iterations btkmin=1.e-7 a floor value for backtracking fractioning trace=0 no information is printed during iterations (1 - print is active) rcond=1.e10 maximal condition number in QR decomposition starting from which a matrix is considered as numerically rank deficient. The inverse of this number is also used as a measure of very small number. ci = list of options relative to confidence interval estimation p=0.95 p-value of confidence interval. If you need a plain sd value, set p-value to 0.6826895 report=T report (or not and hence calculate or not) confidence interval information (in the field hci, as 'half confidence interval' width) history=TRUE report or not the convergence history adaptbt=FALSE use or not adaptive backtracking mnorm=NULL a norm matrix for a case sln==TRUE (cf next parameter) sln=FALSE use or not (default) least norm of the solution (instead of least norm of the increment) maxstep=NULL a maximal norm for increment step (if not NULL), must be positive monotone=FALSE should or not the cost decrease be monotone. If TRUE, then at first non decrease of the cost, the iterations are stopped with a warning message.
e	linear equality constraint matrix in e%*%par = eco (dense)
eco	right hand side vector in equality constraints
flsi	function solving linear least squares problem. For a custom function, see interfaces in lsi or lsi_ln help pages.
...	parameters passed through to r()

Details

Solving method consist in sequential LSI problems globalized by backtracking technique. If e, eco are not NULL, reduce jacobian to basis of e's kernel before lsi() call.
NB. If the function r() returns a list having a field "jacobian" it is supposed to be equal to the jacobian dr/dpar. If not, numerical derivation numDeriv::jacobian() is automatically used for its calculation.
NB2. nlsic() does not call stop() on possible errors. Instead, 'error' field is set to 1 in the returned result. This is done to allow a user to examine the current state of data and possibly take another path then to merely stop the program. So, a user must allways check this field at return from nlsic().
NB3. User should test that field 'mes' is not NULL even when error is 0. It may contain a warning message.

Value

a list with following components (some components can be absent depending on 'control' parameter)

• 'par' estimated values of par

• 'lastp' the last LSI solution during non linear iterations

• 'hci' vector of half-width confidence intervals for par

• 'ci_p' p-value for which CI was calculated

• 'ci_fdeg' freedom degree used for CI calculation

• 'sd_res' standard deviation of residuals

• 'covpar' covariance matrix for par

• 'laststep' the last increment after possible back-tracking iterations

• 'normp' the norm of lastp

• 'res' the last residual vector

• 'prevres' residual vector from previous non linear iteration

• 'jacobian' the last used jacobian

• 'retres' last returned result of r() call

• 'it' non linear iteration number where solution was obtained

• 'btit' back-tracking iteration number done during the last non linear iteration

• 'history' list with convergence history information

• 'error' error code: 0 - normal end, 1 - some error occurred, see message in 'mes'

• 'mes' textual message explaining what problem was in case of error

See Also

lsi, lsi_ln, uplo2uco

Examples

# solve min_{a,b} ||exp(a*x+b)-meas||, a,b>=1
a_true=1; b_true=2; x=0:5
# simulation function
sim=function(par, x) exp(par[["a"]]*x+par[["b"]])
# residual function to be passed to nlsic()
resid=function(par, cjac, ...) {
  dots=list(...)
  s=sim(par, dots$x)
  result=list(res=s-dots$meas)
  if (cjac) {
     result$jacobian=cbind(a=s*dots$x, b=s)
  }
  result
}
# simulated noised measurements for true parameters
set.seed(7) # for reproducible results
meas=sim(c(a=a_true, b=b_true), x)+rnorm(x)
# starting values for par
par=c(a=0, b=0)
# prepare constraints
uco=uplo2uco(par, lower=c(a=1, b=1))
# main call: solve the problem
fit=nlsic(par, resid, uco$u, uco$co, control=list(trace=1), x=x, meas=meas)
if (fit$error == 1) {
   stop(fit$mes)
} else {
   print(fit$par) # a=1.001590, b=1.991194
   if (! is.null(fit$mes)) {
      warning(fit$mes)
   }
}

nlsic documentation built on July 10, 2023, 2:03 a.m.