[R] Non linear optimization with nloptr package fail to produce true optimal result

[Duncan Murdoch]

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 using 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 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 https://www.R-project.org/posting-
>>>> guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
> 
>