Understanding why COBYLA doesn't find optimum here

New to Julia
1

Hi, I’m trying to understand why the COBYLA algorithm doesn’t find the optimal value in the following MWE. By contrast, IPNewton does work. Am I violating some COBYLA assumption that I’m forgetting about? For the problem with non-binding constraints, uncomment the other definition for prob and comment the current one; this gives similar behaviour.

using Random, Optimization, OptimizationMOI, OptimizationNLopt, OptimizationOptimJL, OptimizationPRIMA, DifferentiationInterface, Ipopt
Random.seed!(42)

function RunTest()
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = [0.5, 0.5]
    _p = [1.0, 100.0]
    
    f(x) = [x[1]^2 + x[2]^2, x[1] * x[2]]
    cons(res, x, p) = (res .= f(x))
    
    optprob = OptimizationFunction(rosenbrock, SecondOrder(AutoForwardDiff(), AutoForwardDiff()), cons = cons)
    #prob = OptimizationProblem(optprob, x0, _p, lb = [0., 0.], ub = [1., 1.], lcons = [-1., -1.0], ucons = [2., 2.0])
    prob = OptimizationProblem(optprob, x0, _p, lb = [0., 0.], ub = [1., 1.], lcons = [-1., -1.0], ucons = [0.8, 2.0])

    opt_vec = [COBYLA(), IPNewton()]
    opt_names = ["COBYLA", "IPNewton"]

    for (opt, opt_name) in zip(opt_vec, opt_names)
        sol = solve(prob, opt)
        print("====\n$opt_name\n$sol")
    end
end

RunTest()

I get for the problem with constraints:

====
COBYLA
retcode: Success
u: [1.58113883008419e-6, 4.743416490252569e-6]
Final objective value:     0.9999968399748375
====
IPNewton
retcode: Success
u: [0.7247018404470709, 0.5242206047577943]
Final objective value:     0.07588358475785834
2

I think it has to do with the way Optimization.jl wraps PRIMA, because at first sight using PRIMA.jl directly seems to lead to fine results (although the API is a bit different, of course…)

3

Can you share the PRIMA script with a side-by-side of how you ran it? That would be good for debugging why it’s different.

4

@ChrisRackauckas can you specify what you mean precisely? In the MWE I wrote up the call I made to PRIMA (through Optimization.jl). I’ve included a version with the cobyla function from PRIMA.jl directly, which does give the right solution:

using Random, Optimization, OptimizationMOI, OptimizationNLopt, OptimizationOptimJL, OptimizationPRIMA, DifferentiationInterface, Ipopt, PRIMA
Random.seed!(42)

function RunTest()
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = [0.5, 0.5]
    _p = [1.0, 100.0]
    
    f(x) = [x[1]^2 + x[2]^2, x[1] * x[2]]
    cons(res, x, p) = (res .= f(x))
    
    optprob = OptimizationFunction(rosenbrock, SecondOrder(AutoForwardDiff(), AutoForwardDiff()), cons = cons)
    #prob = OptimizationProblem(optprob, x0, _p, lb = [0., 0.], ub = [1., 1.], lcons = [-1., -1.0], ucons = [2., 2.0])
    prob = OptimizationProblem(optprob, x0, _p, lb = [0., 0.], ub = [1., 1.], lcons = [-.8, -2.0], ucons = [0.8, 2.0])

    opt_vec = [COBYLA(), IPNewton()]
    opt_names = ["COBYLA", "IPNewton"]

    for (opt, opt_name) in zip(opt_vec, opt_names)
        sol = solve(prob, opt)
        print("====\n$opt_name\n$(sol)")
    end

    x, info = cobyla(x -> rosenbrock(x, _p), x0; nonlinear_ineq = x -> f(x) .- [0.8, 2.])
    print("====\nusing cobyla directly: $x\n$info\n")
end

RunTest()

If I had to guess it may be that PRIMA.jl requires constraints to be of the “f <= 0” type rather than of the “lc <= f <= uc” type, perhaps there is some error in translation here when Optimization.jl wraps COBYLA() and calls cobyla() in PRIMA.jl? I’m not too much aware of the exact implementation here.

====
COBYLA
retcode: Success
u: [1.58113883008419e-6, 4.743416490252569e-6]
Final objective value:     0.9999968399748375
====
IPNewton
retcode: Success
u: [0.7247018404470709, 0.5242206047577943]
Final objective value:     0.07588358475785834
====
using cobyla directly: [0.7247016181219436, 0.5242209192620982]
PRIMA.Info(0.07588358340753011, 365, PRIMA.SMALL_TR_RADIUS, 7.500562615447848e-9, Float64[], [7.500562615447848e-9, -1.6200962515573847])
5

Can you open an issue?

6

Sure, this is now Wrapping cobyla from PRIMA.jl seems to handle constraints incorrectly · Issue #854 · SciML/Optimization.jl · GitHub.