Adding Equality constraint breaks variables / gradients computation

[JeremyBois]

Hi,

Description

Adding an equality constraint silently break the optimization. Independent variables becomes NaN for all iteration steps (excepting the first one) breaking gradient computation too.

Notes

• I can add inequality constraints
• I test writing equality constraints with delegate, lambda, class method --> Same result
• I'm using latest version of NLoptNet
• Tried with NLopt 2.6.2 and 2.6.1 --> Same result

Minimal breaking example

Below the minimal code that can be used to reproduce the behaviour.

public void TestSLSQP()
{
    var gradientsHistory = new List<double[]>();
    var variablesHistory = new List<double[]>();
    double? finalScore;
    NloptResult result;

    using (var solver = new NLoptSolver(NLoptAlgorithm.LD_SLSQP, 2, maximumIterations: 20))
    {
        var initialValue = new[] { 1.234, 5.678 };
        solver.SetLowerBounds(new[] { double.NegativeInfinity, 0.0000000000001 });

        int count = 0;
        solver.SetMinObjective((variables, gradient) =>
        {
            if (gradient != null)
            {
                gradient[0] = 0.0;
                gradient[1] = variables[1];
            }
            count++;

            var iterationGradients = new double[2];
            Array.Copy(gradient, iterationGradients, 2);
            gradientsHistory.Add(iterationGradients);

            var iterationValues = new double[2];
            Array.Copy(variables, iterationValues, 2);
            variablesHistory.Add(iterationValues);

            return variables[1];
        });

        // @BUG This is the line that broke everything
        solver.AddEqualZeroConstraint(TestEq);
        // @BUG This is the line that broke everything

        result = solver.Optimize(initialValue, out finalScore);
    }
}

private double TestEq(double[] values, double[] grads)
{
    return values.Length + grads.Length;
}

[BrannonKing]

You're certain that variables[1] is never zero?

Also, without looking into it, I'm unsure why @Martin1994 made the gradient handling on the equality constraint different from the inequality mechanism. See

NLoptNet/NLoptNet/NLoptSolver.cs

Line 249 in 72be440

if (gradient != null)

[JeremyBois]

Actually this is only a dummy example to isolated the issue I have with the problem I want to optimize.
In the complete optimization I use a two-point estimation for gradient calculation so this is not a division by zero error.
I test again the dummy example provided without division and the problem still persist.

I cannot see any differences between implementation for both equality and inequality (at least on the C# side). There is two different interfaces depending on the chosen optimizer (some required gradient to be computed where others not).

[Martin1994]

I tried to reproduce this with libnlopt. Such an equality constraint neither works with C, but behaves differently. Instead of using up 20 iterations, SLSQP actually exits right after the first iteration with NLOPT_ROUNDOFF_LIMITED error.

#include <stdio.h>
#include <nlopt.h>

int count = 0;

double minObj(unsigned int n, const double* variables, double* gradient, void* data) {
    if (gradient)
    {
        gradient[0] = 0.0;
        gradient[1] = variables[1];
    }

    printf("> Iteration #%d: x = [%f, %f], gradient = [%f, %f]\n", count, variables[0], variables[1], gradient[0], gradient[1]);

    count++;

    return variables[1];
}

double testEq(unsigned int n, const double* variables, double* gradient, void* data) {
    return 4;
}

int main() {
    double zero = 0;
    double finalScore;
    nlopt_result result;

    nlopt_opt opt = nlopt_create(NLOPT_LD_SLSQP, 2);
    if (!opt) {
        printf("Unable to initialize the algorithm.\n");
        return 1;
    }

    result = nlopt_set_xtol_rel(opt, 0.0001);
    if (result != NLOPT_SUCCESS) {
        printf("Unable to set primary tolerance. Result: %d\n", result);
        return 1;
    }

    result = nlopt_set_maxeval(opt, 20);
    if (result != NLOPT_SUCCESS) {
        printf("Unable to set max evaluation. Result: %d\n", result);
        return 1;
    }

    double initialValue[] = { 1.234, 5.678 };

    double minimums[] = { -1.0 / 0.0, 0.0000000000001 };
    result = nlopt_set_lower_bounds(opt, minimums);
    if (result != NLOPT_SUCCESS) {
        printf("Unable to set lower bounds. Result: %d\n", result);
        return 1;
    }

    result = nlopt_set_min_objective(opt, minObj, &zero);
    if (result != NLOPT_SUCCESS) {
        printf("Unable to set the objective function. Result: %d\n", result);
    }

    // @BUG This is the line that broke everything
    nlopt_add_equality_constraint(opt, testEq, &zero, 0.001);
    // @BUG This is the line that broke everything

    result = nlopt_optimize(opt, initialValue, &finalScore);

    nlopt_destroy(opt);

    printf("result: %d\n", result);
}

Output:

> Iteration #0: x = [1.234000, 5.678000], gradient = [0.000000, 5.678000]
result: -4

[JeremyBois]

A more complete code example using numerical two point differentiation for gradient calculation and a equality constraint that is always respected (always 0.0).

• Works without the equality constraint
• Breaks with the equality constraint with two different errors

@Martin1994, I also get the ROUNDOFF_LIMITED error code alternating with MAXEVAL_REACHED (see outputs below)
Can you confirm you get the same behaviour with this code in C++. If this is the case, maybe we could open an issue directly on the NLopt repository because it will not be related to marshalling issues.

Source code (C#)

public void TestSLSQP()
{
    var gradientsHistory = new List<double[]>();
    var variablesHistory = new List<double[]>();
    double? finalScore;
    NloptResult result;
    int count = 0;

    using (var solver = new NLoptSolver(NLoptAlgorithm.LD_SLSQP, 2, maximumIterations: 20))
    {
        var initialValue = new[] { 1.234, 5.678 };
        solver.SetLowerBounds(new[] { 1.0, 1.0 });
        solver.SetUpperBounds(new[] { 10.0, 10.0 });

        solver.SetMinObjective((variables, gradient) =>
        {
            var h = 0.00001;
            var f0 = Math.Pow(variables[1] * variables[0], 2);
            var fi0 = Math.Pow(variables[1] * (variables[0] + h), 2);
            var fi1 = Math.Pow((variables[1] + h) * variables[0], 2);

            if (gradient != null)
            {
                // two point numerical differentiation
                gradient[0] = (fi0 - f0) / h;
                gradient[1] = (fi1 - f0) / h;
            }
            count++;

            var iterationGradients = new double[2];
            Array.Copy(gradient, iterationGradients, 2);
            gradientsHistory.Add(iterationGradients);

            var iterationValues = new double[2];
            Array.Copy(variables, iterationValues, 2);
            variablesHistory.Add(iterationValues);

            return f0;
        });

        // @BUG This is the line that broke everything
        solver.AddEqualZeroConstraint(TestEq);
        // @BUG This is the line that broke everything

        result = solver.Optimize(initialValue, out finalScore);
        Debug.Log($"Target value = {finalScore.Value} (nloptCode={result})");
    }

    for (int i = 0; i < count; i++)
    {
        Debug.Log($"#{i} --> variables=({variablesHistory[i][0]}, {variablesHistory[i][1]}) - gradient=({gradientsHistory[i][0]}, {gradientsHistory[i][1]})");
    }
}

private double TestEq(double[] values, double[] grads)
{
    return 0.0;
}

Outputs

// Without equality constraint
Target value = 1 (nloptCode=XTOL_REACHED)
#0 --> variables=(1.234, 5.678) - gradient=(79.5678625088669, 17.2924323628365)
#1 --> variables=(1, 1) - gradient=(2.00001000001393, 2.00001000001393)

// With equality constraint (happens from time to time ...)
Target value = 49.093172249104 (nloptCode=MAXEVAL_REACHED)
#0 --> variables=(1.234, 5.678) - gradient=(79.5678625088669, 17.2924323628365)
#1 --> variables=(NaN, NaN) - gradient=(NaN, NaN)
...
#19 --> variables=(NaN, NaN) - gradient=(NaN, NaN)

// With equality constraint (happens from time to time ...)
Target value = 49.093172249104 (nloptCode=ROUNDOFF_LIMITED)
#0 --> variables=(1.234, 5.678) - gradient=(79.5678625088669, 17.2924323628365)

[JeremyBois]

Open an issue directly to nlopt repository. See stevengj/nlopt#400.