src/library/stats/tests/nlm.R

#  File src/library/stats/tests/nlm.R
#  Part of the R package, https://www.R-project.org
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

## nlm()  testing  ---- partly same as in ../demo/nlm.R
## NB: Strict regression tests -- output not "looked at"
library(stats)

## "truly 64 bit platform"
##  {have seen "x86-64" (instead of "x86_64") on Windows 2008 server}
(b.64 <- grepl("^x86.64", Sys.info()[["machine"]]) && .Machine$sizeof.pointer == 8)
(Lb64 <- b.64 && Sys.info()[["sysname"]] == "Linux")

## Example 1: Rosenbrock banana valley function (2 D)
##
f.Rosenb <- function(x1, x2) 100*(x2 - x1*x1)^2 + (1-x1)^2
grRosenb <- function(x1, x2) c(-400*x1*(x2 - x1*x1) - 2*(1-x1),
                               200*(x2 - x1*x1))
hessRosenb <- function(x1, x2) {
  a11 <- 2 - 400*x2 + 1200*x1*x1
  a21 <- -400*x1
  matrix(c(a11, a21, a21, 200), 2, 2)
}

fg <- function(x) {   # analytic gradient only
  x1 <- x[1]; x2 <- x[2]
  structure(f.Rosenb(x1, x2), "gradient" = grRosenb(x1, x2))
}
##
fgh <- function(x) {  # analytic gradient and Hessian
  x1 <- x[1];  x2 <- x[2]
  structure(f.Rosenb(x1, x2),
            "gradient" =  grRosenb(x1, x2),
            "hessian" = hessRosenb(x1, x2))
}

nlm3 <- function(x0, ...) {
    stopifnot(length(x0) == 2, is.numeric(x0))
    list(nl.f  = nlm(function(x) f.Rosenb(x[1],x[2]), x0, ...),
         nl.fg = nlm(fg , x0, ...),
         nl.fgh= nlm(fgh, x0, ...))
}

utils::str(l3.0 <- nlm3(x0 = c(-1.2, 1)))

chkNlm <- function(nlL, estimate, tols, codes.wanted = 1:2)
{
    stopifnot(is.list(nlL), ## nlL = list(<nlm>, <nlm>, <nlm>,...)
              sapply(nlL, is.list), lengths(nlL) == 5, # nlm(.) like
              is.numeric(estimate),
              is.list(tols), names(tols) %in% c("min","est","grad"),
              sapply(tols, is.numeric), unlist(tols) > 0)
    p <- length(estimate)
    n <- length(nlL)
    tols <- lapply(tols, rep_len, length.out = n)
    stopifnot(
        vapply(nlL, `[[`, pi, "minimum") <= tols$min,
        ##----
        abs(vapply(nlL, `[[`, estimate, "estimate") - estimate) <=
        rep(tols$est, each=p),
        ##----
        abs(vapply(nlL, `[[`, c(0,0), "gradient")) <= rep(tols$grad, each=p),
        ##----
        vapply(nlL, `[[`, 0L, "code") %in% codes.wanted)
}

chkNlm(l3.0, estimate = c(1,1),
       ##                  nl.f  nl.fg nl.fgh
       tols = list(min = c(1e-11,1e-17,1e-16),
                   est = c(4e-5, 1e-9, 1e-8),
                   grad= c(1e-6, 9e-9, 7e-7)))



## nl.fgh, the one with the Hessian had failed in R <= 3.4.0
## ------- and still is less accurate here than the gradient-only version

## all converge here, too,  fgh now being best
utils::str(l3.10 <- nlm3(x0 = c(-10, 10), ndigit = 14, gradtol = 1e-8))

## These tolerances were plucked from thin air: reduced for 32-bit Linux
chkNlm(l3.10, estimate = c(1,1), # lower tolerances now, notably for fgh:
       ##                  nl.f   nl.fg  nl.fgh
       tols = list(min = c(1e-9, 1e-20, 1e-16),
                   est = c(2e-5,  1e-10, 1e-14),
                   grad= c(1e-3,  6e-9 , 1e-12)),
       codes.wanted = if(Lb64) 1:2 else 1:3)

## all 3 fail to converge here
utils::str(l3.1c <- nlm3(x0 = c(-100, 100), iterlim = 1000))
## i.e., all convergence codes > 1:
sapply(l3.1c, `[[`, "code")
## nl.f  nl.fg nl.fgh  (seen on 32-bit and 64-bit)
##    2      2      4