Equality constraint with COBYLA

[fangzhou-xie]

Hi, the NLopt documentation mentions that "Only some of the NLopt algorithms (AUGLAG, SLSQP, COBYLA, and ISRES) currently support nonlinear equality constraints".

However, the following example will run into error:

f <- function(x) {
  return(sum(x^2))
}

eq <- function(x) {
  return(sum(x)-1)
}

sol <- nloptr::nloptr(
  x0 = c(.1, .4),
  eval_f = f,
  lb = c(0, -Inf),
  ub = c(Inf, Inf),
  eval_g_eq = eq,
  opts = list("algorithm" = "NLOPT_LN_COBYLA")
)
#> Warning in nloptr.add.default.options(opts.user = opts, x0 = x0,
#> num_constraints_ineq = num_constraints_ineq, : No termination criterium
#> specified, using default (relative x-tolerance = 1e-04)
#> Error in is.nloptr(ret): If you want to use equality constraints, then you should use one of these algorithms 
# NLOPT_LD_AUGLAG, NLOPT_LN_AUGLAG, NLOPT_LD_AUGLAG_EQ, NLOPT_LN_AUGLAG_EQ, NLOPT_GN_ISRES, NLOPT_LD_SLSQP

^{Created on 2023-03-14 with reprex v2.0.2}

I am not sure why this is the case and I greatly appreciate any help. Thank you very much!

[fangzhou-xie]

For my naive example, I could pack the COBYLA algorithm as the local algo for AUGLAG to make it work. But I am still curious why directly using COBYLA wouldn't work.

f <- function(x) {
  return(sum(x^2))
}

eq <- function(x) {
  return(sum(x)-1)
}

sol <- nloptr::nloptr(
  x0 = c(.1, .4),
  eval_f = f,
  lb = c(0, -Inf),
  ub = c(Inf, Inf),
  eval_g_eq = eq,
  opts = list(
    "algorithm" = "NLOPT_LN_AUGLAG_EQ",
    "xtol_rel" = 1.0e-8,
    "maxeval"    = 1000,
    "local_opts" = list(
      "algorithm" = "NLOPT_LN_COBYLA",
      "maxeval"    = 1000,
      "xtol_rel" = 1.0e-8
    )
  )
)

sol$solution
#> [1] 0.4999605 0.4999605

^{Created on 2023-03-14 with reprex v2.0.2}

[aadler]

It's kicked out by lines 193–203 in is.nlopt. While it isn't in the help there, it is in the help for cobyla. Specifically "COBYLA supports equality constraints by transforming them into two inequality constraints. As this does not give full satisfaction with the implementation in NLOPT, it has not been made available here." I presume @jyypma and @hwborchers agreed not to implement it using nlopt as well.

[astamm]

This is worth looking into. I'll take a look soon. If it is possible to use nonlinear equality constraints in COBYLA from nlopt, I do not see why this would not be supported by nloptr.

[jyypma]

The NLopt docs (https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#cobyla-constrained-optimization-by-linear-approximations) mention the following: "The underlying COBYLA code only supports inequality constraints. Equality constraints are automatically transformed into pairs https://nlopt.readthedocs.io/en/latest/NLopt_Introduction#equality-constraints of inequality constraints, which in the case of this algorithm seems not to cause problems."

When looking at the Cobyla code in NLopt, it does seem to add two inequality constraints when equality constraints are added. I'm not sure why I excluded equality constraints for Cobyla in nloptr, but I'm assuming that this feature in NLopt was not yet available when I implemented the first version of the wrapper. As a workaround, it's quite straightforward to convert the equality constraint to an inequality constraint and pass that to the Cobyla algorithm. That could solve the problem in the original post in the short-run. An example of this work-around is in the code below.

As @aadler noted the more permanent fix can be very simple, by adding "NLOPT_LN_COBYLA" to the list of admissible algorithms in lines 193-203 of the is.nloptr function. The example at the bottom of this message, shows that adapting is.nloptr (and temporariliy by-passing all the checks) seems sufficient to let the Cobyla algorithm use equality constraints.

Code example:

library('nloptr')

f <- function(x) {
  return(sum(x^2))
}

eq <- function(x) {
  return(sum(x)-1)
}

# Re-write equality constraint using
# two inequality constraints.
eq_as_ineq <- function(x) {
  eqx <- eq(x)
  return(c(eqx, -eqx))
}

# Solve problem with equality constraint re-written as
# inequality constraint as a work-around.
nloptr::nloptr(
  x0 = c(.1, .4),
  eval_f = f,
  lb = c(0, -Inf),
  ub = c(Inf, Inf),
  eval_g_ineq = eq_as_ineq,
  opts = list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-6))

# Ugly hack to temporarliy by-pass running is.nloptr.
is.nloptr_original <- nloptr::is.nloptr
assignInNamespace('is.nloptr', function(x) { return(TRUE) }, ns = "nloptr")

# Solve problem using equality constraints.
nloptr::nloptr(
  x0 = c(.1, .4),
  eval_f = f,
  lb = c(0, -Inf),
  ub = c(Inf, Inf),
  eval_g_eq = eq,
  opts = list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-6))

# Return is.nloptr to original implementation.
assignInNamespace('is.nloptr', is.nloptr_original, ns = "nloptr")

[astamm]

Thanks @jyypma and @aadler. I will make the fix in is.nloptr.