R/gHgenb.R In optextras: Tools to Support Optimization Possibly with Bounds and Masks

Documented in gHgenb

```gHgenb <- function(par, fn, gr = NULL, hess = NULL, bdmsk = NULL,
lower = NULL, upper = NULL, control = list(ktrace=0), ...) {
# Generate the gradient and Hessian for a given function at the parameters
#   par
#  WITH recognition of bounds and masks.
#################################################################
#  NOTE: This routine generates the gradient and Hessian ignoring
#    bounds and masks and THEN applies a projection. This can cause
#    difficulties if the small steps taken to compute derivatives
#    give inadmissible parameter sets to functions.
#    We also apply symmetry check BEFORE the bounds/masks projection.
#    JN 2011-6-25
#    ?? How should that be dealt with generally?? 110625
#    ?? specifically can do 1 parameter at a time
#################################################################
#
# Input:
#  par = a vector for function parameters
#  fn = a user function (assumed to be sufficeintly differentiable)
# gr = name of a function to compute the (analytic) gradient of the user
#   function
# hess = name of a function to compute the (analytic) hessian of the user
#   function
#         This will rarely be available, but is included for completeness.
#  bdmsk = index of bounds and masks (see ?? for explanation)
#  lower = vector of lower bounds on the parameters par
#  upper = vector of upper bounds on the parameters par
#  control=list of controls:
#      asymtol = Tolerance to decide if we should force symmetry
#                Default is 1.0e-7.
#      ktrace = integer flag (default 0) to cause allow output if >0
#      stoponerror = logical flag (default=FALSE) to stop if we cannot
#            compute gradient or Hessian approximations. Otherwise return
#            with flags gradOK and/or hessOK set FALSE.
#   functions
#
# Output:
# A list containing
#    g = a vector containing the gradient
#    H = a matrix containing the Hessian
#    hessOK = logical indicator that Hessian is OK. TRUE/FALSE
#
#  Author: John Nash
#  Date:  January 14, 2011; updated June 25, 2011
#
#################################################################
#  IMPORTANT: There is a serious question of how to deal with Hessian
#    and gradient at active bounds, since these have meaning on one
#    side of the constraint only. JN 2011-6-25
#
#################################################################
ctrl <- list(asymtol = 1e-07, ktrace = 0, dowarn=TRUE, stoponerror = FALSE)
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
if (ctrl\$ktrace>0) ctrl\$dowarn<-TRUE # force TRUE for trace
# For bounds constraints, we need to 'project' the gradient and Hessian
# For masks there is no possibility of movement in parameter.
bmset <- sort(unique(c(which(par <= lower), which(par >= upper), c(which(bdmsk ==
0)))))
nbm <- length(bmset)  # number of elements nulled by bounds
#  Note: we assume that we are ON, not over boundary, but use <= and >=.
# No tolerance is used. ?? Do we want to consider parameters 'close' to
#   bounds?
if (ctrl\$ktrace > 0)
if (is.null(gr)) {
gn <- try(grad(fn, par, ...), silent = TRUE)  # change 20100711
}
else {
gn <- try(gr(par, ...), silent = TRUE)  # Gradient at solution # change 20100711
}
if (inherits(gn,"try-error"))
if (ctrl\$ktrace > 0)
print(gn)
if (nbm > 0) {
gn[bmset] <- 0  # 'project' the gradient
if (ctrl\$ktrace > 0) {
print(gn)
}
}
if (ctrl\$ktrace > 0)
cat("Compute Hessian approximation\n")
hessOK <- FALSE
if (is.null(hess)) {
if (ctrl\$ktrace > 0)
cat("is.null(hess) is TRUE\n")  # ???
if (is.null(gr)) {
if (ctrl\$ktrace > 0)
cat("is.null(gr) is TRUE use numDeriv hessian()\n")  # ???
Hn <- try(hessian(fn, par, ...), silent = TRUE)  # change 20100711
if (inherits(Hn, "try-error")) {
if (ctrl\$stoponerror)
stop("Unable to compute Hessian using numDeriv::hessian")
# hessOK still FALSE
}
else hessOK <- TRUE  # Do not need to check for symmetry either.
}
else {
if (ctrl\$ktrace > 0)
cat("is.null(gr) is FALSE use numDeriv jacobian()\n")  # ???
tHn <- try(Hn <- jacobian(gr, par, ...), silent = TRUE)  # change 20100711
if (inherits(tHn,"try-error")) {
if (ctrl\$stoponerror)
stop("Unable to compute Hessian using numderiv::jacobian")
if (ctrl\$ktrace > 0)
cat("Unable to compute Hessian using numderiv::jacobian \n")
}
else hessOK <- TRUE
if (ctrl\$ktrace > 0) {
cat("Hessian from Jacobian:")
# print(Hn) # Printed below
}
if (!isSymmetric(Hn, tol=0.1*sqrt(.Machine\$double.eps)))
{
asym <- sum(abs(t(Hn) - Hn))/sum(abs(Hn))
asw <- paste("Hn from jacobian is reported non-symmetric with asymmetry ratio ",
asym, sep = "")
if (ctrl\$ktrace > 0)
cat(asw, "\n")
if (ctrl\$dowarn) warning(asw)
if (asym > ctrl\$asymtol) {
if (ctrl\$stoponerror)
stop("Hessian too asymmetric")
}
else hessOK <- TRUE
if (ctrl\$ktrace > 0)
cat("Force Hessian symmetric\n")
else if (ctrl\$dowarn) warning("Hessian forced symmetric", call. = FALSE)
Hn <- 0.5 * (t(Hn) + Hn)
}  # end if ! isSymmetric
}  # numerical hessian at 'solution'
}
else {
if (ctrl\$ktrace > 0)
cat("is.null(hess) is FALSE -- trying hess()\n")  # ???
Hn <- try(hess(par, ...), silent = TRUE)  # change 20110222
if (inherits(Hn, "try-error")) {
if (ctrl\$stoponerror)
stop("Hessian evaluation with function hess() failed")
if (ctrl\$ktrace > 0)
cat("Hessian evaluation with function hess() failed \n")
}
else hessOK <- TRUE
print(Hn)
if (!isSymmetric(Hn))
{
asym <- sum(abs(t(Hn) - Hn))/sum(abs(Hn))
asw <- paste("Hn from hess() is reported non-symmetric with asymmetry ratio ",
asym, sep = "")
if (ctrl\$ktrace > 0)
cat(asw, "\n")
else if (ctrl\$dowarn) warning(asw)
if (asym > ctrl\$asymtol) {
if (ctrl\$stoponerror)
stop("Hessian too asymmetric")
}
else hessOK <- TRUE
if (ctrl\$ktrace > 0)
cat("Force Hessian symmetric\n")
else if (ctrl\$dowarn) warning("Hessian forced symmetric", call. = FALSE)
Hn <- 0.5 * (t(Hn) + Hn)
}  # end if ! isSymmetric
}  # end hessian computation
if (ctrl\$ktrace > 0)
print(Hn)
# apply the bounds and masks
if (nbm > 0) {
Hn[bmset, ] <- 0
Hn[, bmset] <- 0  # and the Hessian
if (ctrl\$ktrace > 0) {
print(Hn)
}
}
ansout <- list(gn, Hn, gradOK, hessOK, nbm)
names(ansout) <- c("gn", "Hn", "gradOK", "hessOK", "nbm")
return(ansout)
}  ## end gHgenb
```

Try the optextras package in your browser

Any scripts or data that you put into this service are public.

optextras documentation built on Dec. 31, 2019, 3:01 a.m.