[R] Constrained Optimization in R (alabama)

[Ravi Varadhan]

Berend,

I had pointed out this solution to Axel already, but he aparently wants to solve more general problems like this, so he wanted a solution from constrained optimimzers. But, if the more general problems have the similarly discrete set of candidate parameter values, using constrained optimizers like alabama would still NOT be good approaches. It would be easiest to just go through the discrete values.

Ravi

________________________________________
From: Berend Hasselman [bhh at xs4all.nl]
Sent: Monday, February 11, 2013 12:45 AM
To: Axel Urbiz
Cc: R-help at r-project.org; Ravi Varadhan
Subject: Re: [R] Constrained Optimization in R (alabama)

On 10-02-2013, at 21:16, Axel Urbiz <axel.urbiz at gmail.com> wrote:

> Dear List,
>
> I'm trying to solve this simple optimization problem in R. The parameters
> are the exponents to the matrix mm. The constraints specify that each row
> of the parameter matrix should sum to 1 and their product to 0. I don't
> understand why the constraints are not satisfied at the solution. I must be
> misinterpreting how to specify the constrains somehow.
>
>
> library(alabama)
>
> ff <- function (x) {
>  mm <- matrix(c(10, 25, 5, 10), 2, 2)
>  matx <- matrix(x, 2, 2)
>  -sum(apply(mm ^ matx, 1, prod))
> }
>
>
> ### constraints
>
> heq <- function(x) {
>  h <- rep(NA, 1)
>  h[1] <- x[1] + x[3] -1
>  h[2] <- x[2] + x[4] -1
>  h[3] <- x[1] * x[3]
>  h[4] <- x[2] * x[4]
>  h
> }
>
> res <- constrOptim.nl(par = c(1, 1, 1, 1), fn = ff,
>                      heq = heq)
> res$convergence #why NULL?
> matrix(round(res$par, 2), 2, 2)  #why constraints are not satisfied?


There is no need to use an optimization function in this case. You can easily find the optimum by hand.

>From the constraints :

x[3] equals (1-x[1])
x[4] equals (1-x[2])

Therefore

x[1] * (1-x[1]) equals 0
x[2] * (1-x[2]) equals 0

So you only need to calculate the function value for  4 vectors:

x1 <- c(1,1,0,0)
x2 <- c(0,0,1,1)
x3 <- c(1,0,0,1)
x4 <- c(0,1,1,0)

Modifying your ff function in the way Ravi suggested:

ff <- function (x) {
 mm <- matrix(c(10, 25, 5, 10), 2, 2)
 matx <- matrix(x, 2, 2)
 sum(apply(mm ^ matx, 1, prod))
}

and running ff with x{1,2,3,4} as argument, one gets

> ff(x1)
[1] 35
> ff(x2)
[1] 15
> ff(x3)
[1] 20
> ff(x4)
[1] 30

So ff(x2) achieves the minimum, which corresponds to Ravi's result with auglag.

Berend