GradientDescent.cxx
Go to the documentation of this file.
2 
3 #include "Minuit2/FCNGradientBase.h"
4 #include "Minuit2/FunctionMinimum.h"
5 #include "Minuit2/MnUserParameters.h"
6 
7 #include <iostream>
8 
9 // Globals set to the default settings
10 ROOT::Minuit2::MnUserParameterState gState;
11 ROOT::Minuit2::MnStrategy gStrat;
12 
13 namespace ana
14 {
15  //----------------------------------------------------------------------
17  GradientDescent(const ROOT::Minuit2::FCNGradientBase& func,
18  const ROOT::Minuit2::MnUserParameters& pars)
19  : ROOT::Minuit2::MnApplication(func, gState, gStrat),
20  fFunc(func), fPars(pars)
21  {
22  }
23 
24  //----------------------------------------------------------------------
25  ROOT::Minuit2::FunctionMinimum GradientDescent::
26  operator()(unsigned int /*maxfcn*/, double /*tolerance*/)
27  {
28  // Initialize at the input parameters
29  std::vector<double> pt = fPars.Params();
30  const unsigned int N = pt.size();
31 
32  // Use the user errors as a crude way to set the initial step size
33  double step = Magnitude(fPars.Errors());
34 
35  // Usually we reuse the chisq from the last iteration, so need to
36  // initialize it once out here.
37  double chi = fFunc(pt);
38  int ncalls = 1;
39 
40  while(true){
41  // std::cout << "New direction, chisq = " << chi << std::endl;
42 
43  // Evaluate gradient to figure out a direction to step in
44  std::vector<double> grad = fFunc.Gradient(pt);
45  ++ncalls;
46  // Make gradient into a unit vector but keep the magnitude around
47  const double gradmag = Magnitude(grad);
48  MakeUnit(grad);
49 
50  // Take step in that direction
51  std::vector<double> trialpt = pt;
52  for(unsigned int i = 0; i < N; ++i) trialpt[i] -= grad[i]*step;
53 
54  // And evaluate the chisq there
55  const double trialchi = fFunc(trialpt);
56  ++ncalls;
57 
58  // std::cout << "step = " << step << " -> chisq " << trialchi << std::endl;
59 
60  // Estimate the second derivative from the two points and one
61  // gradient
62  const double d2 = (trialchi+gradmag*step-chi)/(step*step);
63  // We expect the function to be convex
64  if(d2 > 0){
65  // std::cout << " d2 = " << d2;
66  // If so, we can compute a better step size, one that would place us
67  // right at the minimum if the function only had quadratic terms.
68  step = gradmag/(2*d2);
69  // std::cout << " new step = " << step << std::endl;
70  }
71 
72  while(true){
73  // Keep trying steps until we find one that reduces the chisq
74  std::vector<double> newpt = pt;
75  for(unsigned int i = 0; i < N; ++i) newpt[i] -= grad[i]*step;
76 
77  const double newchi = fFunc(newpt);
78  ++ncalls;
79  // std::cout << " step = " << step << " -> chisq = " << newchi << std::endl;
80 
81  if(newchi > chi){
82  // If the chisq went up instead, try again with smaller step
83  step /= 2;
84  if(step < 1e-8){
85  // std::cout << "Step too small!" << std::endl;
86  // Maybe it's just because we're already in the minimum. If we
87  // can't improve even with very tiny steps, call it a day.
88  return Package(pt, chi, ncalls);
89  }
90  continue;
91  }
92  else{
93  if(chi-newchi < 1e-5){
94  // std::cout << "Small chisq step" << std::endl;
95  // If the improvement we did get is really tiny we're also likely
96  // already at the minimum
97  return Package(newpt, newchi, ncalls);
98  }
99 
100  // In all other cases (ie we took a reasonable step and found a
101  // reasonably better chisq) we want to update our state, take another
102  // look at the gradient vector to figure out which direction to go
103  // next, and preserve our step size, which is some kind of good
104  // estimate of the scale of the function.
105  pt = newpt;
106  chi = newchi;
107  break;
108  }
109  } // end line search
110  } // end while (searching for better pt)
111  }
112 
113  //----------------------------------------------------------------------
114  ROOT::Minuit2::FunctionMinimum GradientDescent::
115  Package(const std::vector<double>& pt, double chi, int ncalls) const
116  {
117  // In practice we only use the chisq and the parameter values, so don't be
118  // too careful about setting this all up.
119  const unsigned int N = pt.size();
120  ROOT::Minuit2::MnAlgebraicVector vec(N);
121  for(unsigned int i = 0; i < N; ++i) vec(i) = pt[i];
122  ROOT::Minuit2::MinimumParameters params(vec, chi);
123  ROOT::Minuit2::MinimumState state(params, 0, ncalls);
124  ROOT::Minuit2::MnUserTransformation trans(pt, std::vector<double>(N));
125  ROOT::Minuit2::MinimumSeed seed(state, trans);
126  return ROOT::Minuit2::FunctionMinimum(seed, {state}, chi);
127  }
128 
129  //----------------------------------------------------------------------
130  double GradientDescent::Magnitude(const std::vector<double>& xs) const
131  {
132  double ret = 0;
133  for(double x: xs) ret += x*x;
134  return sqrt(ret);
135  }
136 
137  //----------------------------------------------------------------------
138  void GradientDescent::MakeUnit(std::vector<double>& xs) const
139  {
140  const double mag = Magnitude(xs);
141  for(double& x: xs) x /= mag;
142  }
143 }
const ROOT::Minuit2::MnUserParameters & fPars
double Magnitude(const std::vector< double > &xs) const
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ROOT::Minuit2::MnUserParameterState gState
T sqrt(T number)
Definition: d0nt_math.hpp:156
GradientDescent(const ROOT::Minuit2::FCNGradientBase &func, const ROOT::Minuit2::MnUserParameters &pars)
ROOT::Minuit2::MnStrategy gStrat
static void grad(vari *vi)
Definition: grad.hpp:30
ROOT::Minuit2::FunctionMinimum Package(const std::vector< double > &pt, double chi, int ncalls) const
unsigned int seed
Definition: runWimpSim.h:102
double func(double x, double y)
Eigen::VectorXd vec
const ROOT::Minuit2::FCNGradientBase & fFunc
std::string pars("Th23Dmsq32")
virtual ROOT::Minuit2::FunctionMinimum operator()(unsigned int maxfcn, double tolerance) override
void MakeUnit(std::vector< double > &xs) const
Float_t e
Definition: plot.C:35