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