Hello,

I'm trying to understand how to use Levenberg Marquardt's algorithm within Eigen.
I already have an Octave working code, I'm now trying to write it in C++.

I don't have e a Jacobian estimation so I want to use the Numerical differentiation module of Eigen

I have found some examples of code but I really don't get how it works. This topic relates the same problem: viewtopic.php?f=74&t=97683&p=206972 (the user can't share the code )
I started from an example code on StackOverflow and I want to fully understand it before continuing.

I came to this code (my questions are in the comments)

Code: Select all

#include <iostream>
#include <Eigen/Dense>

#include <unsupported/Eigen/NonLinearOptimization>
#include <unsupported/Eigen/NumericalDiff>

std::string
printExitStatus (int status);

// Generic functor
// http://stackoverflow.com/questions/356950/c-functors-and-their-uses
template <typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic>
struct Functor
{
    typedef _Scalar Scalar;
    enum
    {
      InputsAtCompileTime = NX, ValuesAtCompileTime = NY
    };
    typedef Eigen::Matrix<Scalar, InputsAtCompileTime, 1> InputType;
    typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, 1> ValueType;
    typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, InputsAtCompileTime> JacobianType;

    int m_inputs, m_values;

    Functor () :
        m_inputs (InputsAtCompileTime),
        m_values (ValuesAtCompileTime)
    {
    }

    Functor (int inputs,
             int values) :
        m_inputs (inputs),
        m_values (values)
    {
    }

    int
    inputs () const  // What is this ? It never gets called
    {
      //std::cout << "m_inputs:" << m_inputs << std::endl;
      return m_inputs;
    }

    int
    values () const  // What is this ?
    {
      //std::cout << "m_values:" << m_values << std::endl;
      return m_values;
    }
};

struct my_functor : Functor<double>
{
    my_functor (void) :
        Functor<double> (3, 3)  // 3 inputs, 3 values (constraints?). How do I "choose" these values ?!
    {
    }

    int
    operator () (const Eigen::VectorXd &x,     // x is the actual estimation
                 Eigen::VectorXd &fvec) const  // a reference to the current vector function, what is fvec ?
    {
      std::cout << "x:\n" << x << "\n" << std::endl;
      std::cout << "fvec:\n" << fvec << " | norm: " << fvec.norm() << "\n" << std::endl;
      //std::cout << "fvec: " << &fvec << "\nx: " << &x << std::endl;

      // Random equations just for fun
      fvec (0) = 10.0 * pow (x (0) + 3.0, 2) + pow (x (1) - 5.0, 2);
      fvec (1) = 0;
      fvec (2) = 0.5 * x (2) - x (0) + 5;

      return (0);  // is this return value only here to state operator() went fine ?
    }
};

int
main (int argc,
      char *argv[])
{
  Eigen::VectorXd wpr (3);
  //wpr << -2.0, 5.6;
  wpr << -3.5, 2.0, -25.0;
  std::cout << "Initial guess of the solution:\n" << wpr << std::endl;  // The answer to the problem is: -3, 5, -16

  my_functor functor;  // Instantiate a Functor structure
  Eigen::NumericalDiff<my_functor> numDiff (functor);  // Instantiate a numerical diff for Jacobian estimation
  Eigen::LevenbergMarquardt<Eigen::NumericalDiff<my_functor>, double> lm (numDiff);  // Instantiate a LM object
  lm.parameters.maxfev = 2000;    // Maximum number of evaluations
  lm.parameters.xtol = 1.0e-10;  // Tolerance

  int ret = lm.minimize (wpr);  // Minimize using LM algorithm; initial guess is "x"
  std::cout << "\nIterations: " << lm.iter << std::endl;
  std::cout << "Exit code: " << ret << " = " << printExitStatus (ret) << std::endl;

  std::cout << "\nx that minimizes the function:\n" << wpr << std::endl;  // x has been minimized
  return (0);
}

std::string
printExitStatus (int status)
{
  switch (status)
  {
    case -2:
      return ("NotStarted");
      break;
    case -1:
      return ("Running");
      break;
    case 0:
      return ("ImproperInputParameters");
      break;
    case 1:
      return ("RelativeReductionTooSmall");
      break;
    case 2:
      return ("RelativeErrorTooSmall");
      break;
    case 3:
      return ("RelativeErrorAndReductionTooSmall");
      break;
    case 4:
      return ("CosinusTooSmall");
      break;
    case 5:
      return ("TooManyFunctionEvaluation");
      break;
    case 6:
      return ("FtolTooSmall");
      break;
    case 7:
      return ("XtolTooSmall");
      break;
    case 8:
      return ("GtolTooSmall");
      break;
    case 9:
      return ("UserAsked");
      break;
    default:
      return ("Unknown exit code");
      break;
  }
  return ("error");
}



Bye