Monday, September 07, 2015

Nonlinear constraints with a modified constrOptim

I often use constrOptim to quickly solve nonlinear optimization problems. constrOptim works well as a general tool to tackle constrained problems like
\min_{Ax -b \geq0} f(x)
There are many other options and packages for specific problems but constrOptim is likely to be the first choice when little is known on the problem at hand or for exploratory/quick-and-dirty analysis.

The main problem I have with it is the nature of the constraints: they must be linear (written in matrix form as \(Ax-b≥ 0\)). Often I would like to compute (local) solutions of more general problems of the form
\min_{g(x) -b \geq0} f(x)
where g(x) takes vector values, one component for each scalar constraint, and is in general nonlinear, i.e., \(g(x)\) cannot be written as \(Ax\).

Hence, I tweaked the code of constrOptim replacing all the occurrences of the linear constraints with my "new" non linear constraints. (In order to do that, I had to discard the chance to use "BFGS" as a solving method, for now...)  The result is a function, called constrOptimNL, see the code below, which can handle non linear constraints if a vector function g(x) and b are provided. All the drawbacks of the old function are inherited, that's life, but now non linear problems can be solved as shown by the following simple but non entirely trivial example in two dimensions.
 f <- function(x,y) x^2*y+x-x*y  
 fb <- function(x) f(x[1],x[2]) #vectorial version of f  
 x <- seq(0,3,length=31)  
 y <- seq(0,3,length=31)  
 z <- outer(x,y,f)  
 #draws a graph  
Heatmap and contour lines of the objective function, the feasible domain is filled in white.

Consider the constraints \(x\geq0,y\geq 1, x+y\leq 3\), the feasible domain is a triangle with vertices (0,0), (0,3), (3,0): due to linearity, the problem can be solved with

 A <- matrix(c(1,0,-1,0,1,-1),3,2)  
 b <- c(0,1,-3)  
 res <- constrOptim(c(0.5,1.5),fb,NULL,A,b)  
[1] 0.2792393 2.7207607
[1] -0.2683538
function gradient
246 NA 

The same problems can be solved with constrOptimNL, which trivially can handle linear constraints as well:
 g <- function(x,y) c(x,y,-x-y)  
 gb <- function(x) g(x[1],x[2])  
 resNL <- constrOptimNL(c(0.5,1.5),fb,NULL,gb,b)  
[1] 0.2792393 2.7207607
[1] -0.2683538
function gradient
246 NA

Assume now that some non linear constraint is involved in the optimization problem: \(x\geq0,y\geq1,(x-1)^2+(y-1)^2\leq1\), the feasible domain is now the upper half-circle centered at (1,1) with unit radius. This problem has nonlinearities in the constraints and cannot be solved with the standard constrOptim.
 g <- function(x,y) c(x,y,-(x-1)^2-(y-1)^2)  
 b <- c(0,1,-1)  
 resNL <- constrOptimNL(c(0.5,1.5),fb,NULL,gb,b)  
[1] 0.2649931 1.6780596 

A graph shows that this is indeed the correct solution:
 alpha <- seq(0,pi,len=31)  
 image(x,y,z); contour(x,y,z,add=T)  
 polygon(cos(alpha)+1,sin(alpha)+1,col="white",density=20) #draws the half-circle  
Contour levels and a nonlinear feasible domain in white (the solution is plotted with a filled circle).

You can check the solution more formally substituting the nonlinear constraint in polar coordinates, \(x=\cos\alpha+1,y=\sin\alpha+1\):
 f2 <- function(alpha) f(cos(alpha)+1,sin(alpha)+1)  
 sol <- optimize(f2,c(0,3))  
 #"same" as resNL  

Objective function restricted to the half-circle

The code

The code for constrOptimNL is given below (in this version the parameter grad is kept for compatibility with constrOptim but will never be used: Nelder-Mead and SANN are derivatives free and invoking BFGS stops the function). Hope this helps.
 constrOptimNL <- function (theta, f, grad, g, ci, mu = 1e-04, control = list(),   
   method = if (is.null(grad)) "Nelder-Mead" else "BFGS", outer.iterations = 100,   
   outer.eps = 1e-05, ..., hessian = FALSE)   
   if(method=="BFGS") stop("method BFGS not available, use Nelder-Mead or SANN")  
   if(!is.null(control$fnscale) && control$fnscale < 0)   
     mu <- -mu  
   R <- function(theta, theta.old, ...) {  
     ui.theta <- g(theta,...)  
     gi <- ui.theta - ci  
     if (any(gi < 0))   
     gi.old <- g(theta.old,...) - ci  
     bar <- sum(gi.old * log(gi) - ui.theta)  
     if (!is.finite(bar))   
       bar <- -Inf  
     f(theta, ...) - mu * bar  
   dR <- function(theta, theta.old, ...) {  
     ui.theta <- g(theta,...)  
     gi <- drop(ui.theta - ci)  
     gi.old <- drop(g(theta.old,...) - ci)  
     dbar <- colSums(ui * gi.old/gi - ui)  
     grad(theta, ...) - mu * dbar  
   if (any(g(theta,...) - ci <= 0))   
     stop("initial value is not in the interior of the feasible region")  
   obj <- f(theta, ...)  
   r <- R(theta, theta, ...)  
   fun <- function(theta, ...) R(theta, theta.old, ...)  
   gradient <- if (method == "SANN") {  
     if (missing(grad))   
     else grad  
   else function(theta,...) dR(theta, theta.old, ...)  
   totCounts <- 0 <- sign(mu)  
   for (i in seq_len(outer.iterations)) {  
     obj.old <- obj  
     r.old <- r  
     theta.old <- theta  
     a <- optim(theta.old, fun, gradient, control = control,   
       method = method, hessian = hessian, ...)  
     r <- a$value  
     if (is.finite(r) && is.finite(r.old) && abs(r - r.old) <   
       (0.001 + abs(r)) * outer.eps)   
     theta <- a$par  
     totCounts <- totCounts + a$counts  
     obj <- f(theta, ...)  
     if ( * obj > * obj.old)   
   if (i == outer.iterations) {  
     a$convergence <- 7  
     a$message <- gettext("Barrier algorithm ran out of iterations and did not converge")  
   if (mu > 0 && obj > obj.old) {  
     a$convergence <- 11  
     a$message <- gettextf("Objective function increased at outer iteration %d",   
   if (mu < 0 && obj < obj.old) {  
     a$convergence <- 11  
     a$message <- gettextf("Objective function decreased at outer iteration %d",   
   a$outer.iterations <- i  
   a$counts <- totCounts  
   a$barrier.value <- a$value  
   a$value <- f(a$par, ...)  
   a$barrier.value <- a$barrier.value - a$value