Does not handle infinities correctly about lbfgspp

[yixuan]

Thanks. I'll take a look then.

[mad-s]

Here's a simple MWE that gives fx = inf while the initial value is feasible (f(x0) != inf): https://gist.github.com/mad-s/c097680d0c083c1df3a5c5d6e71c0856

[yixuan]

This has been fixed, but note that your code to compute the gradient needs to be modified. These two lines lack two multipliers a1 and a2:

grad += -mu/c1 * Vector2d(-3*a1*pow(a1*x(0)+b1, 2), 1);
grad += -mu/c2 * Vector2d(-3*a2*pow(a2*x(0)+b2, 2), 1);

The example itself is a bit ill-conditioned, so to get the global optimum you need a better initial value, for example (0.0, 10.0). You also need to increase max_linesearch whenever it is necessary. Below is a working example modified from yours:

#include <Eigen/Core>
#include <iostream>
#include <LBFGS.h>

using namespace std;
using namespace LBFGSpp;
using Eigen::Vector2d;
using Eigen::VectorXd;

int main() {

    LBFGSParam<double> param;
    param.max_iterations = 100;
    param.max_linesearch = 30;
    LBFGSSolver<double> solver(param);
    
    auto fun = [] (const VectorXd &x, VectorXd &grad) -> double {
        double mu = 0.0001;

        double fx = x(1);//sqrt(x(1));
        double a1 = 2;
        double b1 = 0;
        double a2 = -1;
        double b2 = 1;

        // c_i >= 0
        double c1 = x(1) - pow(a1*x(0)+b1, 3);
        double c2 = x(1) - pow(a2*x(0)+b2, 3);
        double c3 = x(1);

        fx += c1 > 0 ? -mu*log(c1) : std::numeric_limits<double>::infinity();
        fx += c2 > 0 ? -mu*log(c2) : std::numeric_limits<double>::infinity();
        fx += c3 > 0 ? -mu*log(c3) : std::numeric_limits<double>::infinity();

        grad = Vector2d(0., 1/* / sqrt(x(1))*/);

        grad += -mu/c1 * Vector2d(-3*a1*pow(a1*x(0)+b1, 2), 1);
        grad += -mu/c2 * Vector2d(-3*a2*pow(a2*x(0)+b2, 2), 1);
        grad += -mu/c3 * Vector2d(0.,1.);

        return fx;
    };

    VectorXd x = Vector2d(0., 10.);
    double fx;
    solver.minimize(fun, x, fx);

    cout << "fx = " << fx << endl;
    cout << "x = " << x << endl;
}

Output is

fx = 0.29831
x = 0.333296
0.296496