How to set limits using constrOptim in R?

Question

I am using constrOptim to minimize a log likelihood function for maximum likelihood estimation of parameters.

I wish to set the bounds on my parameters, but to not understand the constrOptim definition of the feasibility region.

The feasible region is defined by ui %*% theta - ci >= 0

I have a set of parameters with bounds [lower, upper]

a[0,5]  ie 0<a<5
b[0,Inf]
c[0,Inf]
e[0,1]

theta (starting values) = c(1, 1, 0.01,0.1)

What are the ui (constraint matrix (k x p)) and ci (constraint vector of length k) for these parameter bounds?

Is there a straightforward way to get from a list of upper and lower bounds to a ui and ci value?

Answer 1

Here's an example that we can use to illustrate ui and ci, with some extraneous output removed for brevity. It's maximizing the log likelihood of a normal distribution. In the first part, we use the optim function with box constraints, and in the second part, we use the constrOptim function with its version of the same box constraints.

# function to be optimized
> foo.unconstr <- function(par, x) -sum(dnorm(x, par[1], par[2], log=TRUE))

> x <- rnorm(100,1,1)

> optim(c(1,1), foo.unconstr, lower=c(0,0), upper=c(5,5), method="L-BFGS-B", x=x)
$par
[1] 1.147652 1.077654

$value
[1] 149.3724

> 
> # constrOptim example
> 
> ui <- cbind(c(1,-1,0,0),c(0,0,1,-1))
> ui
     [,1] [,2]
[1,]    1    0
[2,]   -1    0
[3,]    0    1
[4,]    0   -1
> ci <- c(0, -5, 0, -5)
> 
> constrOptim(c(1,1), foo.unconstr, grad=NULL, ui=u1, ci=c1, x=x)
$par
[1] 1.147690 1.077712

$value
[1] 149.3724

... blah blah blah ...

outer.iterations
[1] 2

$barrier.value
[1] -0.001079475

> 

If you look at the ui matrix and imagine multiplying by the parameter vector to be optimized, call it $\theta$, you'll see that the result has four rows, the first of which is $\theta_1$, the second $-\theta_1$, the third $\theta_2$, and the fourth $-\theta_2$. Subtracting off the ci vector and enforcing the $\ge 0$ constraint on each row results in $\theta_1 \ge 0$, $-\theta_1 + 5 \ge 0$, $\theta_2 \ge 0$ and $-\theta_2 + 5 \ge 0$. Obviously, multiplying the second and fourth constraints by -1 and moving the constant to the right hand side gets you to $\theta_1 \le 5$ and $\theta_2 \le 5$, the upper bound constraints.

Just substitute your own values into the ci vector and add appropriate columns (if any) to the ui vector to get the box constraint set you want.

Answer 2

Your constraints are of two types, either $\theta_i \geq a_i$, or $\theta_i \leq b_i$. The first ones are already in the right form (and the matrix ui is just the identity matrix), while the others can be written as $-\theta_i \geq - b_i$: ui is then $-I_n$ and ci is $-b$.

# Constraints
bounds <- matrix(c(
  0,5,
  0,Inf,
  0,Inf,
  0,1
), nc=2, byrow=TRUE)
colnames(bounds) <- c("lower", "upper")

# Convert the constraints to the ui and ci matrices
n <- nrow(bounds)
ui <- rbind( diag(n), -diag(n) )
ci <- c( bounds[,1], - bounds[,2] )

# Remove the infinite values
i <- as.vector(is.finite(bounds))
ui <- ui[i,]
ci <- ci[i]

# Constrained minimization
f <- function(u) sum((u+1)^2)
constrOptim(c(1,1,.01,.1), f, grad=NULL, ui=ui, ci=ci)

We can check how the constraint matrices ci and ui are interpreted:

# Print the constraints
k <- length(ci)
n <- dim(ui)[2]
for(i in seq_len(k)) {
  j <- which( ui[i,] != 0 )
  cat(paste( ui[i,j], " * ", "x[", (1:n)[j], "]", sep="", collapse=" + " ))
  cat(" >= " )
  cat( ci[i], "\n" )
}
# 1 * x[1] >= 0 
# 1 * x[2] >= 0 
# 1 * x[3] >= 0 
# 1 * x[4] >= 0 
# -1 * x[1] >= -5 
# -1 * x[4] >= -1 

Some of the algorithms in optim allow you to specify the lower and upper bounds directly: that is probably easier to use.