stevengj / nlopt

library for nonlinear optimization, wrapping many algorithms for global and local, constrained or unconstrained, optimization
Other
1.86k stars 575 forks source link

Getting NLOPT constraint and objective functions to safely return cached values #492

Open olumide-x opened 1 year ago

olumide-x commented 1 year ago

The problem that I am trying to solve is associated with a large data structure that should be accessed as infrequently as possible. This is why I would like to compute and cache my objective and constraints values as shown in Problem::update() below

struct Problem;

struct ConstraintParams 
{
    ConstraintParams(Problem* p)
        : problem(p)
        , constraintId(0){}
    Problem*   problem;
    int        constraintId;
};

double objectiveFunc(const std::vector<double>& x, std::vector<double>& grad, void* data);
double constraintFunc(const std::vector<double>& x, std::vector<double>& grad, void* data);

struct Problem
{
    // update objective and constraints in one go
    // To be called ONCE per iteration
    void update(const std::vector<double>& x)
    {
        // update objective for the current iteration
        // update constraints for the current iteration
    }

    // Sets up solver.
    void solve()
    {
        int numConstraints = constraint.size();
        nlopt::opt opt(nlopt::LD_SLSQP, 2);
        std::vector<ConstraintParams> constraintParams(numConstraints, this);

        opt.set_min_objective(objectiveFunc, this);

        for(int i = 0; i < numConstraints; ++i)
        {
            // Pair constraint index with this pointer
            constraintParams[i].constraintId = i;
            opt.add_equality_constraint(constraintFunc, &constraintParams[i]);
        }
    }

    std::vector<double> constraint;
    double objective;
};

double objectiveFunc(const std::vector<double>& x, std::vector<double>& grad, void* data)
{
    Problem* problem = (Problem*) data;
    problem->update(x); // Recall: do this once per iteration
    return problem->objective;
}

double constraintFunc(const std::vector<double>& x, std::vector<double>& grad, void* data)
{
    ConstraintParams* constraintParams = (ConstraintParams*) data;
    Problem* problem = constraintParams->problem;
    int index = constraintParams->index;
    return problem->constraint[index];
}

What I am not sure about is whether my assumption that:

  1. The objective function is always called first is valid. (If not then a segfault is the best outcome I can hope for.)
  2. The objective function is only called once per iteration regardless of the dimension of the problem is valid. If the objective is called multiple times per iteration, e.g. in order to approximate the gradients of the objective by finite difference, then the last x that is passed to Problem::update() would be used to compute the constraints values.

An idea that's just popped to mind is to use constraint and objective functions of the form

double func(unsigned n, const double* x, double* grad, void* data);

and current and previous pointers x in order to determine whether Problem::update() should be called. This would solve the problem of needing to know whether the constraint or objective functions gets called first. However I am not sure whether my assumption about whether the same pointer x is passed to all functions is valid.

jschueller commented 1 year ago

I'd say dont rely on the pointers, just cache the input values and the associated objective and constraints values and in your case maybe it will be easier to use add_inequality_mconstraint to return all constraints at once

olumide-x commented 1 year ago

I hadn't thought about void add_equality_mconstraint. Thanks.

Problem though is that my x is typically high dimensional vector which be expensive to cache as well as lookup.

If only I could get my hands on the iteration index.

olumide-x commented 1 year ago

Minor question. The signature of the function passed to add_equality_mconstraint is

void constraintsFunc(unsigned m, double* result, unsigned n, const double* x, double* grad, void* data)

Where m is the number of constraints, n the number of variables. What is the size of grad? My guess is m*n. Is this correct? If so n what order should the derivatives of the constraints be written to grad? The derivatives of constraint 1, the the derivatives of constraint 2 etc.

jschueller commented 1 year ago

yes, per constraint then per variable; this is explained in the doc: https://nlopt.readthedocs.io/en/latest/NLopt_Reference/#vector-valued-constraints