Chi2 minimization using multidimensional grid

Dear experts,

I want to fit some datasets to the output of a numerical model, where 5 parameters are left free.
The problem is that this numerical model is very slow, even running it in parallel on a GPU or on multiple CPUs (>= 5 minutes), so it’s not feasible to call it from within a TF1.
In addition, I have multiple datasets to fit, so it’s not efficient to run the model with the same combination of parameters more than once, like it would happen with a TF1.
So I created a 5D grid spanning reasonable ranges in the parameters of interest, run the model for all the grid points (order of 1e6) and for each datasets I computed the chi2 for every grid points and saved them into separate TTrees.
Then I created a class to provide a FCN for ROOT::Fit::Fitter, by linearly interpolating the chi2 values from the 5D grid:

class GridChi2
{
   public:
      GridChi2(TTree *t_chi2, ModelsCollection *mc) :
         t_chi2(t_chi2), mc(mc), npars(5)
      {
         t_chi2->SetBranchAddress("chi2", &chi2);
      }

      // 5D linear interpolation of chi2 from the grid saved in t_chi2
      Double_t operator()(const Double_t *par)
      {
         UShort_t ipars[npars];
         Double_t s[npars];
         mc->GetLinearInterpolatingFactors(par[0], par[1], par[2], par[3], par[4], ipars, s);

         Double_t chi2 = 0.;
         for (UShort_t n = 0; n < 1<<npars; ++n)
         {
            UShort_t shift[npars] = { 0 };
            Double_t k = 1.;
            for (UShort_t ipar = 0; ipar < npars; ++ipar)
            {
               shift[ipar] = (n >> ipar) & 1;
               if (shift[ipar]) k*= s[ipar];
               else k *= 1 - s[ipar];
               ipars[ipar] += shift[ipar];
            }

            Long64_t imodel = mc->iModel(ipars);
            t_chi2->GetEntry(imodel);
            chi2 += k*model_chi2;

            for (UShort_t ipar = 0; ipar < npars; ++ipar)
            {
               ipars[ipar] -= shift[ipar];
            }
         }

         return chi2;
      }

   private:
      TTree *t_chi2;
      ModelsCollection *mc;

      Float_t model_chi2;

      const UShort_t npars;
};

ModelsCollection is a class which contains information about the grid and the models output: its code is not important here.
I tested the multidimensional linear interpolating code and it does what it’s supposed to do.
Then I fit the data in this way:

ROOT::Fit::Fitter fitter;
fitter.Config().MinimizerOptions().SetMaxIterations(5e4);
fitter.Config().MinimizerOptions().SetMaxFunctionCalls(1e5);
fitter.Config().MinimizerOptions().SetTolerance(1);
fitter.Config().MinimizerOptions().SetPrintLevel(3);
fitter.Config().MinimizerOptions().SetStrategy(2);
fitter.Config().ParamsSettings().resize(npars);
for (UShort_t ipar = 0; ipar < npars; ++ipar)
{
   fitter.Config().ParSettings(ipar).Set(par_name[ipar], par_best[ipar], grid_step[ipar]/2., max(0., par_min[ipar] - grid_step[ipar]), par_max[ipar] + grid_step[ipar]);
}

GridChi2 grid_chi2(t_chi2, &mc);
fitter.FitFCN(npars, grid_chi2, 0, dof, true);

where par_best is an array with the values of the parameters for which the chi2 in the grid is minimum.
I tried all combinations of minimizer and algorithm, but none of them is working.
Migrad, Minimizer and MigradImproved fail with status = 4 (in reality because the Hessian matrix is not positively defined), and so does Minos.
Simplex and Scan converge, but they do not move much far from the starting point (even if the starting point is very far from the minimum of the grid…) and the errors are too small to make sense.
All others algorithms (Seek, GSL*, Genetic) either do not converge or seems to run indefinitely.

For now, I’m doing a sort of manual fitting: for each parameter, I marginalize the chi2 vs the parameter, then compute the parabola passing through the lowest chi2 and its neighbor points (one left, one right of the minimum). This way I can extract for each parameter a “best fit” value and the classic 1 sigma error by imposing chi2_err = chi2_min + 1.
The problem is that sometimes:

  1. for some parameter the marginalized chi2 curve is not parabolic at all, so the minimum is ill-defined.
  2. each parameter has a different minimum chi2, because the parabola is different for each parameter: so what is the actual chi2 of the global fit?
  3. the minimum of the parabola can become negative if the marginalized chi2 curve is very narrow: this of course does not make any sense and will give a smaller error than the actual one.

Essentially, I’m reimplementing a very simple version of Minos, with no covariance matrix, nor 2D contours.
That’s why I wanted to use ROOT::Fit::Fitter…

Any idea if my GridChi2 class and the way I use it make sense?
Or how to solve the problems in my manual fitting method?

Thanks.
Claudio

Hi,

I suspect the problem is with the linear interpolation. Doing a linear interpolation you will have discontinuities in the derivative matrix and this could be a problem for Migrad.
Maybe a simple gradient based method, as those in GSL could be better.

Are you 100% sure also that the linear interpolation is correct ?

Lorenzo

Are you 100% sure also that the linear interpolation is correct ?

I tested my interpolation code for 1 parameter (keeping fixed all parameters except one) against TGraph::Eval and for 2 parameters by eye, then it becomes hard to visualize it for more parameters :slight_smile:
Still, I’m assuming that the generalization from 2 to higher dimensions is correct.

Maybe a simple gradient based method, as those in GSL could be better.

All GSLMultiMin methods fails with status = 27 (the iteration do not make any more progress towards the solution).

If I use the Genetic minimizer with PrintLevel=3, it looks like it explores a wider range of parameters with respect to Migrad and then it starts to converge to values closer to my best guess, but it’s painfully slow: 10k iterations in one hour and was still running even if the fitness function values was 3e-6 (starting from 200).
What’s its stopping criteria? Maybe I could ask for less iterations or a higher tolerance?

Anyway, I was curious, how is Minos coping with the problems I had doing the chi2 marginalization separately for each parameter?
In the end, I would be fine with using my “manual fitting”, as long as I’m sure that the error estimation is reasonable (which I’m not, hence the need to use ROOT::Fit::Fitter) and if I was able to build 2D contours for each pair of parameters.

Hi,

I would double check the result with some independent code. For example CERNLIB provided a function FINT for linear interpolation in 5 dimension, but it i in Fortran. Maybe on the Web you can find translations to C/C++.
But it is true, some peculiarities of the data you have can make the convergence very difficult.

In case of Genetic, you need to ask for a much larger tolerance. Eventually you could try to use a gradient based method afterwards.

I guess what you mean with marminalizatuion is keeping all parameter fixed and minimise with respect to the single one varying. This is not really a marginalisation.
Minos works differently, because it computes the profile likelihood/chi2 by scanning one parameter and for each value of that parameter, it will not consider the others fixed, but it will find their best (minimum) value.

Lorenzo

I’ll check my code with the FINT function, thanks for the suggestion.

I managed to get Genetic to work, but it’s not computing any error…

Maybe marginalization is not the correct term, but what I’m doing is to take one parameter and for each value of this parameter I find the minimum chi2, so all the other parameters are not kept fixed, and then I find a curve of minimum chi2 vs parameter value. Then I repeat this for all parameters.
It looks to me that it’s similar to what Minos is doing, right?

Hi,

For each value of the parameter you consider you need to find the chi2 minimum (i.e. you minimise with respect to all other parameters). This is the procedure used by Minos.

Lorenzo

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.