Abby Spurdle
2021-May-22 07:42 UTC
[R] Testing optimization solvers with equality constraints
Sorry, this might sound like a poor question: But by "on the unit sphere", do you mean on the ***surface*** of the sphere? In which case, can't the surface of a sphere be projected onto a pair of circles? Where the cost function is reformulated as a function of two (rather than three) variables. On Sat, May 22, 2021 at 3:01 AM Hans W <hwborchers at 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 at 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.
Yes. "*on* the unit sphere" means on the surface, as you can guess from the equality constraint. And 'auglag()' does find the minimum, so no need for a special approach. I was/am interested in why all these other good solvers get stuck, i.e., do not move away from the starting point. And how to avoid this in general, not only for this specific example. On Sat, 22 May 2021 at 09:44, Abby Spurdle <spurdle.a at gmail.com> wrote:> > Sorry, this might sound like a poor question: > But by "on the unit sphere", do you mean on the ***surface*** of the sphere? > In which case, can't the surface of a sphere be projected onto a pair > of circles? > Where the cost function is reformulated as a function of two (rather > than three) variables. >