Logo ROOT   6.10/09
Reference Guide
MnHesse.cxx
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
10 #include "Minuit2/MnHesse.h"
12 #include "Minuit2/MnUserFcn.h"
13 #include "Minuit2/FCNBase.h"
14 #include "Minuit2/MnPosDef.h"
18 #include "Minuit2/MinimumState.h"
21 
22 //#define DEBUG
23 
24 #if defined(DEBUG) || defined(WARNINGMSG)
25 #include "Minuit2/MnPrint.h"
26 #endif
27 #if defined(DEBUG) && !defined(WARNINGMSG)
28 #define WARNINGMSG
29 #endif
30 
31 #include "Minuit2/MPIProcess.h"
32 
33 namespace ROOT {
34 
35  namespace Minuit2 {
36 
37 
38 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const std::vector<double>& err, unsigned int maxcalls) const {
39  // interface from vector of params and errors
40  return (*this)(fcn, MnUserParameterState(par, err), maxcalls);
41 }
42 
43 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, unsigned int nrow, const std::vector<double>& cov, unsigned int maxcalls) const {
44  // interface from vector of params and covariance
45  return (*this)(fcn, MnUserParameterState(par, cov, nrow), maxcalls);
46 }
47 
48 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const MnUserCovariance& cov, unsigned int maxcalls) const {
49  // interface from vector of params and covariance
50  return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
51 }
52 
53 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameters& par, unsigned int maxcalls) const {
54  // interface from MnUserParameters
55  return (*this)(fcn, MnUserParameterState(par), maxcalls);
56 }
57 
58 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameters& par, const MnUserCovariance& cov, unsigned int maxcalls) const {
59  // interface from MnUserParameters and MnUserCovariance
60  return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
61 }
62 
63 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameterState& state, unsigned int maxcalls) const {
64  // interface from MnUserParameterState
65  // create a new Minimum state and use that interface
66  unsigned int n = state.VariableParameters();
67  MnUserFcn mfcn(fcn, state.Trafo(),state.NFcn());
69  for(unsigned int i = 0; i < n; i++) x(i) = state.IntParameters()[i];
70  double amin = mfcn(x);
72  MinimumParameters par(x, amin);
73  FunctionGradient gra = gc(par);
74  MinimumState tmp = (*this)(mfcn, MinimumState(par, MinimumError(MnAlgebraicSymMatrix(n), 1.), gra, state.Edm(), state.NFcn()), state.Trafo(), maxcalls);
75 
76  return MnUserParameterState(tmp, fcn.Up(), state.Trafo());
77 }
78 
79 void MnHesse::operator()(const FCNBase& fcn, FunctionMinimum& min, unsigned int maxcalls) const {
80  // interface from FunctionMinimum to be used after minimization
81  // use last state from the minimization without the need to re-create a new state
82  // do not reset function calls and keep updating them
83  MnUserFcn mfcn(fcn, min.UserState().Trafo(),min.NFcn());
84  MinimumState st = (*this)( mfcn, min.State(), min.UserState().Trafo(), maxcalls);
85  min.Add(st);
86 }
87 
88 MinimumState MnHesse::operator()(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo, unsigned int maxcalls) const {
89  // internal interface from MinimumState and MnUserTransformation
90  // Function who does the real Hessian calculations
91 
92  const MnMachinePrecision& prec = trafo.Precision();
93  // make sure starting at the right place
94  double amin = mfcn(st.Vec());
95  double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
96 
97  // diagonal Elements first
98 
99  unsigned int n = st.Parameters().Vec().size();
100  if(maxcalls == 0) maxcalls = 200 + 100*n + 5*n*n;
101 
102  MnAlgebraicSymMatrix vhmat(n);
103  MnAlgebraicVector g2 = st.Gradient().G2();
104  MnAlgebraicVector gst = st.Gradient().Gstep();
105  MnAlgebraicVector grd = st.Gradient().Grad();
106  MnAlgebraicVector dirin = st.Gradient().Gstep();
107  MnAlgebraicVector yy(n);
108 
109 
110  // case gradient is not numeric (could be analytical or from FumiliGradientCalculator)
111 
112  if(st.Gradient().IsAnalytical() ) {
113  Numerical2PGradientCalculator igc(mfcn, trafo, fStrategy);
114  FunctionGradient tmp = igc(st.Parameters());
115  gst = tmp.Gstep();
116  dirin = tmp.Gstep();
117  g2 = tmp.G2();
118  }
119 
121 
122 #ifdef DEBUG
123  std::cout << "\nMnHesse " << std::endl;
124  std::cout << " x " << x << std::endl;
125  std::cout << " amin " << amin << " " << st.Fval() << std::endl;
126  std::cout << " grd " << grd << std::endl;
127  std::cout << " gst " << gst << std::endl;
128  std::cout << " g2 " << g2 << std::endl;
129  std::cout << " Gradient is analytical " << st.Gradient().IsAnalytical() << std::endl;
130 #endif
131 
132 
133  for(unsigned int i = 0; i < n; i++) {
134 
135  double xtf = x(i);
136  double dmin = 8.*prec.Eps2()*(fabs(xtf) + prec.Eps2());
137  double d = fabs(gst(i));
138  if(d < dmin) d = dmin;
139 
140 #ifdef DEBUG
141  std::cout << "\nDerivative parameter " << i << " d = " << d << " dmin = " << dmin << std::endl;
142 #endif
143 
144 
145  for(unsigned int icyc = 0; icyc < Ncycles(); icyc++) {
146  double sag = 0.;
147  double fs1 = 0.;
148  double fs2 = 0.;
149  for(unsigned int multpy = 0; multpy < 5; multpy++) {
150  x(i) = xtf + d;
151  fs1 = mfcn(x);
152  x(i) = xtf - d;
153  fs2 = mfcn(x);
154  x(i) = xtf;
155  sag = 0.5*(fs1+fs2-2.*amin);
156 
157 #ifdef DEBUG
158  std::cout << "cycle " << icyc << " mul " << multpy << "\t sag = " << sag << " d = " << d << std::endl;
159 #endif
160  // Now as F77 Minuit - check taht sag is not zero
161  if (sag != 0) goto L30; // break
162  if(trafo.Parameter(i).HasLimits()) {
163  if(d > 0.5) goto L26;
164  d *= 10.;
165  if(d > 0.5) d = 0.51;
166  continue;
167  }
168  d *= 10.;
169  }
170 
171 L26:
172 #ifdef WARNINGMSG
173 
174  // get parameter name for i
175  // (need separate scope for avoiding compl error when declaring name)
176  {
177  const char * name = trafo.Name( trafo.ExtOfInt(i));
178  MN_INFO_VAL2("MnHesse: 2nd derivative zero for Parameter ", name);
179  MN_INFO_MSG("MnHesse fails and will return diagonal matrix ");
180  }
181 #endif
182 
183  for(unsigned int j = 0; j < n; j++) {
184  double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
185  vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
186  }
187 
188  return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
189 
190 L30:
191  double g2bfor = g2(i);
192  g2(i) = 2.*sag/(d*d);
193  grd(i) = (fs1-fs2)/(2.*d);
194  gst(i) = d;
195  dirin(i) = d;
196  yy(i) = fs1;
197  double dlast = d;
198  d = sqrt(2.*aimsag/fabs(g2(i)));
199  if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
200  if(d < dmin) d = dmin;
201 
202 #ifdef DEBUG
203  std::cout << "\t g1 = " << grd(i) << " g2 = " << g2(i) << " step = " << gst(i) << " d = " << d
204  << " diffd = " << fabs(d-dlast)/d << " diffg2 = " << fabs(g2(i)-g2bfor)/g2(i) << std::endl;
205 #endif
206 
207 
208  // see if converged
209  if(fabs((d-dlast)/d) < Tolerstp()) break;
210  if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
211  d = std::min(d, 10.*dlast);
212  d = std::max(d, 0.1*dlast);
213  }
214  vhmat(i,i) = g2(i);
215  if(mfcn.NumOfCalls() > maxcalls) {
216 
217 #ifdef WARNINGMSG
218  //std::cout<<"maxcalls " << maxcalls << " " << mfcn.NumOfCalls() << " " << st.NFcn() << std::endl;
219  MN_INFO_MSG("MnHesse: maximum number of allowed function calls exhausted.");
220  MN_INFO_MSG("MnHesse fails and will return diagonal matrix ");
221 #endif
222 
223  for(unsigned int j = 0; j < n; j++) {
224  double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
225  vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
226  }
227 
228  return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
229  }
230 
231  }
232 
233 #ifdef DEBUG
234  std::cout << "\n Second derivatives " << g2 << std::endl;
235 #endif
236 
237  if(fStrategy.Strategy() > 0) {
238  // refine first derivative
239  HessianGradientCalculator hgc(mfcn, trafo, fStrategy);
240  FunctionGradient gr = hgc(st.Parameters(), FunctionGradient(grd, g2, gst));
241  // update gradient and step values
242  grd = gr.Grad();
243  gst = gr.Gstep();
244  }
245 
246  //off-diagonal Elements
247  // initial starting values
248  MPIProcess mpiprocOffDiagonal(n*(n-1)/2,0);
249  unsigned int startParIndexOffDiagonal = mpiprocOffDiagonal.StartElementIndex();
250  unsigned int endParIndexOffDiagonal = mpiprocOffDiagonal.EndElementIndex();
251 
252  unsigned int offsetVect = 0;
253  for (unsigned int in = 0; in<startParIndexOffDiagonal; in++)
254  if ((in+offsetVect)%(n-1)==0) offsetVect += (in+offsetVect)/(n-1);
255 
256  for (unsigned int in = startParIndexOffDiagonal;
257  in<endParIndexOffDiagonal; in++) {
258 
259  int i = (in+offsetVect)/(n-1);
260  if ((in+offsetVect)%(n-1)==0) offsetVect += i;
261  int j = (in+offsetVect)%(n-1)+1;
262 
263  if ((i+1)==j || in==startParIndexOffDiagonal)
264  x(i) += dirin(i);
265 
266  x(j) += dirin(j);
267 
268  double fs1 = mfcn(x);
269  double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
270  vhmat(i,j) = elem;
271 
272  x(j) -= dirin(j);
273 
274  if (j%(n-1)==0 || in==endParIndexOffDiagonal-1)
275  x(i) -= dirin(i);
276 
277  }
278 
279  mpiprocOffDiagonal.SyncSymMatrixOffDiagonal(vhmat);
280 
281  //verify if matrix pos-def (still 2nd derivative)
282 
283 #ifdef DEBUG
284  std::cout << "Original error matrix " << vhmat << std::endl;
285 #endif
286 
287  MinimumError tmpErr = MnPosDef()(MinimumError(vhmat,1.), prec);
288 
289 #ifdef DEBUG
290  std::cout << "Original error matrix " << vhmat << std::endl;
291 #endif
292 
293  vhmat = tmpErr.InvHessian();
294 
295 #ifdef DEBUG
296  std::cout << "PosDef error matrix " << vhmat << std::endl;
297 #endif
298 
299 
300  int ifail = Invert(vhmat);
301  if(ifail != 0) {
302 
303 #ifdef WARNINGMSG
304  MN_INFO_MSG("MnHesse: matrix inversion fails!");
305  MN_INFO_MSG("MnHesse fails and will return diagonal matrix.");
306 #endif
307 
308  MnAlgebraicSymMatrix tmpsym(vhmat.Nrow());
309  for(unsigned int j = 0; j < n; j++) {
310  double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
311  tmpsym(j,j) = tmp < prec.Eps2() ? 1. : tmp;
312  }
313 
314  return MinimumState(st.Parameters(), MinimumError(tmpsym, MinimumError::MnInvertFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
315  }
316 
317  FunctionGradient gr(grd, g2, gst);
319 
320  // if matrix is made pos def returns anyway edm
321  if(tmpErr.IsMadePosDef()) {
323  double edm = estim.Estimate(gr, err);
324 #ifdef WARNINGMSG
325  MN_INFO_MSG("MnHesse: matrix was forced pos. def. ");
326 #endif
327  return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
328  }
329 
330  //calculate edm for good errors
331  MinimumError err(vhmat, 0.);
332  double edm = estim.Estimate(gr, err);
333 
334 #ifdef DEBUG
335  std::cout << "\nHesse is ACCURATE. New state from MnHesse " << std::endl;
336  std::cout << "Gradient " << grd << std::endl;
337  std::cout << "Second Deriv " << g2 << std::endl;
338  std::cout << "Gradient step " << gst << std::endl;
339  std::cout << "Error " << vhmat << std::endl;
340  std::cout << "edm " << edm << std::endl;
341 #endif
342 
343 
344  return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
345 }
346 
347 /*
348  MinimumError MnHesse::Hessian(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo) const {
349 
350  const MnMachinePrecision& prec = trafo.Precision();
351  // make sure starting at the right place
352  double amin = mfcn(st.Vec());
353  // if(fabs(amin - st.Fval()) > prec.Eps2()) std::cout<<"function Value differs from amin by "<<amin - st.Fval()<<std::endl;
354 
355  double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
356 
357  // diagonal Elements first
358 
359  unsigned int n = st.Parameters().Vec().size();
360  MnAlgebraicSymMatrix vhmat(n);
361  MnAlgebraicVector g2 = st.Gradient().G2();
362  MnAlgebraicVector gst = st.Gradient().Gstep();
363  MnAlgebraicVector grd = st.Gradient().Grad();
364  MnAlgebraicVector dirin = st.Gradient().Gstep();
365  MnAlgebraicVector yy(n);
366  MnAlgebraicVector x = st.Parameters().Vec();
367 
368  for(unsigned int i = 0; i < n; i++) {
369 
370  double xtf = x(i);
371  double dmin = 8.*prec.Eps2()*fabs(xtf);
372  double d = fabs(gst(i));
373  if(d < dmin) d = dmin;
374  for(int icyc = 0; icyc < Ncycles(); icyc++) {
375  double sag = 0.;
376  double fs1 = 0.;
377  double fs2 = 0.;
378  for(int multpy = 0; multpy < 5; multpy++) {
379  x(i) = xtf + d;
380  fs1 = mfcn(x);
381  x(i) = xtf - d;
382  fs2 = mfcn(x);
383  x(i) = xtf;
384  sag = 0.5*(fs1+fs2-2.*amin);
385  if(sag > prec.Eps2()) break;
386  if(trafo.Parameter(i).HasLimits()) {
387  if(d > 0.5) {
388  std::cout<<"second derivative zero for Parameter "<<i<<std::endl;
389  std::cout<<"return diagonal matrix "<<std::endl;
390  for(unsigned int j = 0; j < n; j++) {
391  vhmat(j,j) = (g2(j) < prec.Eps2() ? 1. : 1./g2(j));
392  return MinimumError(vhmat, 1., false);
393  }
394  }
395  d *= 10.;
396  if(d > 0.5) d = 0.51;
397  continue;
398  }
399  d *= 10.;
400  }
401  if(sag < prec.Eps2()) {
402  std::cout<<"MnHesse: internal loop exhausted, return diagonal matrix."<<std::endl;
403  for(unsigned int i = 0; i < n; i++)
404  vhmat(i,i) = (g2(i) < prec.Eps2() ? 1. : 1./g2(i));
405  return MinimumError(vhmat, 1., false);
406  }
407  double g2bfor = g2(i);
408  g2(i) = 2.*sag/(d*d);
409  grd(i) = (fs1-fs2)/(2.*d);
410  gst(i) = d;
411  dirin(i) = d;
412  yy(i) = fs1;
413  double dlast = d;
414  d = sqrt(2.*aimsag/fabs(g2(i)));
415  if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
416  if(d < dmin) d = dmin;
417 
418  // see if converged
419  if(fabs((d-dlast)/d) < Tolerstp()) break;
420  if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
421  d = std::min(d, 10.*dlast);
422  d = std::max(d, 0.1*dlast);
423  }
424  vhmat(i,i) = g2(i);
425  }
426 
427  //off-diagonal Elements
428  for(unsigned int i = 0; i < n; i++) {
429  x(i) += dirin(i);
430  for(unsigned int j = i+1; j < n; j++) {
431  x(j) += dirin(j);
432  double fs1 = mfcn(x);
433  double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
434  vhmat(i,j) = elem;
435  x(j) -= dirin(j);
436  }
437  x(i) -= dirin(i);
438  }
439 
440  return MinimumError(vhmat, 0.);
441  }
442  */
443 
444  } // namespace Minuit2
445 
446 } // namespace ROOT
void Add(const MinimumState &state)
virtual double Up() const =0
Error definition of the function.
#define MN_INFO_VAL2(loc, x)
Definition: MnPrint.h:130
double par[1]
Definition: unuranDistr.cxx:38
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
double Estimate(const FunctionGradient &, const MinimumError &) const
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
Force the covariance matrix to be positive defined by adding extra terms in the diagonal.
Definition: MnPosDef.h:26
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int Strategy() const
Definition: MnStrategy.h:39
double Up() const
Definition: MnFcn.cxx:35
unsigned int ExtOfInt(unsigned int internal) const
unsigned int EndElementIndex() const
Definition: MPIProcess.h:60
const std::vector< double > & IntParameters() const
determines the relative floating point arithmetic precision.
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
LASymMatrix MnAlgebraicSymMatrix
Definition: MnMatrix.h:41
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
double sqrt(double)
MnStrategy fStrategy
Definition: MnHesse.h:92
Double_t x[n]
Definition: legend1.C:17
class performing the numerical gradient calculation
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
const MnMachinePrecision & Precision() const
forwarded interface
unsigned int StartElementIndex() const
Definition: MPIProcess.h:56
const MnAlgebraicVector & G2() const
const MinimumState & State() const
const char * Name(unsigned int) const
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
const MnUserTransformation & Trafo() const
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:47
HessianGradientCalculator: class to calculate Gradient for Hessian.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
unsigned int Nrow() const
Definition: LASymMatrix.h:239
class which holds the external user and/or internal Minuit representation of the parameters and error...
unsigned int size() const
Definition: LAVector.h:198
const MnAlgebraicVector & Gstep() const
TGraphErrors * gr
Definition: legend1.C:25
Wrapper used by Minuit of FCN interface containing a reference to the transformation object...
Definition: MnUserFcn.h:26
class dealing with the transformation between user specified parameters (external) and internal param...
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
double TolerG2() const
Definition: MnHesse.h:88
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
const MnAlgebraicVector & Grad() const
double Eps2() const
eps2 returns 2*sqrt(eps)
const MinuitParameter & Parameter(unsigned int) const
MinimumError keeps the inv.
Definition: MinimumError.h:26
MnUserParameterState operator()(const FCNBase &, const std::vector< double > &, const std::vector< double > &, unsigned int maxcalls=0) const
low-level API
Definition: MnHesse.cxx:38
double Tolerstp() const
Definition: MnHesse.h:87
bool SyncSymMatrixOffDiagonal(ROOT::Minuit2::MnAlgebraicSymMatrix &mnmatrix)
Definition: MPIProcess.cxx:201
unsigned int Ncycles() const
forward interface of MnStrategy
Definition: MnHesse.h:86
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const MnUserParameterState & UserState() const
const Int_t n
Definition: legend1.C:16
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...