PatWie / CppNumericalSolvers

a lightweight header-only C++17 library of numerical optimization methods for nonlinear functions based on Eigen
MIT License
878 stars 201 forks source link

L-BFGS-B doesn't work with simple example #116

Closed dandoh closed 5 years ago

dandoh commented 5 years ago

I tried to minimize norm2square(x[4]) given 20 <= x <= 40.

#include <iostream>
#include "../../include/cppoptlib/meta.h"
#include "../../include/cppoptlib/boundedproblem.h"
#include "../../include/cppoptlib/solver/lbfgsbsolver.h"

namespace cppoptlib {

template<typename T>
class NonNegativeLeastSquares : public BoundedProblem<T> {
  public:
    using Superclass = BoundedProblem<T>;
    using typename Superclass::TVector;
    using TMatrix = typename Superclass::THessian;

  public:
    NonNegativeLeastSquares(int dim) :
        Superclass(dim) {}

    T value(const TVector &beta) {
        return beta.dot(beta);
    }

    void gradient(const TVector &beta, TVector &grad) {
        grad = beta * 2;
    }
};

}
int main(int argc, char const *argv[]) {

    const size_t DIM = 4;
    typedef double T;
    typedef cppoptlib::NonNegativeLeastSquares<T> TNNLS;
    typedef typename TNNLS::TVector TVector;

    // perform non-negative least squares
    TNNLS f(4);
    f.setLowerBound(TVector::Ones(DIM) * 20);
    f.setUpperBound(TVector::Ones(DIM) * 40);
    // create initial guess
    TVector beta = TVector::Ones(DIM) * 21;
    cppoptlib::LbfgsbSolver<TNNLS> solver;
    solver.minimize(f, beta);
    std::cout << "final b = " << beta.transpose() << "\tloss:" << f(beta) << std::endl;

    return 0;
}

But I got the result:

final b = 0 0 0 0  loss:0

It should be:

final b = 20 20 20 20   loss:1600
xbh2015 commented 5 years ago

I got this problem too and I just add the bound check after line search step which returns right answer.In paper, I found this sentence "We then truncate the path toward the solution so as to satisfy the bounds".

PatWie commented 5 years ago

Could you please elaborate what you exactly changed?

xbh2015 commented 5 years ago
  const Scalar rate = MoreThuente<TProblem, 1>::linesearch(x,  SubspaceMin-x ,  problem, alpha_init);
  // update current guess and function information
  x = x - rate*(x-SubspaceMin);

After these two step, I just check whether x in bound. If any of variable out of bound, I just set the varibale upper or lower bound. I am not sure whether it is right explanation to that sentence. At first, I think equation 5.8 or 5.9 should be considered as the explanation of that sentence.

PatWie commented 5 years ago

This is fixed now. I added your example in src/examples/nonnegls2.cpp