Logo ROOT  
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#include "Minuit2/MnPrint.h"
22#include "Minuit2/MPIProcess.h"
23
24namespace ROOT {
25
26namespace Minuit2 {
27
28MnUserParameterState MnHesse::operator()(const FCNBase &fcn, const std::vector<double> &par,
29 const std::vector<double> &err, unsigned int maxcalls) const
30{
31 // interface from vector of params and errors
32 return (*this)(fcn, MnUserParameterState(par, err), maxcalls);
33}
34
35MnUserParameterState MnHesse::operator()(const FCNBase &fcn, const std::vector<double> &par, unsigned int nrow,
36 const std::vector<double> &cov, unsigned int maxcalls) const
37{
38 // interface from vector of params and covariance
39 return (*this)(fcn, MnUserParameterState(par, cov, nrow), maxcalls);
40}
41
42MnUserParameterState MnHesse::operator()(const FCNBase &fcn, const std::vector<double> &par,
43 const MnUserCovariance &cov, unsigned int maxcalls) const
44{
45 // interface from vector of params and covariance
46 return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
47}
48
49MnUserParameterState MnHesse::operator()(const FCNBase &fcn, const MnUserParameters &par, unsigned int maxcalls) const
50{
51 // interface from MnUserParameters
52 return (*this)(fcn, MnUserParameterState(par), maxcalls);
53}
54
56 unsigned int maxcalls) const
57{
58 // interface from MnUserParameters and MnUserCovariance
59 return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
60}
61
63MnHesse::operator()(const FCNBase &fcn, const MnUserParameterState &state, unsigned int maxcalls) const
64{
65 // interface from MnUserParameterState
66 // create a new Minimum state and use that interface
67 unsigned int n = state.VariableParameters();
68 MnUserFcn mfcn(fcn, state.Trafo(), state.NFcn());
70 for (unsigned int i = 0; i < n; i++)
71 x(i) = state.IntParameters()[i];
72 double amin = mfcn(x);
74 MinimumParameters par(x, amin);
75 FunctionGradient gra = gc(par);
76 MinimumState tmp =
77 (*this)(mfcn, MinimumState(par, MinimumError(MnAlgebraicSymMatrix(n), 1.), gra, state.Edm(), state.NFcn()),
78 state.Trafo(), maxcalls);
79
80 return MnUserParameterState(tmp, fcn.Up(), state.Trafo());
81}
82
83void MnHesse::operator()(const FCNBase &fcn, FunctionMinimum &min, unsigned int maxcalls) const
84{
85 // interface from FunctionMinimum to be used after minimization
86 // use last state from the minimization without the need to re-create a new state
87 // do not reset function calls and keep updating them
88 MnUserFcn mfcn(fcn, min.UserState().Trafo(), min.NFcn());
89 MinimumState st = (*this)(mfcn, min.State(), min.UserState().Trafo(), maxcalls);
90 min.Add(st);
91}
92
94 unsigned int maxcalls) const
95{
96 // internal interface from MinimumState and MnUserTransformation
97 // Function who does the real Hessian calculations
98 MnPrint print("MnHesse");
99
100 const MnMachinePrecision &prec = trafo.Precision();
101 // make sure starting at the right place
102 double amin = mfcn(st.Vec());
103 double aimsag = std::sqrt(prec.Eps2()) * (std::fabs(amin) + mfcn.Up());
104
105 // diagonal Elements first
106
107 unsigned int n = st.Parameters().Vec().size();
108 if (maxcalls == 0)
109 maxcalls = 200 + 100 * n + 5 * n * n;
110
111 MnAlgebraicSymMatrix vhmat(n);
112 MnAlgebraicVector g2 = st.Gradient().G2();
113 MnAlgebraicVector gst = st.Gradient().Gstep();
114 MnAlgebraicVector grd = st.Gradient().Grad();
115 MnAlgebraicVector dirin = st.Gradient().Gstep();
117
118 // case gradient is not numeric (could be analytical or from FumiliGradientCalculator)
119
120 if (st.Gradient().IsAnalytical()) {
122 FunctionGradient tmp = igc(st.Parameters());
123 gst = tmp.Gstep();
124 dirin = tmp.Gstep();
125 g2 = tmp.G2();
126 }
127
129
130 print.Debug("Gradient is", st.Gradient().IsAnalytical() ? "analytical" : "numerical", "\n point:", x,
131 "\n fcn :", amin, "\n grad :", grd, "\n step :", gst, "\n g2 :", g2);
132
133 for (unsigned int i = 0; i < n; i++) {
134
135 double xtf = x(i);
136 double dmin = 8. * prec.Eps2() * (std::fabs(xtf) + prec.Eps2());
137 double d = std::fabs(gst(i));
138 if (d < dmin)
139 d = dmin;
140
141 print.Debug("Derivative parameter", i, "d =", d, "dmin =", dmin);
142
143 for (unsigned int icyc = 0; icyc < Ncycles(); icyc++) {
144 double sag = 0.;
145 double fs1 = 0.;
146 double fs2 = 0.;
147 for (unsigned int multpy = 0; multpy < 5; multpy++) {
148 x(i) = xtf + d;
149 fs1 = mfcn(x);
150 x(i) = xtf - d;
151 fs2 = mfcn(x);
152 x(i) = xtf;
153 sag = 0.5 * (fs1 + fs2 - 2. * amin);
154
155 print.Debug("cycle", icyc, "mul", multpy, "\tsag =", sag, "d =", d);
156
157 // Now as F77 Minuit - check that sag is not zero
158 if (sag != 0)
159 goto L30; // break
160 if (trafo.Parameter(i).HasLimits()) {
161 if (d > 0.5)
162 goto L26;
163 d *= 10.;
164 if (d > 0.5)
165 d = 0.51;
166 continue;
167 }
168 d *= 10.;
169 }
170
171 L26:
172 // get parameter name for i
173 // (need separate scope for avoiding compl error when declaring name)
174 print.Warn("2nd derivative zero for parameter", trafo.Name(trafo.ExtOfInt(i)),
175 "; MnHesse fails and will return diagonal matrix");
176
177 for (unsigned int j = 0; j < n; j++) {
178 double tmp = g2(j) < prec.Eps2() ? 1. : 1. / g2(j);
179 vhmat(j, j) = tmp < prec.Eps2() ? 1. : tmp;
180 }
181
183 mfcn.NumOfCalls());
184
185 L30:
186 double g2bfor = g2(i);
187 g2(i) = 2. * sag / (d * d);
188 grd(i) = (fs1 - fs2) / (2. * d);
189 gst(i) = d;
190 dirin(i) = d;
191 yy(i) = fs1;
192 double dlast = d;
193 d = std::sqrt(2. * aimsag / std::fabs(g2(i)));
194 if (trafo.Parameter(i).HasLimits())
195 d = std::min(0.5, d);
196 if (d < dmin)
197 d = dmin;
198
199 print.Debug("g1 =", grd(i), "g2 =", g2(i), "step =", gst(i), "d =", d, "diffd =", std::fabs(d - dlast) / d,
200 "diffg2 =", std::fabs(g2(i) - g2bfor) / g2(i));
201
202 // see if converged
203 if (std::fabs((d - dlast) / d) < Tolerstp())
204 break;
205 if (std::fabs((g2(i) - g2bfor) / g2(i)) < TolerG2())
206 break;
207 d = std::min(d, 10. * dlast);
208 d = std::max(d, 0.1 * dlast);
209 }
210 vhmat(i, i) = g2(i);
211 if (mfcn.NumOfCalls() > maxcalls) {
212
213 // std::cout<<"maxcalls " << maxcalls << " " << mfcn.NumOfCalls() << " " << st.NFcn() << std::endl;
214 print.Warn("Maximum number of allowed function calls exhausted; will return diagonal matrix");
215
216 for (unsigned int j = 0; j < n; j++) {
217 double tmp = g2(j) < prec.Eps2() ? 1. : 1. / g2(j);
218 vhmat(j, j) = tmp < prec.Eps2() ? 1. : tmp;
219 }
220
222 st.Edm(), mfcn.NumOfCalls());
223 }
224 }
225
226 print.Debug("Second derivatives", g2);
227
228 if (fStrategy.Strategy() > 0) {
229 // refine first derivative
230 HessianGradientCalculator hgc(mfcn, trafo, fStrategy);
231 FunctionGradient gr = hgc(st.Parameters(), FunctionGradient(grd, g2, gst));
232 // update gradient and step values
233 grd = gr.Grad();
234 gst = gr.Gstep();
235 }
236
237 // off-diagonal Elements
238 // initial starting values
239 if (n > 0) {
240 MPIProcess mpiprocOffDiagonal(n * (n - 1) / 2, 0);
241 unsigned int startParIndexOffDiagonal = mpiprocOffDiagonal.StartElementIndex();
242 unsigned int endParIndexOffDiagonal = mpiprocOffDiagonal.EndElementIndex();
243
244 unsigned int offsetVect = 0;
245 for (unsigned int in = 0; in < startParIndexOffDiagonal; in++)
246 if ((in + offsetVect) % (n - 1) == 0)
247 offsetVect += (in + offsetVect) / (n - 1);
248
249 for (unsigned int in = startParIndexOffDiagonal; in < endParIndexOffDiagonal; in++) {
250
251 int i = (in + offsetVect) / (n - 1);
252 if ((in + offsetVect) % (n - 1) == 0)
253 offsetVect += i;
254 int j = (in + offsetVect) % (n - 1) + 1;
255
256 if ((i + 1) == j || in == startParIndexOffDiagonal)
257 x(i) += dirin(i);
258
259 x(j) += dirin(j);
260
261 double fs1 = mfcn(x);
262 double elem = (fs1 + amin - yy(i) - yy(j)) / (dirin(i) * dirin(j));
263 vhmat(i, j) = elem;
264
265 x(j) -= dirin(j);
266
267 if (j % (n - 1) == 0 || in == endParIndexOffDiagonal - 1)
268 x(i) -= dirin(i);
269 }
270
271 mpiprocOffDiagonal.SyncSymMatrixOffDiagonal(vhmat);
272 }
273
274 // verify if matrix pos-def (still 2nd derivative)
275
276 print.Debug("Original error matrix", vhmat);
277
278 MinimumError tmpErr = MnPosDef()(MinimumError(vhmat, 1.), prec);
279 vhmat = tmpErr.InvHessian();
280
281 print.Debug("PosDef error matrix", vhmat);
282
283 int ifail = Invert(vhmat);
284 if (ifail != 0) {
285
286 print.Warn("Matrix inversion fails; will return diagonal matrix");
287
288 MnAlgebraicSymMatrix tmpsym(vhmat.Nrow());
289 for (unsigned int j = 0; j < n; j++) {
290 double tmp = g2(j) < prec.Eps2() ? 1. : 1. / g2(j);
291 tmpsym(j, j) = tmp < prec.Eps2() ? 1. : tmp;
292 }
293
295 mfcn.NumOfCalls());
296 }
297
298 FunctionGradient gr(grd, g2, gst);
300
301 // if matrix is made pos def returns anyway edm
302 if (tmpErr.IsMadePosDef()) {
304 double edm = estim.Estimate(gr, err);
305 return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
306 }
307
308 // calculate edm for good errors
309 MinimumError err(vhmat, 0.);
310 double edm = estim.Estimate(gr, err);
311
312 print.Debug("Hessian is ACCURATE. New state:", "\n First derivative:", grd, "\n Second derivative:", g2,
313 "\n Gradient step:", gst, "\n Covariance matrix:", vhmat, "\n Edm:", edm);
314
315 return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
316}
317
318/*
319 MinimumError MnHesse::Hessian(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo) const {
320
321 const MnMachinePrecision& prec = trafo.Precision();
322 // make sure starting at the right place
323 double amin = mfcn(st.Vec());
324 // if(std::fabs(amin - st.Fval()) > prec.Eps2()) std::cout<<"function Value differs from amin by "<<amin -
325 st.Fval()<<std::endl;
326
327 double aimsag = std::sqrt(prec.Eps2())*(std::fabs(amin)+mfcn.Up());
328
329 // diagonal Elements first
330
331 unsigned int n = st.Parameters().Vec().size();
332 MnAlgebraicSymMatrix vhmat(n);
333 MnAlgebraicVector g2 = st.Gradient().G2();
334 MnAlgebraicVector gst = st.Gradient().Gstep();
335 MnAlgebraicVector grd = st.Gradient().Grad();
336 MnAlgebraicVector dirin = st.Gradient().Gstep();
337 MnAlgebraicVector yy(n);
338 MnAlgebraicVector x = st.Parameters().Vec();
339
340 for(unsigned int i = 0; i < n; i++) {
341
342 double xtf = x(i);
343 double dmin = 8.*prec.Eps2()*std::fabs(xtf);
344 double d = std::fabs(gst(i));
345 if(d < dmin) d = dmin;
346 for(int icyc = 0; icyc < Ncycles(); icyc++) {
347 double sag = 0.;
348 double fs1 = 0.;
349 double fs2 = 0.;
350 for(int multpy = 0; multpy < 5; multpy++) {
351 x(i) = xtf + d;
352 fs1 = mfcn(x);
353 x(i) = xtf - d;
354 fs2 = mfcn(x);
355 x(i) = xtf;
356 sag = 0.5*(fs1+fs2-2.*amin);
357 if(sag > prec.Eps2()) break;
358 if(trafo.Parameter(i).HasLimits()) {
359 if(d > 0.5) {
360 std::cout<<"second derivative zero for Parameter "<<i<<std::endl;
361 std::cout<<"return diagonal matrix "<<std::endl;
362 for(unsigned int j = 0; j < n; j++) {
363 vhmat(j,j) = (g2(j) < prec.Eps2() ? 1. : 1./g2(j));
364 return MinimumError(vhmat, 1., false);
365 }
366 }
367 d *= 10.;
368 if(d > 0.5) d = 0.51;
369 continue;
370 }
371 d *= 10.;
372 }
373 if(sag < prec.Eps2()) {
374 std::cout<<"MnHesse: internal loop exhausted, return diagonal matrix."<<std::endl;
375 for(unsigned int i = 0; i < n; i++)
376 vhmat(i,i) = (g2(i) < prec.Eps2() ? 1. : 1./g2(i));
377 return MinimumError(vhmat, 1., false);
378 }
379 double g2bfor = g2(i);
380 g2(i) = 2.*sag/(d*d);
381 grd(i) = (fs1-fs2)/(2.*d);
382 gst(i) = d;
383 dirin(i) = d;
384 yy(i) = fs1;
385 double dlast = d;
386 d = std::sqrt(2.*aimsag/std::fabs(g2(i)));
387 if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
388 if(d < dmin) d = dmin;
389
390 // see if converged
391 if(std::fabs((d-dlast)/d) < Tolerstp()) break;
392 if(std::fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
393 d = std::min(d, 10.*dlast);
394 d = std::max(d, 0.1*dlast);
395 }
396 vhmat(i,i) = g2(i);
397 }
398
399 //off-diagonal Elements
400 for(unsigned int i = 0; i < n; i++) {
401 x(i) += dirin(i);
402 for(unsigned int j = i+1; j < n; j++) {
403 x(j) += dirin(j);
404 double fs1 = mfcn(x);
405 double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
406 vhmat(i,j) = elem;
407 x(j) -= dirin(j);
408 }
409 x(i) -= dirin(i);
410 }
411
412 return MinimumError(vhmat, 0.);
413 }
414 */
415
416} // namespace Minuit2
417
418} // namespace ROOT
#define d(i)
Definition: RSha256.hxx:102
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:45
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, Status status=MnValid)
add latest minimization state (for example add Hesse result after Migrad)
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:45
unsigned int Nrow() const
Definition: LASymMatrix.h:274
unsigned int size() const
Definition: LAVector.h:225
unsigned int StartElementIndex() const
Definition: MPIProcess.h:55
bool SyncSymMatrixOffDiagonal(ROOT::Minuit2::MnAlgebraicSymMatrix &mnmatrix)
Definition: MPIProcess.cxx:193
unsigned int EndElementIndex() const
Definition: MPIProcess.h:61
MinimumError keeps the inv.
Definition: MinimumError.h:28
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:50
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:27
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:30
double Up() const
Definition: MnFcn.cxx:39
unsigned int NumOfCalls() const
Definition: MnFcn.h:39
double Tolerstp() const
Definition: MnHesse.h:89
unsigned int Ncycles() const
forward interface of MnStrategy
Definition: MnHesse.h:88
MnUserParameterState operator()(const FCNBase &, const std::vector< double > &, const std::vector< double > &, unsigned int maxcalls=0) const
low-level API
Definition: MnHesse.cxx:28
MnStrategy fStrategy
Definition: MnHesse.h:93
double TolerG2() const
Definition: MnHesse.h:90
Sets the relative floating point (double) 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:25
void Debug(const Ts &... args)
Definition: MnPrint.h:138
void Warn(const Ts &... args)
Definition: MnPrint.h:126
unsigned int Strategy() const
Definition: MnStrategy.h:38
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:25
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< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:21
LASymMatrix MnAlgebraicSymMatrix
Definition: MnMatrix.h:27
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.