Using Minuit2 to minimize a Mathematica function

Question

FindMinimum is great, but sometimes it is difficult to get an idea of the uncertainties and intermediate evaluation diagnostics. How to define a function in Mathematica (9.0) and use its evaluations in another minimization routine, like Minuit2? I am specially interested in using mathlink. For example, Minuit2 spurts an "estimated distance to minimum" and a covariance matrix between the parameters.

Answer 1

This is my try in creating a minimal example with only one variable. There are a couple problems to solve:

1. Minuit2 is in C++, so we need a way to integrate it with mathematica through mathlink. This is a bit cumbersome, but in principle straightforward. This first example uses a mathematica function with minimum user input, that is, it does not uses a user-defined interpolation or anything like that. It would be awesome to define your functions from the front end and then somehow "connect" to this current kernel from the minuit routine.

2. A crucial requisite is not initializing a kernel each time the function is evaluated.

3. The compilation and linking, libraries and all that: I will give everything here.

A few words about Minuit2 (documentation and download): this minimization routine minimizes a function which is inside some classes with fixed names: "FCN" is the class where a parentheses overloading defines the function evaluation. "FCN" is a derived class of "FCNBase".

So, the idea is to start mathlink in the "FCN" constructor, overload the parentheses with the usual mathlink calls (MLPutFunction, MLGetReal etc...) , and then close the link in the "FCN" destructor. In his way we keep most of the Mathematica jargon outside the main (not that it is really needed, but I was fond of this idea).

So, first the header of FCN ("cabeceraFCN.hh")

// cabeceraFCN.hh
#ifndef MN_CABEFCN_H_ 
#define MN_CABEFCN_H_  

#include "Minuit2/FCNBase.h"  
#include <vector> 
#include "mathlink.h"

namespace ROOT {
  namespace Minuit2 {
   class FCN : public FCNBase { 

   private:
    double fErrorDef;
    MLEnvironment m_env;
    MLINK m_lp;
    int m_ac;
    char* m_av[];  

   public:
    FCN(int ac, char* av[]){
     for(int i=0;i<ac;i++){
      m_av[i]=av[i];
     }
     m_ac=ac;
     m_lp=MLOpen(m_ac,m_av);
     m_env = MLInitialize(NULL);
     if(m_env==NULL) exit(EXIT_FAILURE);
     if(m_lp==NULL) exit(EXIT_FAILURE);
    }
  ~FCN(){
    MLClose(m_lp);
    MLDeinitialize(m_env);
   }  
  virtual double Up() const {return fErrorDef;} 
  virtual double operator()(const std::vector<double>&) const;
  void SetErrorDef(double def) {fErrorDef = def;}

    };
  }
}
#endif

And then, the main routine. In this case, we are calculating $\min_c\int_0^1(\exp(-u^2)-c)^2du$, so not to keep it too simplistic.

// file mainML.cpp
#include <iostream>
#include "mathlink.h"

#include "cabeceraFCN.hh"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"

namespace ROOT { 

   namespace Minuit2 {

     double FCN::operator()(const std::vector<double>& parametro) const {

       assert(parametro.size() == 1);
       double x=parametro[0];
       MLINK lp=m_lp;

    {MLPutFunction(lp,"NIntegrate",2);
      {MLPutFunction(lp,"Power",2);
        {MLPutFunction(lp,"Plus",2);
          {MLPutFunction(lp,"Exp",1);
            {MLPutFunction(lp,"Times",2);
              MLPutReal(lp,-1.0);
              {MLPutFunction(lp,"Power",2);
                MLPutSymbol(lp,"u");
                MLPutInteger(lp,2);
              }
            }
          }
          {MLPutFunction(lp,"Times",2);
            MLPutReal(lp,-1.0); 
            MLPutReal(lp,x); 
          }
        }
        MLPutInteger(lp,2);
      }
      {MLPutFunction(lp,"List",3);
        MLPutSymbol(lp,"u");
        MLPutReal(lp,0);
        MLPutReal(lp,1);
      }
    }

    MLEndPacket(lp);
    while(MLNextPacket(lp) !=RETURNPKT) MLNewPacket(lp);

    double chi2;
    MLGetReal(lp,&chi2);
    return chi2;

     }
   }
}

using namespace ROOT::Minuit2;

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

  FCN funFCN(argc,argv);
  //create initial starting values for parameters
  double b0 = 0.;

  // create Minuit parameters with names
  MnUserParameters upar;
  upar.Add("b", b0, 0.1);

  // create MIGRAD minimizer
  MnMigrad migrad(funFCN, upar);

  // Minimize
  FunctionMinimum min = migrad();

  // output
  std::cout<<min<<std::endl;

  return 0;
}

Note that the parameter to minimize is in a vector, so it should hopefully be easy to include more free parameters. For the moment, the only error I get is that when I only create and close the mathlink, without evaluating anything (for example, if above the main function would end in the line "FCN funFCN(argc,argv);"), I get:

 LinkConnect::linkc: -- Message text not found -- (LinkObject[*link name*, 3, 1])

but I do not get anything if I evaluate the function (through the parentheses overloading) at least one time. So this problem is not serious I think.

And now, compilation and linking. I am using the neat trick given in this answer to define local variables for the libraries and include directories. I do not know anything about makefiles. Minuit2 is built in my "$HOME/local" folder.

MATH_INC_C=$(dirname $(readlink -f $(which math)))/../SystemFiles/IncludeFiles/C
MATH_C_ADD=$(dirname $(readlink -f $(which math)))/../SystemFiles/Links/MathLink/DeveloperKit/Linux-x86-64/CompilerAdditions
MATH_L_LIB=$(dirname $(readlink -f $(which math)))/../SystemFiles/Libraries/Linux-x86-64
MIN2_INC=$HOME/local/Minuit2-5.34.14/inc
MIN2_LIB=$HOME/local/Minuit2-5.34.14/src/.libs/libMinuit2.a

g++ -Wall -std=c++11 -m64 -fPIC -O2 -o test   -I $MIN2_INC -I $MATH_INC_C -I $MATH_C_ADD  mainML.cpp $MIN2_LIB -L $MATH_C_ADD -L $MATH_L_LIB -l"ML64i3" -lm -lpthread -lrt -lstdc++ -fopenmp

Before executing, the only issue is that library that mathlink needs at runtime ("libML64i3.so" in my case). You can give the location in "/etc/ld.so.conf.d/" and run "ldconfig", but if you do not have administrative privileges in your computer, then you will have to define the location in "LD_LIBRARY_PATH" (see the previous answer I quote). According to the internet, this is bad, so define it like as a local variable. Then you can run

./test  -linkname 'math -mathlink' -linkmode launch

and it works. In my case, it writes:

Minuit did successfully converge.
# of function calls: 15
minimum function Value: 0.04039772131027
minimum edm: 7.523163849222e-25
minimum internal state vector: LAVector parameters:
     0.7468241328117

minimum internal covariance matrix: LASymMatrix parameters:
              1


# ext. ||   Name    ||   type  ||     Value     ||  Error +/- 

   0   ||         b ||  free   ||     0.7468241328117 ||2.63063186581e-155

At this point, I think I am happy with an unstructured communication with the front end to get the minimum and the covariance matrix.

Any comments or suggestions are deeply appreciated, since I do not know much about C++ or mathlink.

---

EDIT and UPDATE

I said that it should be easy to include more than one free parameter, and of course I jynxed it. I have not been able to make this thing work for more than one free parameter and I cannot find the reason. A minimum effort "not working example " could be just add a second parameter squared to the return of the FCN function above. That is, no second parameter is passed to mathlink and that is enough, in my case (Math 9.0, Linux 64 bit Opensuse 13.1), to make the program crash.

Answer 2

Maybe in the age of Google not necessary, but there is way going through Python to add minuit to Mathematica : Github page of the project.

The reason for using minuit is that it is trusted by some of the physics communities. If you use something else, people might ask questions that are hard to answer.