astamm / nloptr

nloptr provides an R interface to NLopt, a free/open-source library for nonlinear optimization providing a common interface to a number of different optimization routines which can handle nonlinear constraints and lower and upper bounds for the controls.
https://astamm.github.io/nloptr/
Other
106 stars 34 forks source link

Output for two auglag() examples in package reference manual is not reproducible #151

Closed adcascone closed 6 months ago

adcascone commented 9 months ago

I used the example from the alabama::auglag help page that's included on p. 8 of the nloptr reference manual to test the function, and I can't reproduce the $par values that are associated with the following code:

## Example from the alabama::auglag help page
fn <- function(x) (x[1] + 3*x[2] + x[3])^2 + 4 * (x[1] - x[2])^2
heq <- function(x) x[1] + x[2] + x[3] - 1
hin <- function(x) c(6*x[2] + 4*x[3] - x[1]^3 - 3, x[1], x[2], x[3])

auglag(runif(3), fn, hin = hin, heq = heq, localsolver="lbfgs")

The output I generated is:

# $par: 4.310718e-08 5.334520e-08 9.999999e-01

Similarly, I used the Powell problem example from the Rsolnp::solnp help page that's included on p. 8 of the nloptr reference manual, and I can't reproduce the $par values that are associated with the following code::

x0 <- c(-2, 2, 2, -1, -1)
fn1 <- function(x) exp(x[1]*x[2]*x[3]*x[4]*x[5])
eqn1 <-function(x)
c(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]+x[4]*x[4]+x[5]*x[5], x[2]*x[3]-5*x[4]*x[5], x[1]*x[1]*x[1]+x[2]*x[2]*x[2])

auglag(x0, fn1, heq = eqn1, localsolver = "mma")

The output I generated is:

# $par: 9.503717e-10 -9.744403e-09 -9.612198e-09 -1.082243e-10 -1.085369e-10
aadler commented 9 months ago

Hi.

Regarding the first question, you did effectively recover the results. First, we don't know the seed used in that example, but there is a call to runif so you are almost certainly not using the same starting values. That being said, you are getting a result of effectively c(0, 0, 1), which is the answer. You appear to be using the default tolerances, which are sqrt(.Machine$double.eps) or roughly 1.49e-8. Try running auglag(runif(3), fn, hin = hin, heq = heq, localsolver="lbfgs", control = list(xtol_rel = 1e-15)) and you will see your results for the first two parameters are within xtol_rel of 0, as they should be.

Regarding your second question, it is the same issue. You are getting a vector of five zeros within your requested tolerance. The exact value of your solution that close to zero depends not only on the processor you are using but also the BLAS. However, you are getting the correct answer within machine tolerance. This 1991 article by David Goldberg is very important in understanding the limits of how machines interpret very small real numbers.

Hope that helps.

aadler commented 6 months ago

@astamm, I think this one can be closed.

astamm commented 6 months ago

Thanks for the discussion.