[R] Testing optimization solvers with equality constraints
Gabor Grothendieck
ggrothend|eck @end|ng |rom gm@||@com
Fri May 28 01:13:50 CEST 2021
In case it is of interest this problem can be solved with an
unconstrained optimizer,
here optim, like this:
proj <- function(x) x / sqrt(sum(x * x))
opt <- optim(c(0, 0, 1), function(x) f(proj(x)))
proj(opt$par)
## [1] 5.388907e-09 7.071068e-01 7.071068e-01
On Fri, May 21, 2021 at 11:01 AM Hans W <hwborchers using gmail.com> wrote:
>
> Just by chance I came across the following example of minimizing
> a simple function
>
> (x,y,z) --> 2 (x^2 - y z)
>
> on the unit sphere, the only constraint present.
> I tried it with two starting points, x1 = (1,0,0) and x2 = (0,0,1).
>
> #-- Problem definition in R
> f = function(x) 2 * (x[1]^2 - x[2]*x[3]) # (x,y,z) |-> 2(x^2 -yz)
> g = function(x) c(4*x[1], 2*x[3], 2*x[2]) # its gradient
>
> x0 = c(1, 0, 0); x1 = c(0, 0, 1) # starting points
> xmin = c(0, 1/sqrt(2), 1/sqrt(2)) # true minimum -1
>
> heq = function(x) 1-x[1]^2-x[2]^2-x[3]^2 # staying on the sphere
> conf = function(x) { # constraint function
> fun = x[1]^2 + x[2]^2 + x[3]^2 - 1
> return(list(ceq = fun, c = NULL))
> }
>
> I tried all the nonlinear optimization solvers in R packages that
> allow for equality constraints: 'auglag()' in alabama, 'solnl()' in
> NlcOptim, 'auglag()' in nloptr, 'solnp()' in Rsolnp, or even 'donlp2()'
> from the Rdonlp2 package (on R-Forge).
>
> None of them worked from both starting points:
>
> # alabama
> alabama::auglag(x0, fn = f, gr = g, heq = heq) # right (inaccurate)
> alabama::auglag(x1, fn = f, gr = g, heq = heq) # wrong
>
> # NlcOptim
> NlcOptim::solnl(x0, objfun = f, confun = conf) # wrong
> NlcOptim::solnl(x1, objfun = f, confun = conf) # right
>
> # nloptr
> nloptr::auglag(x0, fn = f, heq = heq) # wrong
> # nloptr::auglag(x1, fn = f, heq = heq) # not returning
>
> # Rsolnp
> Rsolnp::solnp(x0, fun = f, eqfun = heq) # wrong
> Rsolnp::solnp(x1, fun = f, eqfun = heq) # wrong
>
> # Rdonlp2
> Rdonlp2::donlp2(x0, fn = f, nlin = list(heq), # wrong
> nlin.lower = 0, nlin.upper = 0)
> Rdonlp2::donlp2(x1, fn = f, nlin = list(heq), # right
> nlin.lower = 0, nlin.upper = 0) # (fast and exact)
>
> The problem with starting point x0 appears to be that the gradient at
> that point, projected onto the unit sphere, is zero. Only alabama is
> able to handle this somehow.
>
> I do not know what problem most solvers have with starting point x1.
> The fact that Rdonlp2 is the fastest and most accurate is no surprise.
>
> If anyone with more experience with one or more of these packages can
> give a hint of what I made wrong, or how to change calling the solver
> to make it run correctly, please let me know.
>
> Thanks -- HW
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com
More information about the R-help
mailing list