From m@ech|er @end|ng |rom @t@t@m@th@ethz@ch  Fri Feb  1 17:31:42 2019
From: m@ech|er @end|ng |rom @t@t@m@th@ethz@ch (Martin Maechler)
Date: Fri, 1 Feb 2019 17:31:42 +0100
Subject: [Rd] nlminb with constraints failing on some platforms
In-Reply-To: <81BDC75A-ED23-4E70-A4AC-D34ADDCF6293@xs4all.nl>
References: <b83bb3a0e9674e198eb53d3682171804@dtu.dk>
 <23636.2761.752118.1610@stat.math.ethz.ch>
 <81BDC75A-ED23-4E70-A4AC-D34ADDCF6293@xs4all.nl>
Message-ID: <23636.29806.811582.405360@stat.math.ethz.ch>

>>>>> Berend Hasselman 
>>>>>     on Fri, 1 Feb 2019 15:59:58 +0100 writes:

    >> On 1 Feb 2019, at 10:00, Martin Maechler <maechler at stat.math.ethz.ch> wrote:
    >> 
    >> ........
    >>>> 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 ?
    >> 

    > Identical result on R 3.5.2 on macOS 10.14.3 with the CRAN version of R.

    > Berend Hasselman

Thank you, Berend, and Ralf (with gcc 6.2 and OpenBLAS),
for your reports.

In the mean time I also had tried Windows 32bit and 64bit, today's
R-devel version from CRAN (on a Terminal server) and also got
identical results for 64 bit and astonishingly even a bit more
accuracy for the 32 bit version.

So far, Kasper's platform is the only one where there's been a
problem, but I think it would be good to see/hear more here,
notably for people with optimized non-default BLAS/LAPACK or
system libraries...

    >> 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
    >> 
    >> ______________________________________________
    >> R-devel at r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel



[Rd] nlminb with constraints failing on some platforms
Martin Maechler