Minuit2 Changing only 1 out of four parameters in ChiSquared minimization

Dear all,

I am using Minuit2 to find the best values of four parameters through chi-squared minimization. For some reason, Minuit2 seems to change only one of these parameters, namely params[3], regardless of what I set the parameters to be initially. I really don’t understand why this is happening, especially because I have used Minuit2 in the past and I have never encountered this issue before. Is this something that happens when Minuit cannot find a minimum? And is there a way around it?

Below is my code:

#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FCNBase.h"

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <cmath>
#include <Eigen/Dense>

using namespace ROOT::Minuit2;
using namespace Eigen;

const double PI = 4*atan(1);
// Define minimizer class (inherits from FCNBase)
class Minimizer : public FCNBase {

private:

  const std::vector< std::vector<double> > &m_data;
  //  const std::vector<double> &m_mass;

public:

  Minimizer(const std::vector< std::vector<double> > &data) : m_data(data) {}

  ~Minimizer() {}
  // functions to find factorial and double factorial of given number - used in the Vqq matrix element
 //Omitted here (but they are in the actual code)

  double diagonalization(double n, double L, double p0, double p1, double p2, double p3){
    // Declare the number of terms in the expansion and the number of iterations for the numerical integrals
    int N1 = 12;

    // Declare the variational parameters - values taken from the hybrid decay paper: https://arxiv.org/abs/2005.10850
    double alphaInit = 0.005;
    // Define the matrices - the dimension can be set to be any of the above indexes since they are generally the same
    MatrixXd HMX(N1, N1);
    MatrixXd NMX(N1, N1);
    MatrixXd KMX(N1, N1);
    MatrixXd VMX(N1, N1);

    //Declare the variables for the minimum eigenvalue
    double minEvalue0 = 0.0; // ground state eigenvalue
    double minEvalue1 = 0.0; // first excited state eigenvalue
    double minEvalue2 = 0.0; // second excited state eigenvalue
    double minEvalue3 = 0.0; // third excited state eigenvalue

    std::vector<double> eValues0;
    std::vector<double> eValues1;
    std::vector<double> eValues2;
    std::vector<double> eValues3;
    std::vector<double> minEvalues;

    // Loop over the variational parameters and the matrix elements (k1=n', l1=m', k2=n, l2=m).
    for (double alpha = alphaInit; alpha <= 1.000; alpha += 0.005){
      for (int k1 = 1; k1 <= N1; k1++){
        for (int k2 = 1; k2 <= N1; k2++){
          // Define the full Hamiltonian in configuration space - works for all values of L
          // Overlap Matrix
          NMX(k1-1,k2-1) = pow(2.0,L + 1)*sqrt(2.0)*pow(k1*k2,0.75 + L/2.0)/(pow(k1+k2,1.5 + L));
          // Mass and kinetic terms
          KMX(k1-1,k2-1) = (pow(2.0,L + 1)*sqrt(2.0)*pow(k1*k2,0.75 + L/2.0)/(pow(k1+k2,1.5 + L))) * (2.0*p0 + (3.0+2.0*L)*(k1*k2)*alpha*alpha/(p0*(k1+k2)));
          // Potential term
          VMX(k1-1,k2-1) = (powerfactorial(L)*pow(k1*k2,0.75 + L/2.0)/(doublefactorial(2*L+1)*sqrt(PI)*alpha*pow(k1+k2,2.0+L)))*(2*p1*alpha*alpha*(k1+k2)-3*(L+1)*p2);
        }
      }
      // Diagonalize and print the eigenvalues
      HMX = KMX - (4.0/3.0)*VMX;

      std::setprecision(8);
      GeneralizedSelfAdjointEigenSolver<MatrixXd> es(HMX,NMX);
      //pushback eigenvalues into vectors
      eValues0.push_back(es.eigenvalues()[0]);
      eValues1.push_back(es.eigenvalues()[1]);
      eValues2.push_back(es.eigenvalues()[2]);
      eValues3.push_back(es.eigenvalues()[3]);

    }
    // Calculate the minimum eigenvalue and the corresponding parameters
    minEvalue0 = eValues0[0];
    minEvalue1 = eValues1[0];
    minEvalue2 = eValues2[0];
    minEvalue3 = eValues3[0];

    //cout << "\n";
    //cout << "Here is the full list of eigenvalues: " << std::endl;

    for(unsigned int i = 0; i < eValues0.size(); i++){
      if (eValues0[i] < minEvalue0){
        minEvalue0 = eValues0[i];
      }
    }
    for(unsigned int i = 0; i < eValues1.size(); i++){
      if (eValues1[i] < minEvalue1){
        minEvalue1 = eValues1[i];
      }
    }

    for(unsigned int i = 0; i < eValues2.size(); i++){
      if (eValues2[i] < minEvalue2){
        minEvalue2 = eValues2[i];
      }
    }

    for(unsigned int i = 0; i < eValues3.size(); i++){
      if (eValues3[i] < minEvalue3){
        minEvalue3 = eValues3[i];
      }
    }
    //push_back all the minimized eigenvalues into a vector
    minEvalues.push_back(minEvalue0);
    minEvalues.push_back(minEvalue1);
    minEvalues.push_back(minEvalue2);
    minEvalues.push_back(minEvalue3);

    double eigenvalue = minEvalues[n-1]-p3;//return the desired eigenvalue from minEvalues
    return eigenvalue;
  }

  // Define function that selects the state for diagonalization(n, L, p0, p1, p2, p3)
  double mass (int state, double p0, double p1, double p2, double p3){
    double n;
    double L;
    double mesonMass = 0.0;

    if (state == 0){//1S
        n = 1.0;
        L = 0.0;
    }
    else if (state == 1){//1P
        n = 1.0;
        L = 1.0;
    }
    else if (state == 2){//2S
        n = 2.0;
        L = 0.0;
    }
    else if (state == 3){//1D
        n = 1.0;
        L = 2.0;
 }
    else if (state == 4){//2P
        n = 2.0;
        L = 1.0;
    }
    else if (state == 5){//3S
        n = 3.0;
        L = 0.0;
    }
    else if (state == 6){//4S
        n = 4.0;
        L = 0.0;
    }
    mesonMass = diagonalization(n, L, p0, p1, p2, p3);
    return mesonMass;
  }

  double operator() (const std::vector<double> &params) const {
    double X2 = 0.0;
    //std::cout << mass(0,params[0], params[1], params[2], params[3]) << std::endl;
    //std::vector<double> mass;
    for (int i=0; i<m_data[0].size(); i++) {
      double n = m_data[0][i];
      double L = m_data[1][i];
      double f = mass(i, params[0], params[1], params[2], params[3]);
      X2 += (m_data[2][i] - f)*(m_data[2][i] - f)/(m_data[2][i]);
      //std::cout << "data[2][i] = " << m_data[2][i] << std::endl;
      //std::cout << "mass[i]" << f << std::endl;
    }
    return X2;
  }
  //set to 1 for chi-sq and 0.5 for log-likelihood
  double Up() const {return 1.0;}
};

int main(){

  std::vector< std::vector<double> > data(3);
  // (1S) n=1, L=0; N=0
  data[0].push_back(1.0);
  data[1].push_back(0.0);
  data[2].push_back(9.4449);
  // (1P) n=2, L=1; N=0
  data[0].push_back(1.0);
  data[1].push_back(1.0);
  data[2].push_back(9.8999);
  // (2S) n=2, L=0; N=1
  data[0].push_back(2.0);
  data[1].push_back(0.0);
  data[2].push_back(10.0173);
  // (1D) n=3, L=2; N=0
  data[0].push_back(1.0);
  data[1].push_back(2.0);
  data[2].push_back(10.1637);
  // (2P) n=3, L=1; N=1
  data[0].push_back(2.0);
  data[1].push_back(1.0);
  data[2].push_back(10.2602);
  // (3S) n=3, L=0; N=2
  data[0].push_back(3.0);
  data[1].push_back(0.0);
  data[2].push_back(10.3552);
  // (4S) n=4, L=0; N=3
  data[0].push_back(4.0);
  data[1].push_back(0.0);
  data[2].push_back(10.5794);

  std::vector<double> params(4,1.0);
  std::vector<double> err(4,0.1);
  params[0] = 5.00; //mc
  params[1] = 0.40;  //alphaS
  params[2] = 0.14; //b
  params[3] = 0.18; //mass shift
  std::cout << "Initial values of parameters: " << std::endl;
  for (int i = 0; i < 4; i++){
    std::cout << "params[" << i << "] = " << params[i] << std::endl;
  }
  std::cout << "Computing the values of the parameters..." << std::endl;

  //minimize
  Minimizer fcn(data);
  
  MnMigrad migrad(fcn, params, err);
  FunctionMinimum minParams = migrad();
 std::cout << "CHISQUARED MINIMIZATION RESULTS:" << std::endl;
  std::cout << "\n";
  std::cout << "***************************************" << std::endl;
  std::cout << " chi^2 = " << minParams.Fval() << std::endl;
  std::cout << " for the following values:\n " << "mass q: " << minParams.UserState().Value(0) << ", alphaS: " << minParams.UserState().Value(1) << ", b: " \
            << minParams.UserState().Value(2) << ", shift: " << minParams.UserState().Value(3) << std::endl;
  std::cout <<"****************************************" << std::endl;
  std::cout << "\n";
  //std::cout << "Did the minimizer find a minimum?" << minParams.isValid() << std::endl;
  std::cout << "migrad() results: " << minParams << std::endl;
  std::cout << "Process completed." << std::endl;
  return 0;
}

The output of this is the following:

CHISQUARED MINIMIZATION RESULTS:

***************************************
 chi^2 = 0.00278467
 for the following values:
 mass q: 5.1997, alphaS: 0.4, b: 0.14, shift: 0.637963
****************************************

migrad() results:
WARNING: Minuit did not converge.

# of function calls: 124
minimum function Value: 0.002784672250699
minimum edm: 6.678767958345
minimum internal state vector: LAVector parameters:
      5.199699062206
     0.3999996341642
     0.1399999440163
     0.6379634861636

minimum internal covariance matrix: LASymMatrix parameters:
  1.4000616e-10  4.9897517e-11  1.0394972e-11 -6.9624951e-06
  4.9897517e-11  1.9757816e-11  3.9417031e-12 -2.6078556e-06
  1.0394972e-11  3.9417031e-12   8.394459e-13 -5.4326561e-07
 -6.9624951e-06 -2.6078556e-06 -5.4326561e-07     0.59767527


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

   0   ||        p0 ||  free   ||      5.199699062206 ||1.18324199021e-05
   1   ||        p1 ||  free   ||     0.3999996341642 ||4.444976434746e-06
   2   ||        p2 ||  free   ||     0.1399999440163 ||9.162128011203e-07
   3   ||        p3 ||  free   ||     0.6379634861636 ||0.773094606389



WARNING: FunctionMinimum is invalid:
         Edm is above max


Process completed.

As you can see, only params[3] is changed in the process, and this is true regardless of what initial values I give to the parameters. They remain unchanged, even if the end result is obviously off from the data. When I ran the code for a different set of data, Minuit was converging rather quickly, giving me the optimized values of parameters.

Any help/suggestion would be appreciated.

Thanks.


Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Linux Mint 17
_Compiler:_gcc


Hi,

What Minuit2 does looks correct to me, giving your function. I see that your Covariance matrix has only the (3,3) element significantly different than zero. Since the step size, s, is equals to s = - C * g
where C is the covariance matrix (inverse Hessian) and g is the function gradient, the result is a step only in par[3].
Why this is happening you have to see how your chi2 function is implemented. If it is chi2 it’s minimum value should be around to the number of degrees of freedom of your problem.

Best regards

Lorenzo

Hi @moneta,

and thank you for your reply. What is strange to me is that when I ran the same exact code for a different set of data (and different values for only par[0] and par[3], keeping par[1] and par[2] the same) Minuit converges and varies all four parameters. Could Minuit have trouble because I am trying to fit 4 parameters with only seven data points?

Thanks again.

Christian

In principle 7 datapoint and 4 parameters are still ok (you have 3 degree of freedom). I will try to plot your chi2 as function of the parameters , for example chi2((par[0],par[3]) and chi2( par[0],par[1])

Lorenzo