[Rd] nlminb with constraints failing on some platforms

[Martin Maechler]

>>>>> Kasper Kristensen via R-devel 
>>>>>     on Mon, 28 Jan 2019 08:56:39 +0000 writes:

    > I've noticed unstable behavior of nlminb on some Linux
    > systems. The problem can be reproduced by compiling
    > R-3.5.2 using gcc-8.2 and running the following snippet:

    > f <- function(x) sum( log(diff(x)^2+.01) + (x[1]-1)^2 )
    > opt <- nlminb(rep(0, 10), f, lower=-1, upper=3)
    > xhat <- rep(1, 10)
    > abs( opt$objective - f(xhat) ) < 1e-4  ## Must be TRUE

    > The example works perfectly when removing the bounds. However, when bounds are added the snippet returns 'FALSE'.

    > An older R version (3.4.4), compiled using the same gcc-8.2, did not have the problem. Between the two versions R has changed the flags to compile Fortran sources:

    > < SAFE_FFLAGS = -O2 -fomit-frame-pointer -ffloat-store
    > ---
    >> SAFE_FFLAGS = -O2 -fomit-frame-pointer -msse2 -mfpmath=sse

    > Reverting to the old SAFE_FFLAGS 'solves' the problem.

    >> sessionInfo()
    > R version 3.5.2 (2018-12-20)
    > Platform: x86_64-pc-linux-gnu (64-bit)
    > Running under: Scientific Linux release 6.4 (Carbon)

    > Matrix products: default
    > BLAS/LAPACK: /zdata/groups/nfsopt/intel/2018update3/compilers_and_libraries_2018.3.222/linux/mkl/lib/intel64_lin/libmkl_gf_lp64.so

    > locale:
    > [1] C

    > attached base packages:
    > [1] stats     graphics  grDevices utils     datasets  methods   base

    > loaded via a namespace (and not attached):
    > [1] compiler_3.5.2

So you us Intel's MKL library for BLAS/LAPACK ..

I also use gcc 8.2 (on Fedora 28 Linux) and R's own BLAS/LAPACK
and don't see such problems:

The code

f <- function(x) sum( log(diff(x)^2+.01) + (x[1]-1)^2 )
opt <- nlminb(rep(0, 10), f, lower=-1, upper=3)
str(opt)
xhat <- rep(1, 10)
all.equal(opt$par,         xhat,  tol=0) # good: 5.53 e-7
all.equal(opt$objective, f(xhat), tol=0) # good: 1.8 e-12
abs( opt$objective - f(xhat) ) < 1e-4  ## Must be TRUE

gives

> f <- function(x) sum( log(diff(x)^2+.01) + (x[1]-1)^2 )
> opt <- nlminb(rep(0, 10), f, lower=-1, upper=3)
> str(opt)
List of 6
 $ par        : num [1:10] 1 1 1 1 1 ...
 $ objective  : num -41.4
 $ convergence: int 0
 $ iterations : int 66
 $ evaluations: Named int [1:2] 96 830
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ message    : chr "relative convergence (4)"
> xhat <- rep(1, 10)
> all.equal(opt$par,         xhat,  tol=0) # good: 5.53 e-7
[1] "Mean relative difference: 5.534757e-07"
> all.equal(opt$objective, f(xhat), tol=0) # good: 1.8 e-12
[1] "Mean relative difference: 1.816536e-12"
> abs( opt$objective - f(xhat) ) < 1e-4  ## Must be TRUE
[1] TRUE
>

for me. Maybe others can quickly run the above  7 lines and report ?

Maybe there's something else unusual with your Linux
distribution's libraries?

I'm not an expert on these compiler flags; have you seen what
the R-admin manual
    https://cran.r-project.org/doc/manuals/R-admin.html#Linux
says about them?

Best,
Martin