Skip to content

nlminb with constraints failing on some platforms

12 messages · Ralf Stubner, Berend Hasselman, Martin Maechler +5 more

#
> 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
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)"
[1] "Mean relative difference: 5.534757e-07"
[1] "Mean relative difference: 1.816536e-12"
[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
#
On 01.02.19 10:00, Martin Maechler wrote:
Similar but not equal results for me:
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)"
[1] "Mean relative difference: 3.266165e-07"
[1] "Mean relative difference: 6.722005e-13"
[1] TRUE


Setup:

* Debian Stable with gcc 6.2
* R 3.5.2 from CRAN repository
* OpenBLAS

Greetings
Ralf
#
Identical result on R 3.5.2 on macOS 10.14.3 with the CRAN version of R.


Berend Hasselman
#

        
>> 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
#
Hello,

R 3.5.2 on ubuntu 18.04. sessionInfo() at the end.
Works with me, same results, cannot reproduce the error.


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
#[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


Hope this helps,

Rui Barradas


sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.1 LTS

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

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

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

loaded via a namespace (and not attached):
  [1] Rcpp_1.0.0       rstudioapi_0.8   bindr_0.1.1      magrittr_1.5
  [5] tidyselect_0.2.5 munsell_0.5.0    colorspace_1.3-2 lattice_0.20-38
  [9] R6_2.3.0         rlang_0.3.0.1    stringr_1.3.1    plyr_1.8.4
[13] dplyr_0.7.8      tools_3.5.2      grid_3.5.2       yaml_2.2.0
[17] assertthat_0.2.0 tibble_1.4.2     crayon_1.3.4     bindrcpp_0.2.2
[21] purrr_0.2.5      reshape2_1.4.3   glue_1.3.0       stringi_1.2.4
[25] compiler_3.5.2   pillar_1.3.1     scales_1.0.0     lubridate_1.7.4
[29] pkgconfig_2.0.2  zoo_1.8-4




?s 09:00 de 01/02/2019, Martin Maechler escreveu:
#
No error on Windows 10, R.3.5.2 patched, Rblas compiled with OpenBLAS
0.20, Rlapack is base.
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)"
[1] "Mean relative difference: 3.266165e-07"
[1] "Mean relative difference: 6.722005e-13"
[1] TRUE
R version 3.5.2 Patched (2018-12-26 r75909)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

loaded via a namespace (and not attached):
[1] compiler_3.5.2 tools_3.5.2    yaml_2.2.0
On Fri, Feb 1, 2019 at 3:24 PM Rui Barradas <ruipbarradas at sapo.pt> wrote:
#
Hello,
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS/LAPACK:
/opt/intel/compilers_and_libraries_2019.1.144/linux/mkl/lib/intel64_lin/libmkl_rt.so

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

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

loaded via a namespace (and not attached):
[1] compiler_3.5.2
List of 6
 $ par        : num [1:10] 1 1 1 1 1 ...
 $ objective  : num -41.4
 $ convergence: int 0
 $ iterations : int 67
 $ evaluations: Named int [1:2] 98 850
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ message    : chr "relative convergence (4)"
[1] "Mean relative difference: 5.349781e-07"
[1] "Mean relative difference: 3.812222e-12"
[1] TRUE
On 2019-02-01 4:23 p.m., Rui Barradas wrote:
#
Also no error on MacOS 10.13.6, R 3.5.1 with system-supplied VecLib BLAS.
2 days later
#
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

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

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

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.
Here is the output of a very similar case with the gradient:

 > n??? <- 10
 > f??? <- function(x) {
+?????? result <- 0.0
+?????? for( i in 2 : n ) {
+?????????????? result <- result + log( (x[i] - x[i-1])^2 + 0.01 ) + (x[1] - 1.0)^2
+?????? }
+?????? result
+ }
 > g??? <- function(x) {
+?????? result <- rep(0.0, n)
+?????? for( i in 2 : n ) {
+?????????????? result[1]?? <- result[1] + 2.0 * (x[1] - 1.0)
+?????????????? log_arg???? <- ( x[i] - x[i-1] )^2 + 0.01
+?????????????? log_arg_i?? <- 2.0 * (x[i] - x[i-1])
+?????????????? result[i]?? <- result[i]?? + log_arg_i / log_arg
+?????????????? result[i-1] <- result[i-1] - log_arg_i / log_arg
+?????? }
+?????? result
+ }
 > xstart <- rep(0.0, n)
 > opt??? <- nlminb(
+?????? xstart,
+?????? objective=f ,
+?????? gradient=g,
+?????? lower=-3,
+?????? upper=3,
+?????? control=list(trace=1)
+ )
 ? 0:??? -32.446532:? 0.00000? 0.00000? 0.00000? 0.00000? 0.00000 0.00000? 0.00000? 0.00000? 
0.00000? 0.00000
... snip ...
150:??? -37.750217: 0.796764 0.303221 0.285377 0.271175 0.257584 0.248540 0.239230 0.234184 0.229395 
0.227872
 > opt$par
 ?[1] 0.7967636 0.3032208 0.2853775 0.2711747 0.2575838 0.2485403 0.2392304
 ?[8] 0.2341841 0.2293946 0.2278722
 > g(opt$par)
 ?[1]? 0.23427575 -0.43398577 -0.67415314 -0.11550223 -0.87486481 0.05194325
 ?[7] -0.83926642 -0.05100054 -0.65128392 -0.30441806
 >
1 day later
#
Thank you, Brad (and others),
> 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 :
[1] "/usr/lib64/R"
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)"
[1] "Mean relative difference: 2.233284"
[1] "Mean relative difference: 0.0979214"
[1] 3.696533
[1] FALSE
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
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.
    [.......................]
#
.....
R 3.5.2 from CRAN using R's BLAS/Lapack.

Berend

....
#
If it helps, the BLAS I used is compiled to use 6 threads.
On Wed, Feb 6, 2019 at 3:47 AM Berend Hasselman <bhh at xs4all.nl> wrote: