# R/plugin_nlminb.R In ROI: R Optimization Infrastructure

#### Documented in nlminb2

```## ROI plugin: nlminb

################################################################################
## Utility Functions
################################################################################

## get_lb, get_ub taken from Florians nloptr plugins
## get_lb
## ======
##
## get lower bound constraints
get_lb <- function(x) {
##if( !length(bounds(x)\$lower\$val) ) {
##    lb <- NULL
##} else {
lb <- numeric( length(x\$objective) )
lb[ bounds(x)\$lower\$ind ] <- bounds(x)\$lower\$val
##}
return(lb)
}

## get_ub
## ======
##
## get upper bound constraints
get_ub <- function(x) {
##if( !length(bounds(x)\$upper\$val) ) {
##    ub <- NULL
##} else {
ub <- rep.int(Inf, length(x\$objective))
ub[ bounds(x)\$upper\$ind ] <- bounds(x)\$upper\$val
##}
return(ub)
}

################################################################################
## SOLVER METHODS
################################################################################

.solve_NLP_nlminb <- function( x, control ) {

## since this is contributed with ROI we have to specify the solver name directly
solver <- "nlminb"

## solve the NLP
## adjust arguments depending on problem class

## from nlminb2NLP() by Diethelm Wuertz

## Set Box Constraints:
lb <- get_lb(x)
ub <- get_ub(x)

env <- parent.frame()

FC <- constraints(x)\$F
dir <- constraints(x)\$dir
rhs <- constraints(x)\$rhs

idx_eq <- dir == "=="
if( length(dir) > sum(idx_eq) )
warning( "only equality constraints supported with nlminb, ignoring others." )

## Set Linear and Function Equality Constraints:
if ( any(idx_eq) ) {
eqfun <- function(x){
ans <- double(sum(idx_eq))
for (i in 1:sum(idx_eq))
ans[i] <- FC[idx_eq][[i]](x) - rhs[i]
return(as.double(eval(ans, env))) }
} else {
eqfun <- NULL
}

## TODO: Inequality Constraints (we assume that linear and quadratic terms have been coesced to F_constraints)
## if (!is.null(ineqA) || length(ineqFun) > 0) {
## leqfun <- function(x) {
##     ans <- NULL
##     if( length(ineqFun) > 0 )
##         for( i in 1:length(ineqFun) )
##             ans <- c(ans, ineqFun[[i]](x) - ineqFun.upper[i])
##     if (length(ineqFun) > 0)
##         for (i in 1:length(ineqFun))
##             ans <- c(ans, -ineqFun[[i]](x) + ineqFun.lower[i])
##     return( as.double(eval(ans, env)) )
##}
## } else {
leqfun <- NULL
## }

start <- control\$start
hessian <- control\$hessian
control[[ "start" ]] <- NULL
control[[ "hessian" ]] <- NULL

## now run nlminb2 solver
out <- nlminb2( start = start,
objective = objective(x),
eqFun = eqfun,
leqFun = leqfun,
upper = ub,
lower = lb,
hessian = hessian,
control = control )
ROI_plugin_canonicalize_solution( solution  = out\$par,
optimum   = objective(x)(out\$par),
status    = out\$convergence,
solver    = solver,
message   = out )
## ROI_plugin_canonicalize_solution( solution = out\$solution,
##                                    optimum = objective(x)(out\$solution),
##                                    status = out\$convergence,
##                                    solver = solver,
##                                    message = out )
}

################################################################################
## NOTE: this is work by Diethelm Wuertz taken out of the Rnlminb2
## package available at R-Forge:
## https://r-forge.r-project.org/scm/viewvc.php/pkg/Rnlminb2/
################################################################################

################################################################################
# FUNCTION:                DESCRIPTION:
#  nlminb2                  Nonlinear programming with nonlinear constraints
#  .Log                     Returns log taking care of negative values
################################################################################

##' Nonlinear programming with nonlinear constraints.
##'
##' This function was contributed by Diethelm Wuertz.
##' @param start numeric vector of start values.
##' @param objective the function to be minimized \eqn{f(x)}.
##' @param eqFun functions specifying equal constraints of the form
##' \eqn{h_i(x) = 0}. Default: \code{NULL} (no equal constraints).
##' @param leqFun functions specifying less equal constraints of the
##' form \eqn{g_i(x) <= 0}. Default: \code{NULL} (no less equal
##' constraints).
##' @param lower a numeric representing lower variable
##' bounds. Repeated as needed. Default: \code{-Inf}.
##' @param upper a numeric representing upper variable
##' bounds. Repeated as needed. Default: \code{Inf}.
##' @param hessian hessian of \eqn{f(x)}. Default: \code{NULL} (no
##' hessian provided).
##' @param control a list of control parameters. See
##' \code{\link[stats]{nlminb}()} for details. The parameter
##' \code{"scale"} is set here in contrast to
##' @author Diethelm Wuertz
##' @examples
##' ## Equal constraint function
##' eval_g0_eq <- function( x, params = c(1,1,-1)) {
##'        return( params[1]*x^2 + params[2]*x + params[3] )
##'    }
##' eval_f0 <- function( x, ... ) {
##'        return( 1 )
##'    }
##'
##'
##' @return list()
##' @importFrom "stats" "nlminb"
nlminb2 <- function( start, objective, eqFun = NULL, leqFun = NULL,
lower = -Inf, upper = Inf, gradient = NULL,
hessian = NULL, control = list() ) {

## Details:
##                        min f(x)
##    s.t.
##                 lower_i < x_i < upper_i
##                       h_i(x)  = 0
##                       g_i(x) <= 0

# Arguments:
#   start - numeric vector of start values
#   objective - objective function to be minimized f(x)
#   eqFun - equal constraint functions h_i(x) = 0
#   leqFun - less equal constraint functions g_i(x) <= 0
#   lower, upper - lower and upper bounds
#   hessian - optional hessian of f(x)
#   scale - control parameter
#   control - control list
#       eval.max - maximum number of evaluations (200)
#       iter.max - maximum number of iterations (150)
#       trace - value of the objective function and the parameters
#           is printed every trace'th iteration (0)
#       abs.tol - absolute tolerance (1e-20)
#       rel.tol - relative tolerance (1e-10)
#       x.tol - X tolerance (1.5e-8)
#       step.min - minimum step size (2.2e-14)

## TODO:  R, N and alpha should become part of the control list.

## FIXME: is it possible to set default start values (like 1/n)
stopifnot( is.numeric(start) )

## Control list:
ctrl <- .make_nlminb2_control_defaults()
if( length(control) > 0 )
for( name in names(control) )
ctrl[[ name ]] <- control[[ name ]]
control <- ctrl

# Minimization:
steps.tol <- control\$steps.tol
R <- control\$R
beta <- control\$beta
scale <- control\$scale

# Trace:
TRACE <- control\$trace > 0

# Reset Control:
control2 <- control
control2[["R"]] <- NULL
control2[["beta"]] <- NULL
control2[["steps.max"]] <- NULL
control2[["steps.tol"]] <- NULL
control2[["scale"]] <- NULL

## Unconstrained problem (inserted by st):
is_unconstrained <- is.null(eqFun) && is.null(leqFun)

if( is_unconstrained ){
ans <- nlminb( start = start, objective = objective,
scale = scale, control = control2,
lower = lower, upper = upper )
} else {
##
## Composed Objective Function:
## FIXME: we could use default ROI option gradient?
warning( "gradient not applicable for *constrained* NLPs for solver 'nlminb'." )
}
if( is.null(eqFun) ){
## type: "leq"
fun <- function( x, r ){
objective( x ) - r * sum( .Log(-leqFun(x)) ) }
} else if( is.null(leqFun) ){
## type: "eq"
fun <- function( x, r ){
objective( x ) + sum( (eqFun(x))^2 / r ) }
} else {
## type: "both"
fun <- function( x, r ){
objective( x ) + sum( (eqFun(x))^2 / r ) - r * sum( .Log(-leqFun(x)) ) }
}

counts <- 0
test <- 0
while (counts < control\$steps.max && test == 0) {
counts <- counts + 1
ans <- nlminb(
start = start, objective = fun,
scale = scale, control = control2, lower = lower, upper = upper,
r = R )
start <- ans\$par
tol <- abs((fun(ans\$par, R) - objective(ans\$par))/objective(ans\$par))
if (!is.na(tol))
if (tol < steps.tol) test <- 1
if (TRACE) {
cat("Iterations:  ", counts)
cat("\nConvergence: ", ans\$convergence)
cat("\nMessage:     ", ans\$message)
cat("\nSolution:    ", ans\$par)
cat("\nR-Objective: ", fun(start, R))
cat("\nFunction-Obj:", objective(start))
cat("\n")
}
R <- beta * R
}
}
# Return Value (plain return of nlminb()):
ans
}

##' Returns default control list for the nlminb solver.
##'
##' A function contributed by Diethelm Wuertz.
##' @return a list of control parameters.
##' @noRd
.make_nlminb2_control_defaults <- function()
list( eval.max = 500,
iter.max = 400,
trace = 0,
abs.tol =  1.0e-20,
rel.tol = 1.0e-10,
x.tol = 1.5e-8,
step.min = 2.2e-14,
scale = 1,
R = 1,
beta = 0.01,
steps.max = 10,
steps.tol = 1e-6 )

##' Returns log taking care of negative values.
##' @param x a numeric from which to calculate the log.
##' @return a numeric.
##' @noRd
.Log <-
function(x)
{
# Check for negative values:
x[x < 0] <- 0

# Return Value:
log(x)
}
# ------------------------------------------------------------------------------

################################################################################
## FIXME: we need for each problem class a separate solver method
## the following QP function is currently defunct.
.solve_QP_nlminb <- function( x, control ) {
##if( length(bounds(x)) )
##  constraints(x) <- rbind(constraints(x),
##                         .make_box_constraints_from_bounds(bounds(x),
##                                     dim(terms(objective(x))\$Q)[1]) )

solver <- "nlminb"
## solve the QP
## adjust arguments depending on problem class
out <- .nlminb_solve_QP( Q = terms(objective(x))\$Q,
L = terms(objective(x))\$L,
mat = constraints(x)\$L,
dir = constraints(x)\$dir,
rhs = constraints(x)\$rhs,
bounds = bounds(x),
max = x\$maximum,
#hessian = control\$hessian,
control = control )
ROI_plugin_canonicalize_solution( solution = out\$solution,
optimum = objective(x)(out\$solution),
status = out\$convergence,
solver = solver,
message = out )
}

.nlminb_solve_QP <- function(Q, L, mat, dir, rhs, bounds, max, gradient, control = list()) {

# nlminb does not directly support constraints
# we need to translate Ax ~ b constraints to lower, upper bounds
## FIXME: what about variable bounds????
stopifnot( is.null(bounds) )

A <- solve(t(mat))
n_obj <- ifelse( !is.null(Q),
dim(Q)[1],
length(L) )

## start
start <- as.numeric( control\$start )
if( !length(start) )
start <- slam::tcrossprod_simple_triplet_matrix( mat, matrix(rep(1/n_obj, n_obj), nrow = 1))
stopifnot( length(start) == n_obj )

lower <- rhs
upper <- c(Inf, Inf, Inf)

## possibly transformed objective function
foo <- function(x, L, A, Q) {
X = A %*% x
Objective = - slam::tcrossprod_simple_triplet_matrix(L, t(X)) + 0.5 * ( t(X) %*% slam::tcrossprod_simple_triplet_matrix(Q, t(X)))
Objective[[1]]
}
## FIXME: SPARSE!!! control list handling ok? what about "scale" parameter?
L = L, A = A, Q = Q,
control = control, lower = lower, upper = upper)
out\$solution <- as.numeric(A %*% out\$par)

# Return Value:
out
}
################################################################################

################################################################################
## STATUS CODES
################################################################################

## add all status codes generated by the solver to db

## Two examples are listed here:
0L,
"CONVERGENCE",
"Solution is optimal",
0L
)
1L,
"NON_CONVERGENCE",
"No solution."
)
invisible(TRUE)
}

################################################################################
## SIGNATURES
################################################################################

## based on nloptr interface
ROI_make_NLP_FXCV_signatures <- function()
ROI_plugin_make_signature( objective   = c("F"), # FIXME: c("L", "Q", "F"),
constraints = c("X", "F"), # FIXME: c("X", "L", "Q", "F"),
types       = c("C"),
bounds      = c("X", "V"),
cones       = c("X"),
maximum     = c(TRUE, FALSE) )

################################################################################
## SOLVER CONTROLS
################################################################################
## ROI + nlminb
ROI_plugin_register_solver_control( "nlminb", "trace", "verbose" )
ROI_plugin_register_solver_control( "nlminb", "iter.max", "max_iter" )
ROI_plugin_register_solver_control( "nlminb", "rel.tol", "tol" )
ROI_plugin_register_solver_control( "nlminb", "start", "start" )
## nlminb only
ROI_plugin_register_solver_control( "nlminb", "hessian", "X" )
ROI_plugin_register_solver_control( "nlminb", "eval.max", "X" )
ROI_plugin_register_solver_control( "nlminb", "abs.tol", "X" )
ROI_plugin_register_solver_control( "nlminb", "x.tol", "X" )
ROI_plugin_register_solver_control( "nlminb", "xf.tol", "X" )
ROI_plugin_register_solver_control( "nlminb", "step.min", "X" )
ROI_plugin_register_solver_control( "nlminb", "sing.tol", "X" )
ROI_plugin_register_solver_control( "nlminb", "scale.init", "X" )
ROI_plugin_register_solver_control( "nlminb", "diff.g", "X" )

invisible( TRUE )
}
```

## Try the ROI package in your browser

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

ROI documentation built on Jan. 27, 2018, 1:04 a.m.