Rvmmin 2015 version examples and tests

Rvmmin is an all-R version of the Fletcher-Nash variable metric nonlinear parameter optimization code of Fletcher (1970) ?? ref needed as modified by Nash (1979) ??ref needed.

Algorithm implementation

Fletcher's variable metric method attempts to mimic Newton's iteration for function minimization approximately.

Newton's method starts with an original set of parameters x[0]. At a given iteraion, which could be the first, we want to solve

x[k+1] = x[k] - H-1 g

where H is the Hessian and g is the gradient at x[k].

Newton's method is unattractive in general function minimization situations because

While the base Newton algorithm is as given, generally we carry out some sort of line search along the search direction delta from the current iterate x[k]. Indeed, many otherwise highly educated workers try to implement it without paying attention to safeguarding the iterations and ensuring appropriate progress towards a minimum.

Termination nuances

Termination variation with control tolerances

We use Chebyquad n = 4 test with different controls eps and acctol and tabulate the results.

cyq.f <- function (x) {
  rv<-cyq.res(x)
  f<-sum(rv*rv)
}

cyq.res <- function (x) {
# Fletcher's chebyquad function m = n -- residuals 
   n<-length(x)
   res<-rep(0,n) # initialize
   for (i in 1:n) { #loop over resids
     rr<-0.0
     for (k in 1:n) {
  z7<-1.0
  z2<-2.0*x[k]-1.0
        z8<-z2
        j<-1
        while (j<i) {
            z6<-z7
            z7<-z8
            z8<-2*z2*z7-z6 # recurrence to compute Chebyshev polynomial
            j<-j+1
        } # end recurrence loop
        rr<-rr+z8
      } # end loop on k
      rr<-rr/n
      if (2*trunc(i/2) == i) { rr <- rr + 1.0/(i*i - 1) }
      res[i]<-rr
    } # end loop on i
    res
}

cyq.jac<- function (x) {
#  Chebyquad Jacobian matrix
   n<-length(x)
   cj<-matrix(0.0, n, n)
   for (i in 1:n) { # loop over rows
     for (k in 1:n) { # loop over columns (parameters)
       z5<-0.0
       cj[i,k]<-2.0
       z8<-2.0*x[k]-1.0 
       z2<-z8
       z7<-1.0
       j<- 1
       while (j<i) { # recurrence loop
         z4<-z5
         z5<-cj[i,k]
         cj[i,k]<-4.0*z8+2.0*z2*z5-z4
         z6<-z7
         z7<-z8
         z8<-2.0*z2*z7-z6
         j<- j+1
       } # end recurrence loop
       cj[i,k]<-cj[i,k]/n
     } # end loop on k
   } # end loop on i
   cj
}


cyq.g <- function (x) {
   cj<-cyq.jac(x)
   rv<-cyq.res(x)
   gg<- as.vector(2.0* rv %*% cj)
}

require(Rvmmin)
## Loading required package: Rvmmin
nn <- 4
xx0 <- 1:nn
xx0 <- xx0 / (nn+1.0) # Initial value suggested by Fletcher

# cat("aed\n")
# aed <- Rvmminu(xx0, cyq.f, cyq.g, control=list(trace=2, checkgrad=FALSE))
# print(aed)
#================================
veps <- c(1e-3, 1e-5, 1e-7, 1e-9, 1e-11)
vacc <- c(.1, .01, .001, .0001, .00001, .000001)
resdf <- data.frame(eps=NA, acctol=NA, nf=NA, ng=NA, fval=NA, gnorm=NA)
for (eps in veps) {
  for (acctol in vacc) {
    ans <- Rvmminu(xx0, cyq.f, cyq.g, 
          control=list(eps=eps, acctol=acctol, trace=0))
    gn <- as.numeric(crossprod(cyq.g(ans$par)))
    resdf <- rbind(resdf, 
              c(eps, acctol, ans$counts[1], ans$counts[2], ans$value, gn))
  }
}
resdf <- resdf[-1,]
# print(resdf)
xtabs(formula = fval ~ acctol + eps, data=resdf)
##        eps
## acctol         1e-11        1e-09        1e-07        1e-05        0.001
##   1e-06 3.964816e-29 3.964816e-29 3.964816e-29 7.049696e-24 7.486504e-15
##   1e-05 3.964816e-29 3.964816e-29 3.964816e-29 7.049696e-24 7.486504e-15
##   1e-04 3.964816e-29 3.964816e-29 3.964816e-29 7.049696e-24 7.486504e-15
##   0.001 3.964816e-29 3.964816e-29 3.964816e-29 7.049696e-24 7.486504e-15
##   0.01  3.964816e-29 3.964816e-29 3.964816e-29 7.049696e-24 7.486504e-15
##   0.1   3.964816e-29 3.964816e-29 3.964816e-29 7.049696e-24 7.486504e-15
xtabs(formula = gnorm ~ acctol + eps, data=resdf)
##        eps
## acctol         1e-11        1e-09        1e-07        1e-05        0.001
##   1e-06 7.809261e-30 7.809261e-30 7.809261e-30 3.645064e-22 1.089927e-13
##   1e-05 7.809261e-30 7.809261e-30 7.809261e-30 3.645064e-22 1.089927e-13
##   1e-04 7.809261e-30 7.809261e-30 7.809261e-30 3.645064e-22 1.089927e-13
##   0.001 7.809261e-30 7.809261e-30 7.809261e-30 3.645064e-22 1.089927e-13
##   0.01  7.809261e-30 7.809261e-30 7.809261e-30 3.645064e-22 1.089927e-13
##   0.1   7.809261e-30 7.809261e-30 7.809261e-30 3.645064e-22 1.089927e-13
xtabs(formula = nf ~ acctol + eps, data=resdf)
##        eps
## acctol  1e-11 1e-09 1e-07 1e-05 0.001
##   1e-06    22    22    22    17    12
##   1e-05    22    22    22    17    12
##   1e-04    22    22    22    17    12
##   0.001    22    22    22    17    12
##   0.01     22    22    22    17    12
##   0.1      22    22    22    17    12
xtabs(formula = ng ~ acctol + eps, data=resdf)
##        eps
## acctol  1e-11 1e-09 1e-07 1e-05 0.001
##   1e-06    15    15    15    12     9
##   1e-05    15    15    15    12     9
##   1e-04    15    15    15    12     9
##   0.001    15    15    15    12     9
##   0.01     15    15    15    12     9
##   0.1      15    15    15    12     9

Problems of function scale

One of the more difficult aspects of termination decisions is that we need to decide when we have a “nearly” zero gradient. However, this “zero gradient” is relative to the overall scale of the function and its parameters.

ssq.f<-function(x){
   nn<-length(x)
   yy <- 1:nn
   f<-sum((yy-x/10^yy)^2)
   f
}
ssq.g <- function(x){
   nn<-length(x)
   yy<-1:nn
   gg<- 2*(x/10^yy - yy)*(1/10^yy)
}

xy <- c(1, 1/10, 1/100, 1/1000)
# note: checked gradient using numDeriv
veps <- c(1e-3, 1e-5, 1e-7, 1e-9, 1e-11)
vacc <- c(.1, .01, .001, .0001, .00001, .000001)
resdf <- data.frame(eps=NA, acctol=NA, nf=NA, ng=NA, fval=NA, gnorm=NA)
for (eps in veps) {
  for (acctol in vacc) {
    ans <- Rvmminu(xy, ssq.f, ssq.g, 
          control=list(eps=eps, acctol=acctol, trace=0))
    gn <- as.numeric(crossprod(ssq.g(ans$par)))
    resdf <- rbind(resdf, 
              c(eps, acctol, ans$counts[1], ans$counts[2], ans$value, gn))
  }
}
resdf <- resdf[-1,]
# print(resdf)
xtabs(formula = fval ~ acctol + eps, data=resdf)
##        eps
## acctol         1e-11        1e-09        1e-07        1e-05        0.001
##   1e-06 0.000000e+00 0.000000e+00 1.475416e-29 5.767419e-19 8.977439e-11
##   1e-05 0.000000e+00 0.000000e+00 1.475416e-29 5.767419e-19 8.977439e-11
##   1e-04 0.000000e+00 0.000000e+00 1.475416e-29 5.767419e-19 8.977439e-11
##   0.001 0.000000e+00 0.000000e+00 1.475416e-29 5.767419e-19 8.977439e-11
##   0.01  0.000000e+00 0.000000e+00 1.475416e-29 5.767419e-19 8.977439e-11
##   0.1   0.000000e+00 0.000000e+00 1.475416e-29 5.767419e-19 8.977439e-11
xtabs(formula = gnorm ~ acctol + eps, data=resdf)
##        eps
## acctol         1e-11        1e-09        1e-07        1e-05        0.001
##   1e-06 0.000000e+00 0.000000e+00 7.783028e-33 3.430257e-23 3.473135e-14
##   1e-05 0.000000e+00 0.000000e+00 7.783028e-33 3.430257e-23 3.473135e-14
##   1e-04 0.000000e+00 0.000000e+00 7.783028e-33 3.430257e-23 3.473135e-14
##   0.001 0.000000e+00 0.000000e+00 7.783028e-33 3.430257e-23 3.473135e-14
##   0.01  0.000000e+00 0.000000e+00 7.783028e-33 3.430257e-23 3.473135e-14
##   0.1   0.000000e+00 0.000000e+00 7.783028e-33 3.430257e-23 3.473135e-14
xtabs(formula = nf ~ acctol + eps, data=resdf)
##        eps
## acctol  1e-11 1e-09 1e-07 1e-05 0.001
##   1e-06    56    56    55    53    51
##   1e-05    56    56    55    53    51
##   1e-04    56    56    55    53    51
##   0.001    56    56    55    53    51
##   0.01     56    56    55    53    51
##   0.1      56    56    55    53    51
xtabs(formula = ng ~ acctol + eps, data=resdf)
##        eps
## acctol  1e-11 1e-09 1e-07 1e-05 0.001
##   1e-06    56    56    55    53    51
##   1e-05    56    56    55    53    51
##   1e-04    56    56    55    53    51
##   0.001    56    56    55    53    51
##   0.01     56    56    55    53    51
##   0.1      56    56    55    53    51

Weeds problem with random starts

This notorious problem (see Nash, 1979 or Nash, 2014 for details under the Hobbs Weeds problem) is small but generally difficult due to bad scaling and a near-singular Hessian in the original parameterization.

However, the Fletcher variable metric method can solve this problem quite well. It is important to ensure there are enough iterations to allow the method to “grind” at the problem. If one uses default settings for various aspects of the termination criteria (maxit in optim:BFGS or stopbadupdate=TRUE in Rvmminu), then the success rate drops to less than 2/3 of cases tried below.

## hobbstarts.R -- starting points for Hobbs problem 
hobbs.f<- function(x){ # # Hobbs weeds problem -- function
    if (abs(12*x[3]) > 500) { # check computability
       fbad<-.Machine$double.xmax
       return(fbad)
    }
    res<-hobbs.res(x)
    f<-sum(res*res)
##    cat("fval =",f,"\n")
##    f
}
hobbs.res<-function(x){ # Hobbs weeds problem -- residual
# This variant uses looping
    if(length(x) != 3) stop("hobbs.res -- parameter vector n!=3")
    y<-c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443, 
         38.558, 50.156, 62.948, 75.995, 91.972)
    t<-1:12
    if(abs(12*x[3])>50) {
       res<-rep(Inf,12)
    } else {
       res<-x[1]/(1+x[2]*exp(-x[3]*t)) - y
    }
}

hobbs.jac<-function(x){ # Jacobian of Hobbs weeds problem
   jj<-matrix(0.0, 12, 3)
   t<-1:12
    yy<-exp(-x[3]*t)
    zz<-1.0/(1+x[2]*yy)
     jj[t,1] <- zz
     jj[t,2] <- -x[1]*zz*zz*yy
     jj[t,3] <- x[1]*zz*zz*yy*x[2]*t
   return(jj)
}

hobbs.g<-function(x){ # gradient of Hobbs weeds problem
    # NOT EFFICIENT TO CALL AGAIN
    jj<-hobbs.jac(x)
    res<-hobbs.res(x)
    gg<-as.vector(2.*t(jj) %*% res)
    return(gg)
}
require(Rvmmin)
set.seed(12345)
nrun<-100
sstart<-matrix(runif(3*nrun, 0, 5), nrow=nrun, ncol=3)
ustart<-sstart %*% diag(c(100, 10, 0.1))
nsuccR <- 0
nsuccO <- 0
vR <- rep(NA, nrun)
vO <- vR
fR <- vR
gR <- vR
fO <- vR
gO <- vR

for (irun in 1:nrun) {
  us <- ustart[irun,]
  print(us)
#  ans <- Rvmminu(us, hobbs.f, hobbs.g, control=list(trace=1))
#  ans <- optim(us, hobbs.f, hobbs.g, method="BFGS")
  ans <- Rvmminu(us, hobbs.f, hobbs.g, control=list(trace=0))
  ao <- optim(us, hobbs.f, hobbs.g, method="BFGS", 
               control=list(maxit=3000))
# ensure does not max function out

cat(irun,"  Rvmminu value =",ans$value,"  optim:BFGS value=",ao$value,"\n")
  if (ans$value < 2.5879) nsuccR <- nsuccR + 1
  if (ao$value < 2.5879) nsuccO <- nsuccO + 1
#  tmp <- readline()
  vR[irun] <- ans$value
  vO[irun] <- ao$value
  fR[irun] <- ans$counts[1]
  gR[irun] <- ans$counts[2]
  fO[irun] <- ao$counts[1]
  gO[irun] <- ao$counts[2]

}
## [1] 360.4519481  14.7232702   0.2942962
## 1   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 437.8865965  30.8626833   0.4462959
## 2   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 380.49116416  48.71370636   0.06189745
## 3   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 443.0622831  30.9106019   0.2566545
## 4   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 228.2404801  26.0684597   0.3318201
## 5   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 83.185893 45.125190  0.382771
## 6   Rvmminu value = 2.587277   optim:BFGS value= 2.58728 
## [1] 162.54769336  31.87272163   0.04631309
## 7   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 254.61216783  43.21505665   0.03382227
## 8   Rvmminu value = 2.587277   optim:BFGS value= 2.587279 
## [1] 363.8526269  12.5558870   0.2780492
## 9   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 494.8684690  10.7534540   0.3256592
## 10   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 17.2677175 30.4738022  0.1479813
## 11   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 76.1867451 19.1722318  0.3477524
## 12   Rvmminu value = 2.587277   optim:BFGS value= 2.587281 
## [1] 367.842476  37.763552   0.188835
## 13   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1]  0.5682933 18.9868121  0.1304810
## 14   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 195.6016676  39.7486032   0.4700795
## 15   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 231.247327  45.284557   0.215445
## 16   Rvmminu value = 2.587277   optim:BFGS value= 2.58729 
## [1] 194.0719908  49.2013087   0.1635323
## 17   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 201.2425709  29.3973987   0.2674541
## 18   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 89.4817924  0.4732026  0.3129708
## 19   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 475.82937684  16.03960346   0.08521099
## 20   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 226.86403664  27.97283666   0.05465815
## 21   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 163.37620432  25.74587801   0.05363939
## 22   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 482.70766169   4.48559311   0.04637089
## 23   Rvmminu value = 2.587277   optim:BFGS value= 2.587279 
## [1] 353.7409386  34.1625287   0.2672878
## 24   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 322.2713183  35.4986025   0.4481946
## 25   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 194.914242  40.011578   0.355013
## 26   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 349.2718197  47.2039207   0.4532271
## 27   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 272.0289322   9.3323069   0.2183443
## Warning in Rvmminu(us, hobbs.f, hobbs.g, control = list(trace = 0)): Too
## many gradient evaluations
## 28   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 113.2335893  13.9321141   0.2339094
## 29   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 242.2788775  38.7131644   0.2990598
## 30   Rvmminu value = 2.587277   optim:BFGS value= 2.58728 
## [1] 396.5035849   1.6723807   0.1533739
## 31   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1]  2.9938148 33.5959070  0.2860004
## 32   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 93.8562229 40.6383561  0.2396016
## 33   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 340.9168116  33.1580619   0.3483085
## 34   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 185.0520618  47.2237336   0.3676399
## 35   Rvmminu value = 2.587277   optim:BFGS value= 2.587279 
## [1] 180.8127851  43.7422034   0.2582818
## 36   Rvmminu value = 2.587277   optim:BFGS value= 2.58728 
## [1] 434.3974493  34.8685599   0.4926857
## 37   Rvmminu value = 2.587277   optim:BFGS value= 2.587279 
## [1] 452.0773329   9.2520446   0.4589896
## 38   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 308.712283   4.877216   0.236244
## 39   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 67.0158173 40.3879995  0.1903133
## 40   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 391.0966403  26.0888818   0.4694052
## 41   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 214.5994100  34.6770338   0.4722385
## 42   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 463.6369875  10.8164795   0.4534123
## 43   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 386.6216124  31.4971436   0.3921888
## 44   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 129.8406234   7.3503939   0.3131655
## 45   Rvmminu value = 2.587277   optim:BFGS value= 2.587279 
## [1] 160.6123367  47.9535230   0.3378219
## 46   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 30.0975787 49.5055223  0.4685549
## 47   Rvmminu value = 2.587277   optim:BFGS value= 2.587321 
## [1] 21.7282269 24.8219654  0.2672636
## 48   Rvmminu value = 2.587277   optim:BFGS value= 2.587279 
## [1] 27.526909 45.982382  0.168395
## 49   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 312.7713985  18.6065070   0.3953354
## 50   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 482.23514436  41.35287229   0.06163229
## 51   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 413.6514345  27.9521474   0.4795642
## 52   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 157.514118  39.783560   0.016246
## 53   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 106.5127255  17.6010873   0.1298446
## 54   Rvmminu value = 2.587277   optim:BFGS value= 2.587311 
## [1] 366.2480594  10.4263836   0.4382151
## 55   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 249.6205103  39.2390312   0.4228122
## 56   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 3.648860e+02 3.415256e+01 8.426952e-03
## 57   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 40.16802215 21.60055070  0.09182044
## 58   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 217.76524244  42.32169935   0.05792107
## 59   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 118.2902267  30.0971920   0.2392073
## 60   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 395.7838992  37.0330258   0.4748624
## 61   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 129.3421584  29.7451001   0.4128688
## 62   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 492.9919159  42.5779545   0.4057794
## 63   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 378.4368718  40.6950504   0.4578685
## 64   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 489.8891236  19.0211676   0.3500826
## 65   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 109.47391961  24.61541012   0.04151929
## 66   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 474.35359156  49.66887323   0.05903807
## 67   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 74.7289731  8.7044725  0.3968248
## 68   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 300.1784853  41.7140400   0.4980219
## 69   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 473.2153759  19.7907687   0.3455178
## 70   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 344.17667845  48.45280671   0.07366959
## 71   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 252.7668616   9.4103233   0.1589239
## 72   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 186.46812241  37.46108349   0.07798945
## 73   Rvmminu value = 2.587277   optim:BFGS value= 2.58728 
## [1] 167.9025192  15.4234805   0.4976268
## 74   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 24.1256763 18.9529451  0.3741136
## 75   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 309.4737695  37.5382922   0.1965197
## 76   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 480.7236459  35.2100838   0.2756957
## 77   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 327.48025481  46.02965168   0.01470503
## 78   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 255.145996  39.763866   0.168151
## 79   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 75.0491053 41.4538482  0.4291896
## 80   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 435.22394018   8.37738313   0.06354407
## 81   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 257.2208389  20.4979890   0.3957472
## 82   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 4.32395679 9.81263703 0.04894095
## 83   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1]  9.597382741 12.970031216  0.009576733
## 84   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 72.2559508 48.5208133  0.1380112
## 85   Rvmminu value = 2.587277   optim:BFGS value= 2.587411 
## [1] 152.515877  40.582160   0.205576
## 86   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 412.82843135  17.37979198   0.03335091
## 87   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 251.1723212  46.9172764   0.3982658
## 88   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 401.7863129   8.1693447   0.2531626
## 89   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 30.3199911  7.8404580  0.2593119
## 90   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 463.97757  37.70316   0.35273
## 91   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 404.0892755  44.3025201   0.4158095
## 92   Rvmminu value = 2.587277   optim:BFGS value= 18.45938 
## [1] 39.40666991 15.33820682  0.03044822
## 93   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 300.4641408  15.5540874   0.1684179
## 94   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 357.3899403  49.2262595   0.4011205
## 95   Rvmminu value = 2.587277   optim:BFGS value= 2.587281 
## [1] 256.9228266   0.3464615   0.1417513
## 96   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 360.0582525  13.7160289   0.4074386
## 97   Rvmminu value = 2.587277   optim:BFGS value= 2.587277 
## [1] 374.9731350  39.9022107   0.1075416
## 98   Rvmminu value = 2.587277   optim:BFGS value= 2.587278 
## [1] 47.8202619 12.1834859  0.2775764
## 99   Rvmminu value = 2.587277   optim:BFGS value= 2.587316 
## [1] 1.989129e+02 2.291367e+01 7.763851e-03
## 100   Rvmminu value = 2.587277   optim:BFGS value= 2.587279
cat("Rvmminu: number of successes=",nsuccR,"  propn=",nsuccR/nrun,"\n")
## Rvmminu: number of successes= 100   propn= 1
cat("optim:BFGS no. of successes=",nsuccO,"  propn=",nsuccO/nrun,"\n")
## optim:BFGS no. of successes= 99   propn= 0.99
fgc <- data.frame(fR, fO, gR, gO)
# print(fgc)
summary(fgc)
##        fR              fO               gR               gO       
##  Min.   : 41.0   Min.   :  58.0   Min.   : 26.00   Min.   : 16.0  
##  1st Qu.:105.8   1st Qu.: 140.5   1st Qu.: 39.00   1st Qu.: 53.0  
##  Median :155.5   Median : 184.0   Median : 53.00   Median : 68.5  
##  Mean   :205.7   Mean   : 323.5   Mean   : 59.57   Mean   :131.2  
##  3rd Qu.:258.0   3rd Qu.: 453.5   3rd Qu.: 66.00   3rd Qu.:178.8  
##  Max.   :920.0   Max.   :1427.0   Max.   :507.00   Max.   :610.0

bounds and masks

Let us make sure that Rvmminb is doing the right thing with bounds and masks.

Bounds

bt.f<-function(x){
 sum(x*x)
}

bt.g<-function(x){
  gg<-2.0*x
}

lower <- c(0, 1, 2, 3, 4)
upper <- c(2, 3, 4, 5, 6)
bdmsk <- rep(1,5)
xx <- rep(0,5) # out of bounds
ans <- Rvmmin(xx, bt.f, bt.g, lower=lower, upper=upper, bdmsk=bdmsk)
## Warning in Rvmmin(xx, bt.f, bt.g, lower = lower, upper = upper, bdmsk =
## bdmsk): Parameter out of bounds has been moved to nearest bound
ans
## $par
## [1] 0 1 2 3 4
## 
## $value
## [1] 30
## 
## $counts
## function gradient 
##        1        1 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "Rvmminb appears to have converged"
## 
## $bdmsk
## [1]  1 -3 -3 -3 -3

Masks

sq.f<-function(x){
   nn<-length(x)
   yy<-1:nn
   f<-sum((yy-x)^2)
   f
}
sq.g <- function(x){
   nn<-length(x)
   yy<-1:nn
   gg<- 2*(x - yy)
}

xx0 <- rep(pi,3)
bdmsk <- c(1, 0, 1) # Middle parameter fixed at pi
cat("Check final function value (pi-2)^2 = ", (pi-2)^2,"\n")
## Check final function value (pi-2)^2 =  1.303234
require(Rvmmin)
ans <- Rvmmin(xx0, sq.f, sq.g, lower=-Inf, upper=Inf, bdmsk=bdmsk,
               control=list(trace=2))
## Bounds: nolower =  TRUE   noupper =  TRUE  bounds =  TRUE 
## Gradient test with tolerance =  6.055454e-06 
## Analytic gradient uses function  gr 
## function at parameters =  5.909701  with attributes:
## NULL
## Compute analytic gradient
## [1] 4.2831853 2.2831853 0.2831853
## Compute numeric gradient
## [1] 4.2831853 2.2831853 0.2831853
## gradient test tolerance =  6.055454e-06   fval= 5.909701 
##  compare to max(abs(gn-ga))/(1+abs(fval)) =  3.242827e-12 
## admissible =  TRUE 
## maskadded =  FALSE 
## parchanged =  FALSE 
## Bounds: nolower =  FALSE   noupper =  FALSE  bounds =  TRUE 
## Rvmminb -- J C Nash 2009-2015 - an R implementation of Alg 21
## Problem of size n= 3   Dot arguments:
## list()
## Initial fn= 5.909701 
##   1   1   5.909701 
## Gradproj = -18.42587 
## reset steplength= 1 
## *reset steplength= 0.2 
## ig= 2   gnorm= 2.575522     3   2   2.961562 
## Gradproj = -15.04576 
## reset steplength= 1 
## *reset steplength= 0.2 
## ig= 3   gnorm= 0.23879     5   3   1.317489 
## Gradproj = -0.02851034 
## reset steplength= 1 
## ig= 4   gnorm= 0   Small gradient norm
## Seem to be done Rvmminb
ans
## $par
## [1] 1.000000 3.141593 3.000000
## 
## $value
## [1] 1.303234
## 
## $counts
## function gradient 
##        6        4 
## 
## $convergence
## [1] 2
## 
## $message
## [1] "Rvmminb appears to have converged"
## 
## $bdmsk
## [1] 1 0 1
ansnog <- Rvmmin(xx0, sq.f, lower=-Inf, upper=Inf, bdmsk=bdmsk,
               control=list(trace=2))
## Bounds: nolower =  TRUE   noupper =  TRUE  bounds =  TRUE 
## WARNING: forward gradient approximation being used
## admissible =  TRUE 
## maskadded =  FALSE 
## parchanged =  FALSE 
## Bounds: nolower =  FALSE   noupper =  FALSE  bounds =  TRUE 
## Rvmminb -- J C Nash 2009-2015 - an R implementation of Alg 21
## Problem of size n= 3   Dot arguments:
## list()
## WARNING: using gradient approximation ' grfwd '
## Initial fn= 5.909701 
##   1   1   5.909701 
## Gradproj = -18.42587 
## reset steplength= 1 
## *reset steplength= 0.2 
## ig= 2   gnorm= 2.575522     3   2   2.961562 
## Gradproj = -15.04576 
## reset steplength= 1 
## *reset steplength= 0.2 
## ig= 3   gnorm= 0.23879     5   3   1.317489 
## Gradproj = -0.02851034 
## reset steplength= 1 
## ig= 4   gnorm= 2.668644e-08     6   4   1.303234 
## Gradproj = -4.446061e-16 
## reset steplength= 1 
## *reset steplength= 0.2 
## *reset steplength= 0.04 
## *reset steplength= 0.008 
## *reset steplength= 0.0016 
## *reset steplength= 0.00032 
## *reset steplength= 6.4e-05 
## *reset steplength= 1.28e-05 
## *reset steplength= 2.56e-06 
## *reset steplength= 5.12e-07 
## *reset steplength= 1.024e-07 
## Unchanged in step redn 
## No acceptable point
## Reset to gradient search
##   16   4   1.303234 
## Gradproj = -7.121661e-16 
## reset steplength= 1 
## *reset steplength= 0.2 
## *reset steplength= 0.04 
## *reset steplength= 0.008 
## *reset steplength= 0.0016 
## *reset steplength= 0.00032 
## *reset steplength= 6.4e-05 
## *reset steplength= 1.28e-05 
## *reset steplength= 2.56e-06 
## *reset steplength= 5.12e-07 
## *reset steplength= 1.024e-07 
## Unchanged in step redn 
## No acceptable point
## Converged 
## Seem to be done Rvmminb
ansnog
## $par
## [1] 1.000000 3.141593 3.000000
## 
## $value
## [1] 1.303234
## 
## $counts
## function gradient 
##       26        4 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "Rvmminb appears to have converged"
## 
## $bdmsk
## [1] 1 0 1

special tests