Finding global minima with ROOT

[Kristian]

I am trying to perform a test fit on two histograms with 2 free variables using ROOT::Math::Factory::CreateMinimizer, but none of the algorithms I’ve tried can find the global minima. Previously I’ve had some luck with GSL simulated annealing, but that also fails. Are there any tips on what I can do to improve the minimization? I don’t mind if it increases the computation time – currently the annealing already takes 30 minutes to complete.

---

ROOT Version: 6.22/07

[Kristian]

I’ve tried using the Genetic algorithm instead, which provides far better results at the cost of roughly doubling the computation time. I am using it to fit a simulation to my data, and I really want to increase the number of simulated events, which would drastically increase the computation time. Is there any simple way to modify this genetic fitting algorithm to use multiple CPU cores?

[jblomer]

Hi Kristian,

Can you share a sample of your fitting problem (macro and, ideally, data)? That might help us suggesting an approach.

Cheers,
Jakob

[Kristian]

The entire script is around 700 lines long, so I’ve just pasted one of the evaluation methods. I have two others, but they do essentially the same thing.

As for my data, I think it may be easier to just describe what I am actually doing. I start by binning my data into a 2D histogram (Dalitz plot). The simulation is a mix of two nuclear states, where a weight must be calculated individually for each event depending on k, the ratio of the two states, and delta, which describes interference effects between them. I precalculate the Dalitz coordinates for each event, and then only calculate the weights for each call of my evaluation function. I then apply the weights, bin the simulation data, and compare the two Dalitz plots through a ratio likelihood.

As long as I use the genetic algorithm, it works really well, although it takes quite a while to finish. As I already mentioned, I really want to increase the number of simulated events, which would make the fitting time prohibitively long.

    double eval_type1(const double *params) const {
        const double k = params[0];
        const double delta = params[1]*2*M_PI;

        // some of the minimizers do not respect our limits, which results in nan values
        if (k > 1) { 
            return 1e9;
        }

        // calculate the weight
        vector<double> w = vector<double>(sim1x.size());
        for (int i = 0; i < w.size(); i++) {
            w[i] = calc_weight(f1[i], wU1[i], k, delta);
        }

        // generate the dalitz slices
        TH2D* hsim = dalitz_slice(sim1x, sim1y, bins, w);
        double scale = calc_scale(hsim);
        hsim->Scale(scale);

        // calculate chi
        double chi = maximum_likelihood(hsim);
        hsim->Delete();
        return chi;
    }