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

John Fox j|ox @end|ng |rom mcm@@ter@c@
Sat Dec 14 01:48:01 CET 2024 

Hello Duncan,

On 2024-12-13 5:30 p.m., Duncan Murdoch wrote:
> Caution: External email.
> 
> 
> 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

That was one of my first thoughts, but when I tried substituting the 
least-squares objective function, I observed similar behaviour.

> discontinuous derivative which will cause problems for many optimizers.
> It may also have many local optima, though I don't know in this case.

The model being fit is pretty simple -- a constrained linear model -- 
but as I said the columns of X are fairly highly correlated (though not 
so much as to create problems for unconstrained least-squares):

cor(X[, -1])
           x
x 1.0000000 0.9710284 0.9254732
   0.9710284 1.0000000 0.8201191
   0.9254732 0.8201191 1.0000000

Columns 2 and 3 are x^2 and log(x).

Also, the various solutions aren't all that different, and the 
regression residuals are relatively very small:

 > fits <- cbind(X %*% res$solution,
+               X %*% c(0.222, 6.999, 6.17, -19.371),
+               X %*% c(-0.610594, 4.307408, 6.254267, -11.919881),
+               y)
 > colnames(fits) = c("nloptr", "fmincon", "auglag", "y")
 > cor(fits)
            nloptr   fmincon    auglag         y
nloptr  1.0000000 0.9999997 0.9999990 0.9999930
fmincon 0.9999997 1.0000000 0.9999988 0.9999923
auglag  0.9999990 0.9999988 1.0000000 0.9999970
y       0.9999930 0.9999923 0.9999970 1.0000000
 > cor(log(fits))
            nloptr   fmincon    auglag         y
nloptr  1.0000000 0.9999796 0.9989490 0.9957987
fmincon 0.9999796 1.0000000 0.9991778 0.9961371
auglag  0.9989490 0.9991778 1.0000000 0.9986251
y       0.9957987 0.9961371 0.9986251 1.0000000

And take a look at (graph not shown):

plot(x, log(y), col="lightgray")
lines(x, log(fits[, 1]), col="magenta")
lines(x, log(fits[, 2]), col="blue")
lines(x, log(fits[, 3]), col="orange")
legend("bottomright", inset=0.025, col=c("magenta", "blue", "orange"),
        lty=1, legend=c("nloptr", "fmincon", "auglag"))

There's something peculiar going on at the lower-left.

Best,
  John

> 
> 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.
>>>>
>>>>
>>
>>
>

More information about the R-help mailing list