NLopt - Same complex clalculations in objetive function and nonlinear constraints

Usage Performance
#1

Hey guys,
I need a hint on how to effectivley solve a nonlinear optimization problem with NLopt package when my objective function contains some complex calculations and my nonlinear constraint contains the same complex calculations, something like the following:

using NLopt

function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
        # gradient calculation
    end
    Y = SomeComplexCalculations( x)
    return UseYForObjective( Y)
end

function myconstraint(x::Vector, grad::Vector)
    if length(grad) > 0
        # gradient calculation
    end
    Y = SomeComplexCalculations( x)
    return UseYForConstraint( Y)
end

Of course, I could do it as shown above, but that looks very unefficient to me.
I also thought about making Y a global variable, but I am not sure if thats possible at all or I might get into trouble during the optimization process.
So I wonder if there is another way of giving data calculated by the objective function to the constraints.
Thanks for any suggestion or comment in advance!

2 Likes
#2

Yes, it looks like it is the way to go.

#3

Which one is the way to go accroding to you opinion?
The approach that I posted the code example for or the approach with global variables?

I am asking because I am unsure whether global variables will really make it at all. I mean I am not familiar to the internals of NLopt, but I see a potential problem here if I decide to call SomeComplexCalculations( x) in myfunc only since this implies that myconstraint is just called with the same vector x until myfunc is called again. However if myconstrained was called more than once with two different x during the optimiziation process and myfunc was only called after a suitable x that fullfills the condition was found, I would get into trouble.

#4

I am not an expert but I think the way to go is caching the result of the last calculation in two variables, x_cached and y_cached. Then before performing the expensive computation, you check whether the vector x is the same as x_cached for which you have already computed Y and stored it into y_cached.

using NLopt

function SomeComplexCalculations(x)
    println("complex calculation with $x")
    return x.^2
end


function myoptim()
    function myfunc(x::Vector, grad::Vector)
        if length(grad) > 0
            grad[1] = 0
            grad[2] = 0.5/sqrt(x[2])
        end
        if x_cached != x
            println("objfunc: not cached      $x")
            x_cached[:] = x
            y_cached[:] = SomeComplexCalculations(x)
        else
            println("objfunc: using cached $x")
        end
        return sqrt(x[2])
    end

    function myconstraint(x::Vector, grad::Vector, a, b)
        if length(grad) > 0
            grad[1] = 3a * (a*x[1] + b)^2
            grad[2] = -1
        end
        if x_cached != x
            println("constraint: not cached $x")
            x_cached[:] = x
            y_cached[:] = SomeComplexCalculations(x)
        else
            println("constraint: using cached $x")
        end
        (a*x[1] + b)^3 - x[2]
    end

    opt = Opt(:LD_MMA, 2)
    x_cached = Vector{Float64}(undef, 2)
    y_cached = Vector{Float64}(undef, 2)

    opt.lower_bounds = [-Inf, 0.]
    opt.xtol_rel = 1e-4

    opt.min_objective = myfunc
    inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)
    inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)

    (minf,minx,ret) = optimize(opt, [1.234, 5.678])
    numevals = opt.numevals # the number of function evaluations
    println("got $minf at $minx after $numevals iterations (returned $ret)")
end

I was surprised to see that the objective function is computed first and then the constraints, but this may have to do with the choice of the algorithm.

1 Like
#5

That looks like a smart solution.
Thanks for your support!

#6

IIRC, NLopt always evaluates the objective before the constraints.

1 Like
#7

In my ignorance, I assumed that one may want to check whether the solution is feasible before computing the objective value.

#8

Even if the point is infeasible, the objective value is still important to determine the next optimization step.

1 Like
#9

The risk is to sample a point that does not belong to the objective function domain, e.g. a negative number when f(x) = sqrt(x).

2 Likes
#10

NLopt guarantees that the function will never be evaluated outside the feasible set for box constraints, but not for arbitrary nonlinear constraints. (This is true for most nonlinear-programming algorithms that I’m aware of.) So if you have \sqrt{x}, you should use a box constraint x \ge 0 and you will be fine.

(For general constraints, this may not even be possible, e.g. you may have a pair of inequality constraints h(x) \le 0 and h(x) \ge 0, which means the feasible set is in a nonlinear manifold h(x)=0, and it is not practical to stay exactly on such a manifold in general.)

3 Likes
#11

I think you are right, but I am wondering if some simple heuristic could remedy this. Eg if the feasible set D is convex, x \in D, and y \notin D, one could find a z = \alpha x + (1 - \alpha) y \in D using bisection. (This is just a random idea, maybe it has been tried and discarded).

#12

One could certainly devise methods for this. Assuming that the feasible set has a nonempty interior (no equality constraints) and one starts with a feasible point, you could simply shrink the size of the trust region (or whatever method is used to limit step sizes) until the next step is feasible.

However, this may come at a price in convergence rate, since information from function evaluations outside the feasible set may be used by an algorithm to accelerate convergence, e.g. the CCSA algorithm uses this information. Moreover, other algorithms like SLSQP don’t assume that the feasible set has a nonempty interior. (It would be pretty easy to add this as a flag for CCSA, however.) More radical modifications might be needed for penalty methods such as the augmented-Lagrangian approach.

Typically, however, if you have a constraint that must be satisfied in order to evaluate your objective, I would suggest that you impose it as a box constraint or to build it into your parameterization of the problem. (e.g. if you have \sqrt{x}, parameterize x as y^2 so that it can never be negative.) This gives you maximum flexibility in your choice of algorithms.

1 Like
#13

I am not sure how this generalize, when for instance the objective value is f(\sqrt{x}) + g(\sqrt{y}) + h(\sqrt{x-y}), where the domain is the triangle with vertices (0,0), (1,0), (0.1). A change of variables does not help, I guess.

EDIT: The last term should have been h(\sqrt{1-x-y}) to form a triangle…

#14

For your example, I think the domain is x \ge y \ge 0, which can be transformed to (y, z) \ge 0 with x = y + z.

But yes, simplexes (and generally polyhedra) are tricky to transform to a box because you get a singularity at one or more vertices. If you can exclude these vertices (eg if function is continuous, just evaluate it at all vertices and do the transformation so that the singularities are at the ones which cannot be optimal).

Admittedly, these transformations do look a bit tedious and it would be nice to have a method that is clean and pre-packaged and handles everything automatically, but since there isn’t one, it is generally worthwhile to think about transformations.

1 Like