Logo ROOT   6.18/05
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
42void 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
59FunctionMinimum 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
230FunctionMinimum 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 s0 = psdf(s0, prec);
305 step = -1.*s0.Error().InvHessian()*s0.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, s0.Gradient().Grad());
310#ifdef WARNINGMSG
311 MN_INFO_VAL(gdel);
312#endif
313 if(gdel > 0.) {
314 AddResult(result, s0);
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 s0 = psdf(s0, prec);
362 edm = Estimator().Estimate(g, s0.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, s0);
368
369 return FunctionMinimum(seed, result, fcn.Up());
370 }
371 }
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
#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
determines the relative floating point 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 &)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21