Logo ROOT   6.16/01
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"
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
33namespace ROOT {
34
35 namespace Minuit2 {
36
37
38MnUserParameterState 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
43MnUserParameterState 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
48MnUserParameterState 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
53MnUserParameterState 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
58MnUserParameterState 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
63MnUserParameterState 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
79void 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
88MinimumState 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();
108
109
110 // case gradient is not numeric (could be analytical or from FumiliGradientCalculator)
111
112 if(st.Gradient().IsAnalytical() ) {
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
171L26:
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
190L30:
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 if (n > 0) {
249 MPIProcess mpiprocOffDiagonal(n*(n-1)/2,0);
250 unsigned int startParIndexOffDiagonal = mpiprocOffDiagonal.StartElementIndex();
251 unsigned int endParIndexOffDiagonal = mpiprocOffDiagonal.EndElementIndex();
252
253 unsigned int offsetVect = 0;
254 for (unsigned int in = 0; in<startParIndexOffDiagonal; in++)
255 if ((in+offsetVect)%(n-1)==0) offsetVect += (in+offsetVect)/(n-1);
256
257 for (unsigned int in = startParIndexOffDiagonal;
258 in<endParIndexOffDiagonal; in++) {
259
260 int i = (in+offsetVect)/(n-1);
261 if ((in+offsetVect)%(n-1)==0) offsetVect += i;
262 int j = (in+offsetVect)%(n-1)+1;
263
264 if ((i+1)==j || in==startParIndexOffDiagonal)
265 x(i) += dirin(i);
266
267 x(j) += dirin(j);
268
269 double fs1 = mfcn(x);
270 double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
271 vhmat(i,j) = elem;
272
273 x(j) -= dirin(j);
274
275 if (j%(n-1)==0 || in==endParIndexOffDiagonal-1)
276 x(i) -= dirin(i);
277
278 }
279
280 mpiprocOffDiagonal.SyncSymMatrixOffDiagonal(vhmat);
281 }
282
283 //verify if matrix pos-def (still 2nd derivative)
284
285#ifdef DEBUG
286 std::cout << "Original error matrix " << vhmat << std::endl;
287#endif
288
289 MinimumError tmpErr = MnPosDef()(MinimumError(vhmat,1.), prec);
290
291#ifdef DEBUG
292 std::cout << "Original error matrix " << vhmat << std::endl;
293#endif
294
295 vhmat = tmpErr.InvHessian();
296
297#ifdef DEBUG
298 std::cout << "PosDef error matrix " << vhmat << std::endl;
299#endif
300
301
302 int ifail = Invert(vhmat);
303 if(ifail != 0) {
304
305#ifdef WARNINGMSG
306 MN_INFO_MSG("MnHesse: matrix inversion fails!");
307 MN_INFO_MSG("MnHesse fails and will return diagonal matrix.");
308#endif
309
310 MnAlgebraicSymMatrix tmpsym(vhmat.Nrow());
311 for(unsigned int j = 0; j < n; j++) {
312 double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
313 tmpsym(j,j) = tmp < prec.Eps2() ? 1. : tmp;
314 }
315
316 return MinimumState(st.Parameters(), MinimumError(tmpsym, MinimumError::MnInvertFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
317 }
318
319 FunctionGradient gr(grd, g2, gst);
321
322 // if matrix is made pos def returns anyway edm
323 if(tmpErr.IsMadePosDef()) {
325 double edm = estim.Estimate(gr, err);
326#ifdef WARNINGMSG
327 MN_INFO_MSG("MnHesse: matrix was forced pos. def. ");
328#endif
329 return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
330 }
331
332 //calculate edm for good errors
333 MinimumError err(vhmat, 0.);
334 double edm = estim.Estimate(gr, err);
335
336#ifdef DEBUG
337 std::cout << "\nHesse is ACCURATE. New state from MnHesse " << std::endl;
338 std::cout << "Gradient " << grd << std::endl;
339 std::cout << "Second Deriv " << g2 << std::endl;
340 std::cout << "Gradient step " << gst << std::endl;
341 std::cout << "Error " << vhmat << std::endl;
342 std::cout << "edm " << edm << std::endl;
343#endif
344
345
346 return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
347}
348
349/*
350 MinimumError MnHesse::Hessian(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo) const {
351
352 const MnMachinePrecision& prec = trafo.Precision();
353 // make sure starting at the right place
354 double amin = mfcn(st.Vec());
355 // if(fabs(amin - st.Fval()) > prec.Eps2()) std::cout<<"function Value differs from amin by "<<amin - st.Fval()<<std::endl;
356
357 double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
358
359 // diagonal Elements first
360
361 unsigned int n = st.Parameters().Vec().size();
362 MnAlgebraicSymMatrix vhmat(n);
363 MnAlgebraicVector g2 = st.Gradient().G2();
364 MnAlgebraicVector gst = st.Gradient().Gstep();
365 MnAlgebraicVector grd = st.Gradient().Grad();
366 MnAlgebraicVector dirin = st.Gradient().Gstep();
367 MnAlgebraicVector yy(n);
368 MnAlgebraicVector x = st.Parameters().Vec();
369
370 for(unsigned int i = 0; i < n; i++) {
371
372 double xtf = x(i);
373 double dmin = 8.*prec.Eps2()*fabs(xtf);
374 double d = fabs(gst(i));
375 if(d < dmin) d = dmin;
376 for(int icyc = 0; icyc < Ncycles(); icyc++) {
377 double sag = 0.;
378 double fs1 = 0.;
379 double fs2 = 0.;
380 for(int multpy = 0; multpy < 5; multpy++) {
381 x(i) = xtf + d;
382 fs1 = mfcn(x);
383 x(i) = xtf - d;
384 fs2 = mfcn(x);
385 x(i) = xtf;
386 sag = 0.5*(fs1+fs2-2.*amin);
387 if(sag > prec.Eps2()) break;
388 if(trafo.Parameter(i).HasLimits()) {
389 if(d > 0.5) {
390 std::cout<<"second derivative zero for Parameter "<<i<<std::endl;
391 std::cout<<"return diagonal matrix "<<std::endl;
392 for(unsigned int j = 0; j < n; j++) {
393 vhmat(j,j) = (g2(j) < prec.Eps2() ? 1. : 1./g2(j));
394 return MinimumError(vhmat, 1., false);
395 }
396 }
397 d *= 10.;
398 if(d > 0.5) d = 0.51;
399 continue;
400 }
401 d *= 10.;
402 }
403 if(sag < prec.Eps2()) {
404 std::cout<<"MnHesse: internal loop exhausted, return diagonal matrix."<<std::endl;
405 for(unsigned int i = 0; i < n; i++)
406 vhmat(i,i) = (g2(i) < prec.Eps2() ? 1. : 1./g2(i));
407 return MinimumError(vhmat, 1., false);
408 }
409 double g2bfor = g2(i);
410 g2(i) = 2.*sag/(d*d);
411 grd(i) = (fs1-fs2)/(2.*d);
412 gst(i) = d;
413 dirin(i) = d;
414 yy(i) = fs1;
415 double dlast = d;
416 d = sqrt(2.*aimsag/fabs(g2(i)));
417 if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
418 if(d < dmin) d = dmin;
419
420 // see if converged
421 if(fabs((d-dlast)/d) < Tolerstp()) break;
422 if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
423 d = std::min(d, 10.*dlast);
424 d = std::max(d, 0.1*dlast);
425 }
426 vhmat(i,i) = g2(i);
427 }
428
429 //off-diagonal Elements
430 for(unsigned int i = 0; i < n; i++) {
431 x(i) += dirin(i);
432 for(unsigned int j = i+1; j < n; j++) {
433 x(j) += dirin(j);
434 double fs1 = mfcn(x);
435 double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
436 vhmat(i,j) = elem;
437 x(j) -= dirin(j);
438 }
439 x(i) -= dirin(i);
440 }
441
442 return MinimumError(vhmat, 0.);
443 }
444 */
445
446 } // namespace Minuit2
447
448} // namespace ROOT
#define MN_INFO_VAL2(loc, x)
Definition: MnPrint.h:130
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define d(i)
Definition: RSha256.hxx:102
double sqrt(double)
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:47
virtual double Up() const =0
Error definition of the function.
const MnAlgebraicVector & Gstep() const
const MnAlgebraicVector & Grad() const
const MnAlgebraicVector & G2() const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
void Add(const MinimumState &state)
const MnUserParameterState & UserState() const
const MinimumState & State() const
HessianGradientCalculator: class to calculate Gradient for Hessian.
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int Nrow() const
Definition: LASymMatrix.h:239
unsigned int size() const
Definition: LAVector.h:198
unsigned int StartElementIndex() const
Definition: MPIProcess.h:57
bool SyncSymMatrixOffDiagonal(ROOT::Minuit2::MnAlgebraicSymMatrix &mnmatrix)
Definition: MPIProcess.cxx:201
unsigned int EndElementIndex() const
Definition: MPIProcess.h:61
MinimumError keeps the inv.
Definition: MinimumError.h:26
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
double Up() const
Definition: MnFcn.cxx:35
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
double Tolerstp() const
Definition: MnHesse.h:87
unsigned int Ncycles() const
forward interface of MnStrategy
Definition: MnHesse.h:86
MnUserParameterState operator()(const FCNBase &, const std::vector< double > &, const std::vector< double > &, unsigned int maxcalls=0) const
low-level API
Definition: MnHesse.cxx:38
MnStrategy fStrategy
Definition: MnHesse.h:92
double TolerG2() const
Definition: MnHesse.h:88
determines the relative floating point arithmetic precision.
double Eps2() const
eps2 returns 2*sqrt(eps)
Force the covariance matrix to be positive defined by adding extra terms in the diagonal.
Definition: MnPosDef.h:26
unsigned int Strategy() const
Definition: MnStrategy.h:39
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...
Wrapper used by Minuit of FCN interface containing a reference to the transformation object.
Definition: MnUserFcn.h:26
class which holds the external user and/or internal Minuit representation of the parameters and error...
const std::vector< double > & IntParameters() const
const MnUserTransformation & Trafo() const
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
class dealing with the transformation between user specified parameters (external) and internal param...
unsigned int ExtOfInt(unsigned int internal) const
const char * Name(unsigned int) const
const MnMachinePrecision & Precision() const
forwarded interface
const MinuitParameter & Parameter(unsigned int) const
class performing the numerical gradient calculation
double Estimate(const FunctionGradient &, const MinimumError &) const
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TGraphErrors * gr
Definition: legend1.C:25
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
LASymMatrix MnAlgebraicSymMatrix
Definition: MnMatrix.h:41
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21