John Fox
2024-Dec-13 22:11 UTC
[R] Non linear optimization with nloptr package fail to produce true optimal result
Dear Daniel, On 2024-12-13 2:51 p.m., Daniel Lobo wrote:> Caution: External email. > > > Looks like the solution 1.576708 6.456606 6.195305 -19.007996 is > the best solution that nloptr can produce by increasing the iteration > numbers. > > The better set of solution is obtained using pracma package.Not if I read the output correctly. As I showed, the result from pracma:fmincon() produces a larger value of the objective function than the result I obtained from nloptr(). John Nash (who is an expert on optimization -- I'm not) obtained an even lower value of the objective function from alabama::auglag(). As others have pointed out, one can't really draw general conclusions from a particular example, and like others, I don't have the time or inclination to figure out why your problem appears to be ill-conditioned (though note that the columns of X, excluding the constant, are highly correlated). Best, John> > On Sat, 14 Dec 2024 at 01:14, John Fox <jfox at mcmaster.ca> wrote: >> >> Dear Daniel et al., >> >> Following on Duncan's remark and examining the message produced by >> nloptr(), I simply tried increasing the maximum number of function >> evaluations: >> ------ snip ------- >> >> > nloptr(rep(0, 4), f, eval_g_ineq = hin, eval_g_eq = Hx, opts >> + list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-8, >> + maxeval = 1e5) >> + ) >> >> Call: >> >> nloptr(x0 = rep(0, 4), eval_f = f, eval_g_ineq = hin, eval_g_eq = Hx, >> opts = list(algorithm = "NLOPT_LN_COBYLA", xtol_rel = 1e-08, >> maxeval = 1e+05)) >> >> >> Minimization using NLopt version 2.7.1 >> >> NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped >> because xtol_rel or xtol_abs (above) was reached. ) >> >> Number of Iterations....: 46317 >> Termination conditions: xtol_rel: 1e-08 maxeval: 1e+05 >> Number of inequality constraints: 1 >> Number of equality constraints: 1 >> Optimal value of objective function: 1287.71725107671 >> Optimal value of controls: 1.576708 6.456606 6.195305 -19.008 >> >> ---------- snip ---------- >> >> That produces a solution closer to, and better than, the one that you >> suggested (which you obtained how?): >> >> > f(c(0.222, 6.999, 6.17, -19.371)) >> [1] 1325.076 >> >> I hope this helps, >> John >> -- >> John Fox, Professor Emeritus >> McMaster University >> Hamilton, Ontario, Canada >> web: https://www.john-fox.ca/ >> -- >> On 2024-12-13 1:45 p.m., Duncan Murdoch wrote: >>> Caution: External email. >>> >>> >>> You posted a version of this question on StackOverflow, and were given >>> advice there that you ignored. >>> >>> nloptr() clearly indicates that it is quitting without reaching an >>> optimum, but you are hiding that message. Don't do that. >>> >>> Duncan Murdoch >>> >>> On 2024-12-13 12:52 p.m., Daniel Lobo wrote: >>>> library(nloptr) >>>> >>>> set.seed(1) >>>> A <- 1.34 >>>> B <- 0.5673 >>>> C <- 6.356 >>>> D <- -1.234 >>>> x <- seq(0.5, 20, length.out = 500) >>>> y <- A + B * x + C * x^2 + D * log(x) + runif(500, 0, 3) >>>> >>>> #Objective function >>>> >>>> X <- cbind(1, x, x^2, log(x)) >>>> f <- function(theta) { >>>> sum(abs(X %*% theta - y)) >>>> } >>>> >>>> #Constraint >>>> >>>> eps <- 1e-4 >>>> >>>> hin <- function(theta) { >>>> abs(sum(X %*% theta) - sum(y)) - 1e-3 + eps >>>> } >>>> >>>> Hx <- function(theta) { >>>> X[100, , drop = FALSE] %*% theta - (120 - eps) >>>> } >>>> >>>> #Optimization with nloptr >>>> >>>> Sol = nloptr(rep(0, 4), f, eval_g_ineq = hin, eval_g_eq = Hx, opts >>>> list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-8))$solution >>>> # -0.2186159 -0.5032066 6.4458823 -0.4125948 >>> >>> ______________________________________________ >>> 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 https://www.R-project.org/posting- >>> guide.html >>> and provide commented, minimal, self-contained, reproducible code. >> >>
Duncan Murdoch
2024-Dec-13 22:30 UTC
[R] Non linear optimization with nloptr package fail to produce true optimal result
On 2024-12-13 5:11 p.m., John Fox wrote:> Dear Daniel, > > On 2024-12-13 2:51 p.m., Daniel Lobo wrote: >> Caution: External email. >> >> >> Looks like the solution 1.576708 6.456606 6.195305 -19.007996 is >> the best solution that nloptr can produce by increasing the iteration >> numbers. >> >> The better set of solution is obtained using pracma package. > > Not if I read the output correctly. As I showed, the result from > pracma:fmincon() produces a larger value of the objective function than > the result I obtained from nloptr(). > > John Nash (who is an expert on optimization -- I'm not) obtained an even > lower value of the objective function from alabama::auglag(). > > As others have pointed out, one can't really draw general conclusions > from a particular example, and like others, I don't have the time or > inclination to figure out why your problem appears to be ill-conditioned > (though note that the columns of X, excluding the constant, are highly > correlated).I would guess that the main difficulty with the example is the "sum of absolute errors" objective function. That objective function has a discontinuous derivative which will cause problems for many optimizers. It may also have many local optima, though I don't know in this case. Duncan Murdoch> > Best, > John > >> >> On Sat, 14 Dec 2024 at 01:14, John Fox <jfox at mcmaster.ca> wrote: >>> >>> Dear Daniel et al., >>> >>> Following on Duncan's remark and examining the message produced by >>> nloptr(), I simply tried increasing the maximum number of function >>> evaluations: >>> ------ snip ------- >>> >>> > nloptr(rep(0, 4), f, eval_g_ineq = hin, eval_g_eq = Hx, opts >>> + list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-8, >>> + maxeval = 1e5) >>> + ) >>> >>> Call: >>> >>> nloptr(x0 = rep(0, 4), eval_f = f, eval_g_ineq = hin, eval_g_eq = Hx, >>> opts = list(algorithm = "NLOPT_LN_COBYLA", xtol_rel = 1e-08, >>> maxeval = 1e+05)) >>> >>> >>> Minimization using NLopt version 2.7.1 >>> >>> NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped >>> because xtol_rel or xtol_abs (above) was reached. ) >>> >>> Number of Iterations....: 46317 >>> Termination conditions: xtol_rel: 1e-08 maxeval: 1e+05 >>> Number of inequality constraints: 1 >>> Number of equality constraints: 1 >>> Optimal value of objective function: 1287.71725107671 >>> Optimal value of controls: 1.576708 6.456606 6.195305 -19.008 >>> >>> ---------- snip ---------- >>> >>> That produces a solution closer to, and better than, the one that you >>> suggested (which you obtained how?): >>> >>> > f(c(0.222, 6.999, 6.17, -19.371)) >>> [1] 1325.076 >>> >>> I hope this helps, >>> John >>> -- >>> John Fox, Professor Emeritus >>> McMaster University >>> Hamilton, Ontario, Canada >>> web: https://www.john-fox.ca/ >>> -- >>> On 2024-12-13 1:45 p.m., Duncan Murdoch wrote: >>>> Caution: External email. >>>> >>>> >>>> You posted a version of this question on StackOverflow, and were given >>>> advice there that you ignored. >>>> >>>> nloptr() clearly indicates that it is quitting without reaching an >>>> optimum, but you are hiding that message. Don't do that. >>>> >>>> Duncan Murdoch >>>> >>>> On 2024-12-13 12:52 p.m., Daniel Lobo wrote: >>>>> library(nloptr) >>>>> >>>>> set.seed(1) >>>>> A <- 1.34 >>>>> B <- 0.5673 >>>>> C <- 6.356 >>>>> D <- -1.234 >>>>> x <- seq(0.5, 20, length.out = 500) >>>>> y <- A + B * x + C * x^2 + D * log(x) + runif(500, 0, 3) >>>>> >>>>> #Objective function >>>>> >>>>> X <- cbind(1, x, x^2, log(x)) >>>>> f <- function(theta) { >>>>> sum(abs(X %*% theta - y)) >>>>> } >>>>> >>>>> #Constraint >>>>> >>>>> eps <- 1e-4 >>>>> >>>>> hin <- function(theta) { >>>>> abs(sum(X %*% theta) - sum(y)) - 1e-3 + eps >>>>> } >>>>> >>>>> Hx <- function(theta) { >>>>> X[100, , drop = FALSE] %*% theta - (120 - eps) >>>>> } >>>>> >>>>> #Optimization with nloptr >>>>> >>>>> Sol = nloptr(rep(0, 4), f, eval_g_ineq = hin, eval_g_eq = Hx, opts >>>>> list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-8))$solution >>>>> # -0.2186159 -0.5032066 6.4458823 -0.4125948 >>>> >>>> ______________________________________________ >>>> 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 https://www.R-project.org/posting- >>>> guide.html >>>> and provide commented, minimal, self-contained, reproducible code. >>> >>> > >