From m@ech|er @end|ng |rom @t@t@m@th@ethz@ch  Wed Feb  6 10:58:25 2019
From: m@ech|er @end|ng |rom @t@t@m@th@ethz@ch (Martin Maechler)
Date: Wed, 6 Feb 2019 10:58:25 +0100
Subject: [Rd] nlminb with constraints failing on some platforms
In-Reply-To: <7389c349-d3e9-d421-6907-f130149f138b@seanet.com>
References: <23636.2761.752118.1610@stat.math.ethz.ch>
 <7389c349-d3e9-d421-6907-f130149f138b@seanet.com>
Message-ID: <23642.44993.488084.34323@stat.math.ethz.ch>

Thank you, Brad (and others),

>>>>> Brad Bell on Mon, 4 Feb 2019 07:21:18 -0700 writes:

    > I get the failure message. To be specific:

    adcomp.git> R CMD BATCH --quiet test_nlminb.R
    adcomp.git> cat test_nlminb.Rout
    >> 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
    > [1] FALSE

ok... [I gave a version of the above which reveals a bit more ...]

    > My system is described by:
    adcomp.git> uname -a
    > Linux localhost.localdomain 4.17.7-200.fc28.x86_64 #1 SMP Tue Jul 17 16:28:31 UTC 2018 x86_64 x86_64 
    > x86_64 GNU/Linux

That (uname -a) is almost only a description of your kernel + hardware
(including of where the kernel was possibly built), and by the
way, almost equivalent to R's

   Sys.info()

In your case, I guess Fedora 28 Linux (which I also use, ..),
but we'd like a bit more information, notably  sessionInfo()

    > My version of R is described by:
    > Source?????? : R-3.5.2-2.fc28.src.rpm

The above also is not so informative.
I assume it means this is /usr/bin/R you got as a regular Fedora
package, and the above would actually by one of the output lines
of
   dnf inst /usr/bin/R

Indeed, now when I also try Fedora 28's /usr/bin/R,
I do see the same problem....
... and it is *also* *not* using R's own BLAS + LAPACK, but the
OpenBLAS  BLAS+LAPACK combination that comes with Fedora's R, 
 /usr/lib64/R/lib/libRblas.so :

> R.home()
[1] "/usr/lib64/R"
> 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] 0.797 0.303 0.285 0.271 0.258 ...
 $ objective  : num -37.7
 $ convergence: int 1
 $ iterations : int 150
 $ evaluations: Named int [1:2] 155 1611
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ message    : chr "iteration limit reached without convergence (10)"
> xhat <- rep(1, 10)                       # 64b Lnx/Win: -- 32bit even better
> all.equal(opt$par,         xhat,  tol=0) # good: 5.53 e-7
[1] "Mean relative difference: 2.233284"
> all.equal(opt$objective, f(xhat), tol=0) # good: 1.8 e-12
[1] "Mean relative difference: 0.0979214"
> print(abs(opt$objective- f(xhat) )) < 1e-4 # see 7.53 e-11  [Must be TRUE]
[1] 3.696533
[1] FALSE
> sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Fedora 28 (Twenty Eight)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=de_CH.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=de_CH.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=de_CH.UTF-8   
 [7] LC_PAPER=de_CH.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_CH.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] fortunes_1.5-4 sfsmisc_1.1-3 

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

And indeed BLAS/LAPACK  is  /usr/lib64/R/lib/libRblas.so
and that is described as

$ dnf info /usr/lib64/R/lib/libRblas.so
Updating .... repositories.
Last metadata expiration check: 0:08:10 ago on Wed 06 Feb 2019 09:31:40 AM CET.
Installed Packages
Name         : openblas-Rblas
Version      : 0.3.5
Release      : 1.fc28
Arch         : x86_64
Size         : 39 M
Source       : openblas-0.3.5-1.fc28.src.rpm
Repo         : @System
>From repo    : updates
Summary      : A version of OpenBLAS for R to use as libRblas
URL          : https://github.com/xianyi/OpenBLAS/
License      : BSD
Description  : 
             : OpenBLAS is an optimized BLAS library based on GotoBLAS2 1.13 BSD
             : version. The project is supported by the Lab of Parallel Software and
             : Computational Science, ISCAS. http://www.rdcps.ac.cn

---------------------------------------------------------------------------

I summarize what has been reported till:

Failure in these cases
========
1. Kasper K ("Scientific Linux", self compiled R, using Intel's MKL
  	    for BLAS/LAPACK)
2. (By Bill Dunlap): Microsoft R Open (MRO) 3.4.2, also using
  	    MKL with 12 cores
3. (By Brad Bell)  : R 3.5.2 Fedora 28 (x86_64) pkg, OpenBLAS(?)
4. (by MM)         : R 3.5.2 Fedora 28 (x86_64) pkg, BLAS+Lapack = OpenBLAS

Success
=======

- (by MM)	   : R-devel, R 3.5.2 patched on FC28, *self compiled* gcc 8.2,
      		     using R's BLAS/Lapack 
- (by Ralf Stubner): R 3.5.2 from Debian Stable (gcc 6.2) + OpenBLAS
- (by Berend H.)   : R 3.5.2 [from CRAN] on macOS 10.14.3 (BLAS/Lapack ??)
- (by MM)          : R-devel & R 3.5.2 (from CRAN) on Windows 32-bit & 64-b
- (by Rui Barradas): R 3.5.2 on ubuntu 18.04.1, BLAS/Lapack separate, /usr/lib/..
- (by Oliver Dechant)R 3.5.2 on Debian 9(stretch) + Intel MKL
- (by Avraham Adler):R 3.5.2 patched on Windows 10, Rblas=OpenBLAS 0.20, Lapack="base"
- (by Stefan Evert): R 3.5.1 on MacOS 10.13.6, system-supplied VecLib BLAS.

.............

and unfortunately I can not draw any clear conclusion from the
above.  It can't be MKL or OpenBLAS fault, alone.  
On the other hand, we've seen no failure when R's own BLAS was used.
Maybe MKL / OpenBLAS give no problems when only using 1 or 2
threads (and where doing that in the "good" cases), or could it
be that the behavior is even somewhat random on the platforms
with optimzing BLAS, depending on the CPUs' loads etc??

Apropos BLAS/Lapack :
  R's C+Fortran code does use (very little, I think) BLAS, but
  no Lapack routines.

Still, I don't  have much further clues, currently I think.
The only "failure" case, where R was 'self compiled' has been by
Kasper where he even found out that he could solve the problem
by using different F77 SAFE_FLAGS, and indeed these *are* used
when compiling <Rsource>/src/library/stats/portsrc.f ...

It would be great if this could be solved...

Martin


	
    > I have tried passing in the gradient and turning on the trace and it gives nearly the exact same trace with and without the gradient.
    [.......................]



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