[R] Testing optimization solvers with equality constraints

[Gabor Grothendieck]

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