Division by zero in Cauchy.h

[xthan]

Dear Yixuan,

Thanks for sharing the great code. I found that in some rare cases, the variable fpp in Cauchy.h may happen to be zero, which leads to NaN results：

https://github.com/yixuan/LBFGSpp/blob/8dd1b51bdd10cdc67ca17739bfb26cec7a93baa4/include/LBFGSpp/Cauchy.h#L243

[yixuan]

Thanks for the report. I'm a bit busy recently, but I'll look into this issue when I get some time.

[yixuan]

I read the original paper on L-BFGS-B again, and noticed that fpp is guaranteed to be positive mathematically, so in your case there may be some numerical issues. Could you provide a simple example to reproduce the problem?

[xthan]

Here is an example modified from https://github.com/yixuan/LBFGSpp/blob/master/example-quadratic.cpp, which should be able to reproduce the problem during the optimization.

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

using Eigen::VectorXd;
using Eigen::MatrixXd;
using namespace LBFGSpp;

double foo(const VectorXd& x, VectorXd& grad)
{
    const int n = x.size();
    VectorXd d(n);
    for(int i = 0; i < n; i++)
        d[i] = i + 1;
    double f = 20. * (x - d).squaredNorm();
    grad.noalias() = 20. * 2.0 * (x - d);
    std::cout << "fx:" << f << std::endl;

    return f;
}

int main()
{
    const int n = 33;
    LBFGSBParam<double> param;
    LBFGSBSolver<double> solver(param);

    VectorXd lb = VectorXd::Constant(n, -30.0);
    VectorXd ub = VectorXd::Constant(n, 30.0);

    // Initial guess
    VectorXd x = VectorXd::Constant(n, 1.);

    // x will be overwritten to be the best point found
    double fx;
    int niter = solver.minimize(foo, x, fx, lb, ub);

    std::cout << niter << " iterations" << std::endl;
    std::cout << "x = \n" << x.transpose() << std::endl;
    std::cout << "f(x) = " << fx << std::endl;

    return 0;
}

[yixuan]

This issue has been fixed in https://github.com/yixuan/LBFGSpp/commit/880a8da37b41ae104066411a015fde3efdb94762.