Bug 17703 - constrOptim() fails for a well-specified problem

[hturner]

Bug 17703 documents a case where constrOptim() throws an error for a well-spedified constrained optimization problem. The error is thrown because the objective function becomes non-finite due to imprecision caused by floating-point arithmetic. In particular, the error is triggered by a check for negative values in object gi on L32 (prior to a log(gi) computation that would produce NaN for negative values).

@mmaechler indicates that a patch would be welcome along the following lines:

• Avoid false precision by removing noise from gi before checking, e.g. with signif(gi, digits), where digits is a new argument
• Update the corresponding help file constrOptim.Rd

Note: bug report specifies this as an issue in R-3.6.2, I have confirmed this is still present in R-devel (2025/05/15, r88203).

[astamm]

The bug reported on Bugzilla seems to be due to the gradient function taking an extremely large value, ultimately making gradient-based optimisers provide values that might slightly exit the feasible region (up to a floating-point precision error).

A solution could be to issue a warning when this happens suggesting the user to either review his implementation of the gradient or use a derivative-free optimiser (e.g. Nelder-Mead).

We could also use Nelder-Mead (and warn the user about it) for these specific iterations where the optimiser selected by the user would provide an unfeasible solution.

[tdhock]

We made the following figure which shows the function to optimize.

The min is in the bottom right (circle).
The last point in the optimization before the error is nearby (triangle).
The gradients are very strong near the boundaries due to the logs in the objective function.
gradfun(c( 0.953770837623602,0.0462281623763985 ) = c( -1000001.00054123,-999979.36870582 )
With such strong gradients, I'm not sure we should expect gradient-based algorithms to work.
We observed that Nelder-Mead (no gradients necessary) works but BFGS does not, so I would propose fixing this issue by adding a warning that suggests to use Nelder-Mead instead.

Code to reproduce figure:

constrOptimMine <-
    function(theta, f, grad, ui, ci, mu = 0.0001, control = list(),
             method = if(is.null(grad)) "Nelder-Mead" else "BFGS",
             outer.iterations = 100, outer.eps = 0.00001, ..., hessian = FALSE,
             tol.constraint = 1e-6)
{

    if (!is.null(control$fnscale) && control$fnscale < 0)
        mu <- -mu ##maximizing

    R <- function(theta, theta.old, ...) {
        ui.theta <- ui%*%theta
        gi <-  ui.theta-ci
        ##if (any(gi < -tol.constraint)) {
        it.dt.list[[length(it.dt.list)+1]] <<- data.table(
          g1=gi[1],
          g2=gi[2])
        if (any(gi < 0)) {
          ## print(gi)
          ## browser()
          return(NaN)
        }
        ##gi[gi<=0] <- tol.constraint
        gi.old <- ui%*%theta.old-ci
        bar <- sum(gi.old*log(gi)-ui.theta)
        if (!is.finite(bar)) bar <-  -Inf
        f(theta, ...) -mu*bar
    }

    dR <- function(theta, theta.old, ...) {
        ui.theta <- ui%*%theta
        gi <- drop(ui.theta-ci)
        gi.old <- drop(ui%*%theta.old-ci)
        dbar <- colSums(ui*gi.old/gi-ui)
        grad(theta, ...) - mu*dbar
    }

    if (any(ui%*%theta-ci <= 0))
        stop("initial value is not in the interior of the feasible region")
    obj <- f(theta, ...)
    r <- R(theta, theta, ...)
    fun <- function(theta, ...) R(theta, theta.old, ...)
    gradient <- if(method == "SANN") {
        if(missing(grad)) NULL else grad
    } else function(theta, ...) dR(theta, theta.old, ...)
    totCounts <- 0
    s.mu <- sign(mu)

    for(i in seq_len(outer.iterations)) {
        obj.old <- obj
        r.old <- r
        theta.old <- theta

        a <- optim(theta.old, fun, gradient, control = control,
                   method = method, hessian = hessian, ...)
        r <- a$value
        if (is.finite(r) && is.finite(r.old) &&
 	    abs(r - r.old) < (1e-3 + abs(r)) * outer.eps) break
        theta <- a$par
 	totCounts <- totCounts + a$counts
        obj <- f(theta, ...)
        if (s.mu * obj > s.mu * obj.old) break
    }
    if (i == outer.iterations) {
        a$convergence <- 7
        a$message <- gettext("Barrier algorithm ran out of iterations and did not converge")
    }
    if (mu > 0 && obj > obj.old) {
        a$convergence <- 11
        a$message <- gettextf("Objective function increased at outer iteration %d", i)
    }
    if (mu < 0 && obj < obj.old) {
        a$convergence <- 11
        a$message <- gettextf("Objective function decreased at outer iteration %d", i)
    }
    a$outer.iterations <- i
    a$counts <- totCounts
    a$barrier.value <- a$value
    a$value <- f(a$par, ...)
    a$barrier.value <- a$barrier.value - a$value
    a
}

`objective` <- function(p){
    out <- -log(p[1] + p[2]) + log(p[2]) + log(1 - p[1] - p[2])
    ## following line for debugging [not needed to reproduce error]
    cat(paste("objective(c(",paste(p,collapse=","),") =",out,"\n")) 
    return(out)
}
`gradfun` <- function(p){
  out <- c(
  -1/(p[1] + p[2])          - 1/(1 - p[1] - p[2]),
  -1/(p[1] + p[2]) + 1/p[2] - 1/(1 - p[1] - p[2])
  )
    ## following line for debugging [not needed to reproduce error]
  cat(paste("gradfun(c(",paste(p,collapse=","),") = c(",paste(out,collapse=","),")\n"))
  return(out)
}
small <- 1e-6
small <- 1e-6
CI <- c(small,small,small-1)
UI <- rbind(
  c( 1, 0),  # p1>small
  c( 0, 1),  # p2>small
  c(-1,-1)   # p1+p2 < 1-small
)
library(data.table)
library(ggplot2)
grid.vec <- seq(0,1,l=101)
grid.dt <- data.table::CJ(p1=grid.vec, p2=grid.vec)[, let(
  value = -log(p1 + p2) + log(p2) + log(1 - p1 - p2),
  constraints = ifelse(
    p1>small & p2>small & p1+p2 < 1-small,
    "satisfied",
    "violated")
)][]
(best.theta <- c(1-small, small))
best.theta <- c( 0.564012436588663,0.435987563411338 )
objective(best.theta)
gradfun(best.theta)
UI %*% best.theta - CI
point.dt <- rbind(
  data.table(grid.dt[is.finite(value)][which.min(value), .(p1,p2)], point="min"),
  data.table(p1=0.953770837623602,p2=0.0462281623763985,point="last"))
ggplot()+
  theme_bw()+
  scale_color_manual(values=c(
    satisfied="white",
    violated="black"),
    guide=guide_legend(override.aes=list(fill=NA)))+
  geom_tile(aes(
    p1, p2,
    color=constraints,
    fill=log(-value)),
    data=grid.dt)+
  geom_point(aes(
    p1,p2,shape=point),
    size=5,
    data=point.dt)+
  scale_shape_manual(values=c(min=1, last=2))+
  scale_fill_gradient(low="white", high="red")+
  coord_equal()


set.seed(0)
it.dt.list <- list()
for(i in 1:100){
  print(i)  
  startp <- c(1/3,1/3) + rnorm(2)*1e-6 
  cat(paste("startp = c(",paste(startp,collapse=","),")\n"))
  out <- constrOptimMine(
    theta = startp,
    f     = objective,
    grad  = gradfun,
    ui    = UI,
    ci    = CI
  )
  print(out$par)
  print(UI %*% out$par - CI)
}

Output:

...
objective(c( 0.953770837623602,0.0462281623763985 ) = -16.8896756494232 
Erreur dans optim(theta.old, fun, gradient, control = control, method = method,  (depuis bug.R#55) : 
  la valeur initiale dans 'vmin' n'est pas finie
De plus : Messages d'avis :
1: Dans log(1 - p1 - p2) : Production de NaN
2: Dans log(1 - p[1] - p[2]) : Production de NaN

[leopoldguyot]

Another fix would be to add (as mentioned in the initial message) a new parameter tol.contraint. The user will be able to specify a tolerance value to suppress any technical noise that would give a negative value instead of a zero.

Pro:

• Allows to use the "BFGS" method in all cases

Con:

• Needs a new parameter to the function

R <- function(theta, theta.old, ...) {
            ui.theta <- ui%*%theta
            gi <-  ui.theta-ci
            gi[-tol.constraint<gi & gi < 0] <- 0 # suppress noise
            if (any(gi<0)) return(NaN)
            print(gi)
            gi.old <- ui%*%theta.old-ci
            bar <- sum(gi.old*log(gi)-ui.theta)
            if (!is.finite(bar)) bar <-  -Inf
            f(theta, ...) -mu*bar
        }

[tdhock]

@nashjc can you help with this?
What do you think is the best way to deal with this issue?
It involves an objective function with strong gradients near the min.
We wonder if we should add a warning to tell the user to use Nelder-Mead?
Or if we should add some argument to constrOptim(), which would be a numerical threshold to specify the tolerance for verification of constraints: ui %*% theta - ci + tolerance >= 0 ?

[tdhock]

It seems like the solution does not change much between the last few outer iterations, so I wonder if we can solve the problem by playing with the outer iteration parameters, defaults are: outer.iterations = 100, outer.eps = 1e-05,
Here is output that shows the solution as a function of outer iterations.

  for(i in seq_len(outer.iterations)) {
    cat("OUTER ITERATION", i, "\n")
...
> out <- constrOptimMine(
+ theta = startp,
+ f     = objective,
+ grad  = gradfun,
+ ui    = UI,
+ ci    = CI
+ )
objective(c( 0.333334596287618,0.333333007099973 ) = -1.79176466317574 
objective(c( 0.333334596287618,0.333333007099973 ) = -1.79176466317574 
OUTER ITERATION 1 
objective(c( 0.333334596287618,0.333333007099973 ) = -1.79176466317574 
gradfun(c( 0.333334596287618,0.333333007099973 ) = c( -4.50000632289289,-1.50000338678977 )
objective(c( 0.513334849203334,0.393333142571563 ) = -3.20671155050333 
gradfun(c( 0.513334849203334,0.393333142571563 ) = c( -11.817377397968,-9.27500328359227 )
objective(c( 0.532242185496677,0.408172711978915 ) = -3.65498045347121 
gradfun(c( 0.532242185496677,0.408172711978915 ) = c( -17.8460788995925,-15.3961356019336 )
objective(c( 0.560795116853163,0.432805764517313 ) = -5.88264156389459 
gradfun(c( 0.560795116853163,0.432805764517313 ) = c( -157.277961128668,-154.967455858214 )
objective(c( 0.562808208830554,0.434789282252229 ) = -6.86173007978612 
gradfun(c( 0.562808208830554,0.434789282252229 ) = c( -417.233953923582,-414.933989226154 )
objective(c( 0.56387629237121,0.435851477927903 ) = -9.03904583335176 
gradfun(c( 0.56387629237121,0.435851477927903 ) = c( -3674.36874290724,-3672.07438335899 )
objective(c( 0.563951540927358,0.435926679495923 ) = -9.84345732465372 
gradfun(c( 0.563951540927358,0.435926679495923 ) = c( -8212.55765012398,-8210.26368637494 )
objective(c( 0.563985178433308,0.435960307605868 ) = -10.6472032891283 
gradfun(c( 0.563985178433308,0.435960307605868 ) = c( -18344.9248930452,-18342.6311062427 )
objective(c( 0.56400020608556,0.435975333379064 ) = -11.4485946913559 
gradfun(c( 0.56400020608556,0.435975333379064 ) = c( -40883.1796104194,-40880.8859026716 )
objective(c( 0.56400690415293,0.435982031070636 ) = -12.2418869742979 
gradfun(c( 0.56400690415293,0.435982031070636 ) = c( -90377.8825465659,-90375.5888740548 )
objective(c( 0.564009865546864,0.435984992389412 ) = -13.0081983717282 
gradfun(c( 0.564009865546864,0.435984992389412 ) = c( -194475.447187286,-194473.153530354 )
objective(c( 0.564011140008416,0.435986266835932 ) = -13.6927769662221 
gradfun(c( 0.564011140008416,0.435986266835932 ) = c( -385631.534350065,-385629.240699837 )
objective(c( 0.564011645435957,0.435986772260466 ) = -14.1867706051632 
gradfun(c( 0.564011645435957,0.435986772260466 ) = c( -631990.976177834,-631988.682530266 )
objective(c( 0.564011811093593,0.435986937917501 ) = -14.4217179389597 
gradfun(c( 0.564011811093593,0.435986937917501 ) = c( -799368.600738668,-799366.307091971 )
objective(c( 0.564011852996567,0.435986979820355 ) = -14.4910590688428 
gradfun(c( 0.564011852996567,0.435986979820355 ) = c( -856764.620218293,-856762.326571817 )
objective(c( 0.564011897905254,0.435987024728922 ) = -14.5711332667773 
gradfun(c( 0.564011897905254,0.435987024728922 ) = c( -928190.829289087,-928188.535642846 )
objective(c( 0.564011907633523,0.435987034457167 ) = -14.5893576614652 
gradfun(c( 0.564011907633523,0.435987034457167 ) = c( -945261.609606159,-945259.31595997 )
objective(c( 0.564011917539273,0.435987044362893 ) = -14.6082622385021 
gradfun(c( 0.564011917539273,0.435987044362893 ) = c( -963301.343914796,-963299.050268659 )
objective(c( 0.564011927631046,0.435987054454641 ) = -14.6278965251641 
gradfun(c( 0.564011927631046,0.435987054454641 ) = c( -982401.96178309,-982399.668137007 )
objective(c( 0.56401192968739,0.435987056510981 ) = -14.6319450127293 
gradfun(c( 0.56401192968739,0.435987056510981 ) = c( -986387.262231606,-986384.968585534 )
objective(c( 0.564011931750929,0.435987058574515 ) = -14.6360242092018 
gradfun(c( 0.564011931750929,0.435987058574515 ) = c( -990419.144020106,-990416.850374044 )
objective(c( 0.564011933820763,0.435987060644344 ) = -14.640132633919 
gradfun(c( 0.564011933820763,0.435987060644344 ) = c( -994496.573176642,-994494.279530591 )
objective(c( 0.564011935893743,0.43598706271732 ) = -14.6442642922786 
gradfun(c( 0.564011935893743,0.43598706271732 ) = c( -998613.989753819,-998611.696107779 )
objective(c( 0.565846838473174,0.434152160604109 ) = -14.6489474717412 
gradfun(c( 0.565846838473174,0.434152160604109 ) = c( -999079.134023467,-999076.830683553 )
objective(c( 0.565846838877066,0.434152161008001 ) = -14.6497548362884 
gradfun(c( 0.565846838877066,0.434152161008001 ) = c( -999886.080117275,-999883.776777364 )
objective(c( 0.611913279656316,0.388085720326539 ) = -14.7620214489732 
gradfun(c( 0.611913279656316,0.388085720326539 ) = c( -999983.856091194,-999981.279340885 )
objective(c( 0.515216977318553,0.291389154132762 ) = -2.66120221073963 
objective(c( 0.592574019188764,0.368746407087784 ) = -4.21064224837852 
objective(c( 0.608045427562806,0.384217857678788 ) = -5.81055688246542 
objective(c( 0.611139709237614,0.387312147796989 ) = -7.41767413820598 
objective(c( 0.611758565572576,0.387931005820629 ) = -9.0241739112573 
objective(c( 0.611882336839568,0.388054777425357 ) = -10.62073709003 
objective(c( 0.611907091092967,0.388079531746303 ) = -12.1684933326916 
objective(c( 0.611912041943646,0.388084482610492 ) = -13.5163165306671 
objective(c( 0.611913032113782,0.38808547278333 ) = -14.3598437101485 
objective(c( 0.61191323014781,0.388085670817898 ) = -14.6676067428626 
objective(c( 0.611913269754615,0.388085710424811 ) = -14.742411894827 
objective(c( 0.611913277675976,0.388085718346194 ) = -14.758068654726 
objective(c( 0.611913279260248,0.38808571993047 ) = -14.7612296391571 
objective(c( 0.611913279577103,0.388085720247326 ) = -14.761863036845 
objective(c( 0.611913279640474,0.388085720310697 ) = -14.7619897644845 
gradfun(c( 0.611913279640474,0.388085720310697 ) = c( -999952.172638425,-999949.595888116 )
objective(c( 0.869585168855017,0.130413831078691 ) = -15.8524858365984 
gradfun(c( 0.869585168855017,0.130413831078691 ) = c( -999934.712339167,-999927.044440818 )
objective(c( 0.869585168863357,0.130413831087031 ) = -15.852502515714 
gradfun(c( 0.869585168863357,0.130413831087031 ) = c( -999951.390535542,-999943.722637194 )
objective(c( 0.900255893193594,0.0997431067574921 ) = -16.1206179764138 
gradfun(c( 0.900255893193594,0.0997431067574921 ) = c( -999952.088908038,-999942.06315255 )
objective(c( 0.900255893198938,0.0997431067628351 ) = -16.1206286621073 
gradfun(c( 0.900255893198938,0.0997431067628351 ) = c( -999962.774178819,-999952.748423331 )
objective(c( 0.940357726491279,0.0596412734707266 ) = -16.6348790008358 
gradfun(c( 0.940357726491279,0.0596412734707266 ) = c( -999963.006947089,-999946.240034818 )
objective(c( 0.940357726501567,0.0596412734810136 ) = -16.6348995756065 
objective(c( 0.940357726493337,0.059641273472784 ) = -16.6348831157547 
gradfun(c( 0.940357726493337,0.059641273472784 ) = c( -999967.121748489,-999950.354836219 )
objective(c( 0.953770837606663,0.0462281623594607 ) = -16.889641773796 
gradfun(c( 0.953770837606663,0.0462281623594607 ) = c( -999967.124711196,-999945.492875774 )
objective(c( 0.953770837613376,0.0462281623661648 ) = -16.8896551904214 
objective(c( 0.953770837608005,0.0462281623608016 ) = -16.8896444570844 
objective(c( 0.953770837606931,0.0462281623597289 ) = -16.8896423104961 
objective(c( 0.953770837606931,0.0462281623597289 ) = -16.8896423104961 
OUTER ITERATION 2 
objective(c( 0.953770837606931,0.0462281623597289 ) = -16.8896423104961 
gradfun(c( 0.953770837606931,0.0462281623597289 ) = c( -999967.661398565,-999946.029563143 )
objective(c( 0.953770837610287,0.0462281623630842 ) = -16.8896490208056 
gradfun(c( 0.953770837610287,0.0462281623630842 ) = c( -999974.371572696,-999952.739737276 )
objective(c( 0.953770837613642,0.0462281623664395 ) = -16.8896557311809 
gradfun(c( 0.953770837613642,0.0462281623664395 ) = c( -999981.0818577,-999959.450022282 )
objective(c( 0.953770837616997,0.0462281623697948 ) = -16.88966244174 
gradfun(c( 0.953770837616997,0.0462281623697948 ) = c( -999987.792371537,-999966.16053612 )
objective(c( 0.953770837620353,0.0462281623731501 ) = -16.8896691523649 
gradfun(c( 0.953770837620353,0.0462281623731501 ) = c( -999994.502996254,-999972.871160839 )
objective(c( 0.953770837621024,0.0462281623738212 ) = -16.8896704945411 
gradfun(c( 0.953770837621024,0.0462281623738212 ) = c( -999995.845177802,-999974.213342387 )
objective(c( 0.953770837621695,0.0462281623744922 ) = -16.8896718367261 
gradfun(c( 0.953770837621695,0.0462281623744922 ) = c( -999997.187369892,-999975.555534477 )
objective(c( 0.953770837622366,0.0462281623751633 ) = -16.8896731789128 
gradfun(c( 0.953770837622366,0.0462281623751633 ) = c( -999998.529565584,-999976.89773017 )
objective(c( 0.953770837623037,0.0462281623758344 ) = -16.8896745211013 
gradfun(c( 0.953770837623037,0.0462281623758344 ) = c( -999999.87176488,-999978.239929466 )
objective(c( 0.953770837623172,0.0462281623759686 ) = -16.8896747895365 
gradfun(c( 0.953770837623172,0.0462281623759686 ) = c( -1000000.1402024,-999978.508366982 )
objective(c( 0.953770837623306,0.0462281623761028 ) = -16.8896750579717 
gradfun(c( 0.953770837623306,0.0462281623761028 ) = c( -1000000.40864006,-999978.776804642 )
objective(c( 0.95377083762344,0.046228162376237 ) = -16.889675326407 
gradfun(c( 0.95377083762344,0.046228162376237 ) = c( -1000000.67707786,-999979.045242446 )
objective(c( 0.953770837623574,0.0462281623763712 ) = -16.8896755948424 
gradfun(c( 0.953770837623574,0.0462281623763712 ) = c( -1000000.94551581,-999979.313680394 )
objective(c( 0.953770837623601,0.0462281623763981 ) = -16.8896756485489 
objective(c( 0.953770837623601,0.0462281623763981 ) = -16.8896756485489 
OUTER ITERATION 3 
objective(c( 0.953770837623601,0.0462281623763981 ) = -16.8896756485489 
gradfun(c( 0.953770837623601,0.0462281623763981 ) = c( -1000000.99922284,-999979.36738743 )
objective(c( 0.953770837623601,0.0462281623763983 ) = -16.889675648986 
OUTER ITERATION 4 
objective(c( 0.953770837623601,0.0462281623763983 ) = -16.889675648986 
gradfun(c( 0.953770837623601,0.0462281623763983 ) = c( -1000000.99965999,-999979.367824581 )
objective(c( 0.953770837623602,0.0462281623763985 ) = -16.8896756494232 
OUTER ITERATION 5 
Erreur dans optim(theta.old, fun, gradient, control = control, method = method,  (depuis bug.R#53) : 
  la valeur initiale dans 'vmin' n'est pas finie

[mmaechler]

Did anyone try (with some care) to follow my original suggestion (in R's bugzilla, cited by @hturner, top of this page) ?

[tdhock]

We did try your suggestion, but we seem to have found a better solution, involving a warning an an early exit (more details to come).

[leopoldguyot]

Also the "removing noise" method was still giving errors with other optimisation method like CG.