Logo ROOT   6.10/09
Reference Guide
VariableMetricBuilder.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 
12 #include "Minuit2/MinimumState.h"
13 #include "Minuit2/MinimumError.h"
16 #include "Minuit2/MnLineSearch.h"
17 #include "Minuit2/MinimumSeed.h"
18 #include "Minuit2/MnFcn.h"
20 #include "Minuit2/MnPosDef.h"
22 #include "Minuit2/LaSum.h"
23 #include "Minuit2/LaProd.h"
24 #include "Minuit2/MnStrategy.h"
25 #include "Minuit2/MnHesse.h"
26 
27 //#define DEBUG
28 #include "Minuit2/MnPrint.h"
29 
30 // #if defined(DEBUG) || defined(WARNINGMSG)
31 // #endif
32 
33 
34 
35 namespace ROOT {
36 
37  namespace Minuit2 {
38 
39 
40 double inner_product(const LAVector&, const LAVector&);
41 
42 void VariableMetricBuilder::AddResult( std::vector<MinimumState>& result, const MinimumState & state) const {
43  // // if (!store) store = StorageLevel();
44  // // store |= (result.size() == 0);
45  // if (store)
46  result.push_back(state);
47  // else {
48  // result.back() = state;
49  // }
50  if (TraceIter() ) TraceIteration(result.size()-1, result.back() );
51  else {
52  if (PrintLevel() > 1) {
53  MnPrint::PrintState(std::cout, result.back(), "VariableMetric: Iteration # ",result.size()-1);
54  }
55  }
56 }
57 
58 
59 FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientCalculator& gc, const MinimumSeed& seed, const MnStrategy& strategy, unsigned int maxfcn, double edmval) const {
60  // top level function to find minimum from a given initial seed
61  // iterate on a minimum search in case of first attempt is not successful
62 
63  // to be consistent with F77 Minuit
64  // in Minuit2 edm is correct and is ~ a factor of 2 smaller than F77Minuit
65  // There are also a check for convergence if (edm < 0.1 edmval for exiting the loop)
66  // LM: change factor to 2E-3 to be consistent with F77Minuit
67  edmval *= 0.002;
68 
69  int printLevel = PrintLevel();
70 
71 
72 #ifdef DEBUG
73  std::cout<<"VariableMetricBuilder convergence when edm < "<<edmval<<std::endl;
74 #endif
75 
76  if(seed.Parameters().Vec().size() == 0) {
77  return FunctionMinimum(seed, fcn.Up());
78  }
79 
80 
81  // double edm = Estimator().Estimate(seed.Gradient(), seed.Error());
82  double edm = seed.State().Edm();
83 
84  FunctionMinimum min(seed, fcn.Up() );
85 
86  if(edm < 0.) {
87 #ifdef WARNINGMSG
88  MN_INFO_MSG("VariableMetricBuilder: initial matrix not pos.def.");
89 #endif
90  //assert(!seed.Error().IsPosDef());
91  return min;
92  }
93 
94  std::vector<MinimumState> result;
95  if (StorageLevel() > 0)
96  result.reserve(10);
97  else
98  result.reserve(2);
99 
100 
101  // do actual iterations
102  if (printLevel >1) {
103  std::cout << "VariableMetric: start iterating until Edm is < " << edmval << std::endl;
104  MnPrint::PrintState(std::cout, seed.State(), "VariableMetric: Initial state ");
105  }
106 
107  AddResult( result, seed.State());
108 
109 
110  // try first with a maxfxn = 80% of maxfcn
111  int maxfcn_eff = maxfcn;
112  int ipass = 0;
113  bool iterate = false;
114 
115  do {
116 
117  iterate = false;
118 
119 #ifdef DEBUG
120  std::cout << "start iterating... " << std::endl;
121  if (ipass > 0) std::cout << "continue iterating... " << std::endl;
122 #endif
123 
124  min = Minimum(fcn, gc, seed, result, maxfcn_eff, edmval);
125 
126  // if max function call reached exits
127  if ( min.HasReachedCallLimit() ) {
128 #ifdef WARNINGMSG
129  MN_INFO_MSG("VariableMetricBuilder: FunctionMinimum is invalid, reached the function call limit");
130 #endif
131  return min;
132  }
133 
134  // second time check for validity of function Minimum
135  if (ipass > 0) {
136  if(!min.IsValid()) {
137 #ifdef WARNINGMSG
138  MN_INFO_MSG("VariableMetricBuilder: FunctionMinimum is invalid after second try");
139 #endif
140  return min;
141  }
142  }
143 
144  // resulting edm of minimization
145  edm = result.back().Edm();
146  // need to re-coorect for Dcovar ?
147 
148  if( (strategy.Strategy() == 2) ||
149  (strategy.Strategy() == 1 && min.Error().Dcovar() > 0.05) ) {
150 
151 #ifdef DEBUG
152  std::cout<<"MnMigrad will verify convergence and Error matrix. "<< std::endl;
153  std::cout<<"dcov is = "<< min.Error().Dcovar() << std::endl;
154 #endif
155 
156  MinimumState st = MnHesse(strategy)(fcn, min.State(), min.Seed().Trafo(),maxfcn);
157 
158  if (printLevel > 1) {
159  MnPrint::PrintState(std::cout, st, "VariableMetric: After Hessian ");
160  }
161  AddResult( result, st);
162 
163  // check new edm
164  edm = st.Edm();
165 #ifdef DEBUG
166  std::cout << "edm after Hesse calculation " << edm << " requested " << edmval << std::endl;
167 #endif
168  if (edm > edmval) {
169  // be careful with machine precision and avoid too small edm
170  double machineLimit = fabs(seed.Precision().Eps2()*result.back().Fval());
171  if (edm >= machineLimit) {
172  iterate = true;
173 #ifdef WARNINGMSG
174  MN_INFO_MSG("VariableMetricBuilder: Tolerance is not sufficient, continue the minimization");
175  MN_INFO_VAL2("Current Edm is",edm);
176  MN_INFO_VAL2("Required Edm is",edmval);
177 #endif
178  }
179  else {
180 #ifdef WARNINGMSG
181  MN_INFO_MSG("VariableMetricBuilder: Stop the minimization - reached machine accuracy limit");
182  MN_INFO_VAL2("Edm is smaller than machine accuracy",machineLimit);
183  MN_INFO_VAL2("Current Edm is",edm);
184  MN_INFO_VAL2("Required Edm is",edmval);
185 #endif
186  }
187 
188  }
189  }
190 
191 
192  // end loop on iterations
193  // ? need a maximum here (or max of function calls is enough ? )
194  // continnue iteration (re-calculate function Minimum if edm IS NOT sufficient)
195  // no need to check that hesse calculation is done (if isnot done edm is OK anyway)
196  // count the pass to exit second time when function Minimum is invalid
197  // increase by 20% maxfcn for doing some more tests
198  if (ipass == 0) maxfcn_eff = int(maxfcn*1.3);
199  ipass++;
200  } while ( iterate );
201 
202 
203  // Add latest state (Hessian calculation)
204  // and check edm (add a factor of 10 in tolerance )
205  if (edm > 10*edmval) {
206  min.Add( result.back(), FunctionMinimum::MnAboveMaxEdm() );
207 #ifdef WARNINGMSG
208  MN_INFO_VAL2("VariableMetricBuilder: INVALID function minimum - edm is above tolerance,",edm);
209  MN_INFO_VAL2("VariableMetricBuilder: Required tolerance is 10 x edmval ",edmval);
210 #endif
211  }
212  else {
213  // check if minimum has edm above max before
214 #ifdef WARNINGMSG
215  if ( min.IsAboveMaxEdm() ) {
216  MN_INFO_MSG("VariableMetricBuilder: Edm has been re-computed after Hesse");
217  MN_INFO_VAL2("new value is now smaller than the required tolerance,",edm);
218  }
219 #endif
220  min.Add( result.back() );
221  }
222 
223 #ifdef DEBUG
224  std::cout << "Obtained function minimum " << min << std::endl;
225 #endif
226 
227  return min;
228 }
229 
230 FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientCalculator& gc, const MinimumSeed& seed, std::vector<MinimumState>& result, unsigned int maxfcn, double edmval) const {
231  // function performing the minimum searches using the Variable Metric algorithm (MIGRAD)
232  // perform first a line search in the - Vg direction and then update using the Davidon formula (Davidon Error updator)
233  // stop when edm reached is less than required (edmval)
234 
235  // after the modification when I iterate on this functions, so it can be called many times,
236  // the seed is used here only to get precision and construct the returned FunctionMinimum object
237 
238 
239 
240  const MnMachinePrecision& prec = seed.Precision();
241 
242 
243  // result.push_back(MinimumState(seed.Parameters(), seed.Error(), seed.Gradient(), edm, fcn.NumOfCalls()));
244  const MinimumState & initialState = result.back();
245 
246  double edm = initialState.Edm();
247 
248 
249 #ifdef DEBUG
250  std::cout << "\n\nDEBUG Variable Metric Builder \nInitial State: "
251  << " Parameter " << initialState.Vec()
252  << " Gradient " << initialState.Gradient().Vec()
253  << " Inv Hessian " << initialState.Error().InvHessian()
254  << " edm = " << initialState.Edm() << std::endl;
255 #endif
256 
257 
258 
259  // iterate until edm is small enough or max # of iterations reached
260  edm *= (1. + 3.*initialState.Error().Dcovar());
261  MnLineSearch lsearch;
262  MnAlgebraicVector step(initialState.Gradient().Vec().size());
263  // keep also prevStep
264  MnAlgebraicVector prevStep(initialState.Gradient().Vec().size());
265 
266  MinimumState s0 = result.back();
267  assert(s0.IsValid() );
268 
269  do {
270 
271  //MinimumState s0 = result.back();
272 
273  step = -1.*s0.Error().InvHessian()*s0.Gradient().Vec();
274 
275 #ifdef DEBUG
276  std::cout << "\n\n---> Iteration - " << result.size()
277  << "\nFval = " << s0.Fval() << " numOfCall = " << fcn.NumOfCalls()
278  << "\nInternal Parameter values " << s0.Vec()
279  << " Newton step " << step << std::endl;
280 #endif
281 
282  // check if derivatives are not zero
283  if ( inner_product(s0.Gradient().Vec(),s0.Gradient().Vec() ) <= 0 ) {
284 #ifdef DEBUG
285  std::cout << "VariableMetricBuilder: all derivatives are zero - return current status" << std::endl;
286 #endif
287  break;
288  }
289 
290 
291  double gdel = inner_product(step, s0.Gradient().Grad());
292 
293 #ifdef DEBUG
294  std::cout << " gdel = " << gdel << std::endl;
295 #endif
296 
297 
298  if(gdel > 0.) {
299 #ifdef WARNINGMSG
300  MN_INFO_MSG("VariableMetricBuilder: matrix not pos.def, gdel > 0");
301  MN_INFO_VAL(gdel);
302 #endif
303  MnPosDef psdf;
304  MinimumState s0new = psdf(s0, prec);
305  step = -1.*s0new.Error().InvHessian()*s0new.Gradient().Vec();
306  // #ifdef DEBUG
307  // std::cout << "After MnPosdef - Error " << s0.Error().InvHessian() << " Gradient " << s0.Gradient().Vec() << " step " << step << std::endl;
308  // #endif
309  gdel = inner_product(step, s0new.Gradient().Grad());
310 #ifdef WARNINGMSG
311  MN_INFO_VAL(gdel);
312 #endif
313  if(gdel > 0.) {
314  AddResult(result, s0new);
315 
316  return FunctionMinimum(seed, result, fcn.Up());
317  }
318  }
319 
320  MnParabolaPoint pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
321 
322  // <= needed for case 0 <= 0
323  if(fabs(pp.Y() - s0.Fval()) <= fabs(s0.Fval())*prec.Eps() ) {
324 #ifdef WARNINGMSG
325  MN_INFO_MSG("VariableMetricBuilder: no improvement in line search");
326 #endif
327  // no improvement exit (is it really needed LM ? in vers. 1.22 tried alternative )
328  // add new state when only fcn changes
329  if (result.size() <= 1 )
330  AddResult(result, MinimumState(s0.Parameters(), s0.Error(), s0.Gradient(), s0.Edm(), fcn.NumOfCalls()));
331  else
332  // no need to re-store the state
333  AddResult(result, MinimumState(pp.Y(), s0.Edm(), fcn.NumOfCalls()));
334 
335  break;
336 
337 
338  }
339 
340 #ifdef DEBUG
341  std::cout << "Result after line search : \nx = " << pp.X()
342  << "\nOld Fval = " << s0.Fval()
343  << "\nNew Fval = " << pp.Y()
344  << "\nNFcalls = " << fcn.NumOfCalls() << std::endl;
345 #endif
346 
347  MinimumParameters p(s0.Vec() + pp.X()*step, pp.Y());
348 
349 
350  FunctionGradient g = gc(p, s0.Gradient());
351 
352 
353  edm = Estimator().Estimate(g, s0.Error());
354 
355 
356  if(edm < 0.) {
357 #ifdef WARNINGMSG
358  MN_INFO_MSG("VariableMetricBuilder: matrix not pos.def. : edm is < 0. Make pos def...");
359 #endif
360  MnPosDef psdf;
361  MinimumState s0new = psdf(s0, prec);
362  edm = Estimator().Estimate(g, s0new.Error());
363  if(edm < 0.) {
364 #ifdef WARNINGMSG
365  MN_INFO_MSG("VariableMetricBuilder: matrix still not pos.def. : exit iterations ");
366 #endif
367  AddResult(result, s0new);
368 
369  return FunctionMinimum(seed, result, fcn.Up());
370  }
371  }
372  MinimumError e = ErrorUpdator().Update(s0, p, g);
373 
374 #ifdef DEBUG
375  std::cout << "Updated new point: \n "
376  << " Parameter " << p.Vec()
377  << " Gradient " << g.Vec()
378  << " InvHessian " << e.Matrix()
379  << " Hessian " << e.Hessian()
380  << " edm = " << edm << std::endl << std::endl;
381 #endif
382 
383  // update the state
384  s0 = MinimumState(p, e, g, edm, fcn.NumOfCalls());
385  if (StorageLevel() || result.size() <= 1)
386  AddResult(result, s0);
387  else
388  // use a reduced state for not-final iterations
389  AddResult(result, MinimumState(p.Fval(), edm, fcn.NumOfCalls()));
390 
391  // correct edm
392  edm *= (1. + 3.*e.Dcovar());
393 
394 #ifdef DEBUG
395  std::cout << "Dcovar = " << e.Dcovar() << "\tCorrected edm = " << edm << std::endl;
396 #endif
397 
398 
399 
400  } while(edm > edmval && fcn.NumOfCalls() < maxfcn); // end of iteration loop
401 
402  // save last result in case of no complete final states
403  if ( ! result.back().IsValid() )
404  result.back() = s0;
405 
406 
407  if(fcn.NumOfCalls() >= maxfcn) {
408 #ifdef WARNINGMSG
409  MN_INFO_MSG("VariableMetricBuilder: call limit exceeded.");
410 #endif
411  return FunctionMinimum(seed, result, fcn.Up(), FunctionMinimum::MnReachedCallLimit());
412  }
413 
414  if(edm > edmval) {
415  if(edm < fabs(prec.Eps2()*result.back().Fval())) {
416 #ifdef WARNINGMSG
417  MN_INFO_MSG("VariableMetricBuilder: machine accuracy limits further improvement.");
418 #endif
419  return FunctionMinimum(seed, result, fcn.Up());
420  } else if(edm < 10*edmval) {
421  return FunctionMinimum(seed, result, fcn.Up());
422  } else {
423 #ifdef WARNINGMSG
424  MN_INFO_MSG("VariableMetricBuilder: iterations finish without convergence.");
425  MN_INFO_VAL2("VariableMetricBuilder",edm);
426  MN_INFO_VAL2(" requested",edmval);
427 #endif
428  return FunctionMinimum(seed, result, fcn.Up(), FunctionMinimum::MnAboveMaxEdm());
429  }
430  }
431  // std::cout<<"result.back().Error().Dcovar()= "<<result.back().Error().Dcovar()<<std::endl;
432 
433 #ifdef DEBUG
434  std::cout << "Exiting successfully Variable Metric Builder \n"
435  << "NFCalls = " << fcn.NumOfCalls()
436  << "\nFval = " << result.back().Fval()
437  << "\nedm = " << edm << " requested = " << edmval << std::endl;
438 #endif
439 
440  return FunctionMinimum(seed, result, fcn.Up());
441 }
442 
443  } // namespace Minuit2
444 
445 } // namespace ROOT
#define MN_INFO_VAL2(loc, x)
Definition: MnPrint.h:130
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
double Estimate(const FunctionGradient &, const MinimumError &) const
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Force the covariance matrix to be positive defined by adding extra terms in the diagonal.
Definition: MnPosDef.h:26
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
unsigned int Strategy() const
Definition: MnStrategy.h:39
const MnAlgebraicVector & Vec() const
double inner_product(const LAVector &, const LAVector &)
double Up() const
Definition: MnFcn.cxx:35
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:51
determines the relative floating point arithmetic precision.
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
const MinimumState & State() const
Definition: MinimumSeed.h:46
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:47
virtual FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
A point of a parabola.
const DavidonErrorUpdator & ErrorUpdator() const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition: MnHesse.h:40
static void PrintState(std::ostream &os, const MinimumState &state, const char *msg, int iter=-1)
Definition: MnPrint.cxx:58
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
MnAlgebraicSymMatrix Matrix() const
Definition: MinimumError.h:58
#define MN_INFO_VAL(x)
Definition: MnPrint.h:116
unsigned int size() const
Definition: LAVector.h:198
void AddResult(std::vector< MinimumState > &result, const MinimumState &state) const
const MinimumError & Error() const
Definition: MinimumState.h:62
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const
const MnAlgebraicVector & Grad() const
double Eps2() const
eps2 returns 2*sqrt(eps)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
MinimumError keeps the inv.
Definition: MinimumError.h:26
void TraceIteration(int iter, const MinimumState &state) const
double result[121]
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const int strategy
Definition: testNdimFit.cxx:46
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
Implements a 1-dimensional minimization along a given direction (i.e.
Definition: MnLineSearch.h:47
const VariableMetricEDMEstimator & Estimator() const
interface class for gradient calculators
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60