Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
FumiliBuilder.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
18#include "Minuit2/MinimumSeed.h"
19#include "Minuit2/MnFcn.h"
21#include "Minuit2/MnPosDef.h"
22#include "Minuit2/MnStrategy.h"
23#include "Minuit2/MnHesse.h"
24#include "Minuit2/MnPrint.h"
25
26namespace ROOT {
27
28namespace Minuit2 {
29
31 const MnStrategy &strategy, unsigned int maxfcn, double edmval) const
32{
33 // top level function to find minimum from a given initial seed
34 // iterate on a minimum search in case of first attempt is not successful
35
36 MnPrint print("FumiliBuilder", PrintLevel());
37
38 edmval *= 0.0001;
39 // edmval *= 0.1; // use small factor for Fumili
40
41 print.Debug("Convergence when edm <", edmval);
42
43 if (seed.Parameters().Vec().size() == 0) {
44 print.Warn("No variable parameters are defined! - Return current function value ");
45 return FunctionMinimum(seed, fcn.Up());
46 }
47
48 // estimate initial edm value
49 double edm = Estimator().Estimate(seed.Gradient(), seed.Error());
50 print.Debug("initial edm is ", edm);
51 //double edm = seed.State().Edm();
52
53 FunctionMinimum min(seed, fcn.Up());
54
55 if (edm < 0.) {
56 print.Error("Initial matrix not positive defined, edm = ",edm,"\nExit minimization ");
57 return min;
58 }
59
60 std::vector<MinimumState> result;
61 // result.reserve(1);
62 result.reserve(8);
63
64 result.push_back(seed.State());
65
66 print.Info("Start iterating until Edm is <", edmval, '\n', "Initial state", MnPrint::Oneline(seed.State()));
67
68 if (TraceIter())
69 TraceIteration(result.size() - 1, result.back());
70
71 // do actual iterations
72
73 // try first with a maxfxn = 50% of maxfcn
74 // Fumili in principle needs much less iterations
75 int maxfcn_eff = int(0.5 * maxfcn);
76 int ipass = 0;
77 double edmprev = 1;
78
79 do {
80
81 min = Minimum(fcn, gc, seed, result, maxfcn_eff, edmval);
82
83 // second time check for validity of function Minimum
84 if (ipass > 0) {
85 if (!min.IsValid()) {
86
87 print.Warn("FunctionMinimum is invalid");
88
89 return min;
90 }
91 }
92
93 // resulting edm of minimization
94 edm = result.back().Edm();
95
96 print.Debug("Approximate Edm", edm, "npass", ipass);
97
98 // call always Hesse (error matrix from Fumili is never accurate since is approximate)
99
100 if (strategy.Strategy() > 0) {
101 print.Debug("FumiliBuilder will verify convergence and Error matrix; "
102 "dcov",
103 min.Error().Dcovar());
104
105 // // recalculate the gradient using the numerical gradient calculator
106 // Numerical2PGradientCalculator ngc(fcn, min.Seed().Trafo(), strategy);
107 // FunctionGradient ng = ngc( min.State().Parameters() );
108 // MinimumState tmp( min.State().Parameters(), min.State().Error(), ng, min.State().Edm(),
109 // min.State().NFcn() );
110
111 MinimumState st = MnHesse(strategy)(fcn, min.State(), min.Seed().Trafo(), maxfcn);
112 result.push_back(st);
113
114 print.Info("After Hessian");
115 if (TraceIter())
116 TraceIteration(result.size() - 1, result.back());
117
118 // check edm
119 edm = st.Edm();
120
121 print.Debug("Edm", edm, "State", st);
122 }
123
124 // break the loop if edm is NOT getting smaller
125 if (ipass > 0 && edm >= edmprev) {
126 print.Warn("Stop iterations, no improvements after Hesse; current Edm", edm, "previous value", edmprev);
127 break;
128 }
129 if (edm > edmval) {
130 print.Debug("Tolerance not sufficient, continue minimization; Edm", edm, "Requested", edmval);
131 } else {
132 // Case when edm < edmval after Heasse but min is flagged with a bad edm:
133 // make then a new Function minimum since now edm is ok
134 if (min.IsAboveMaxEdm()) {
135 min = FunctionMinimum(min.Seed(), min.States(), min.Up());
136 break;
137 }
138 }
139
140 // end loop on iterations
141 // ? need a maximum here (or max of function calls is enough ? )
142 // continnue iteration (re-calculate function Minimum if edm IS NOT sufficient)
143 // no need to check that hesse calculation is done (if isnot done edm is OK anyway)
144 // count the pass to exit second time when function Minimum is invalid
145 // increase by 20% maxfcn for doing some more tests
146 if (ipass == 0)
148 ipass++;
149 edmprev = edm;
150
151 } while (edm > edmval);
152
153 // Add latest state (Hessian calculation)
154 min.Add(result.back());
155
156 return min;
157}
158
160 std::vector<MinimumState> &result, unsigned int maxfcn, double edmval) const
161{
162
163 // function performing the minimum searches using the FUMILI algorithm
164 // after the modification when I iterate on this functions, so it can be called many times,
165 // the seed is used here only to get precision and construct the returned FunctionMinimum object
166
167 /*
168 Three options were possible:
169
170 1) create two parallel and completely separate hierarchies, in which case
171 the FumiliMinimizer would NOT inherit from ModularFunctionMinimizer,
172 FumiliBuilder would not inherit from MinimumBuilder etc
173
174 2) Use the inheritance (base classes of ModularFunctionMinimizer,
175 MinimumBuilder etc), but recreate the member functions Minimize() and
176 Minimum() respectively (naming them for example minimize2() and
177 minimum2()) so that they can take FumiliFCNBase as Parameter instead FCNBase
178 (otherwise one wouldn't be able to call the Fumili-specific methods).
179
180 3) Cast in the daughter classes derived from ModularFunctionMinimizer,
181 MinimumBuilder.
182
183 The first two would mean to duplicate all the functionality already existent,
184 which is a very bad practice and Error-prone. The third one is the most
185 elegant and effective solution, where the only constraint is that the user
186 must know that they have to pass a subclass of FumiliFCNBase to the FumiliMinimizer
187 and not just a subclass of FCNBase.
188 BTW, the first two solutions would have meant to recreate also a parallel
189 structure for MnFcn...
190 **/
191 // const FumiliFCNBase* tmpfcn = dynamic_cast<const FumiliFCNBase*>(&(fcn.Fcn()));
192
193 const MnMachinePrecision &prec = seed.Precision();
194
195 const MinimumState &initialState = result.back();
196
197 double edm = initialState.Edm();
198
199 MnPrint print("FumiliBuilder");
200
201 print.Info("Iterating FumiliBuilder", maxfcn, edmval);
202
203 print.Debug("Initial State:", "\n Parameter", initialState.Vec(), "\n Gradient", initialState.Gradient().Vec(),
204 "\n Inv Hessian", initialState.Error().InvHessian(), "\n edm", initialState.Edm(), "\n maxfcn",
205 maxfcn, "\n tolerance", edmval);
206
207 // iterate until edm is small enough or max # of iterations reached
208 edm *= (1. + 3. * initialState.Error().Dcovar());
209 MnAlgebraicVector step(initialState.Gradient().Vec().size());
210
211 // initial lambda Value
212
213
214 const bool doLineSearch = (fMethodType == kLineSearch);
216 const bool scaleTR = (fMethodType == kTrustRegionScaled);
217 // trust region parameter
218 // use as initial value 0.3 * || x0 ||
219 auto & x0 = initialState.Vec();
220 double normX0 = std::sqrt(inner_product(x0,x0));
221 double delta = 0.3 * std::max(1.0 , normX0);
222 const double eta = 0.1;
223 // use same values as used in GSL implementation
224 const double tr_factor_up = 3;
225 const double tr_factor_down = 0.5;
226 bool acceptNewPoint = true;
227
228 if (doLineSearch)
229 print.Info("Using Fumili with a line search algorithm");
230 else if (doTrustRegion && scaleTR)
231 print.Info("Using Fumili with a scaled trust region algorithm with factors up/down",tr_factor_up,tr_factor_down);
232 else
233 print.Info("Using Fumili with a trust region algorithm with factors up/down",tr_factor_up,tr_factor_down);
234
235
236 double lambda = (doLineSearch) ? 0.001 : 0;
237
239
240 do {
241
242 // const MinimumState& s0 = result.back();
243 MinimumState s0 = result.back();
244
245 step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
246
247 print.Debug("Iteration -", result.size(), "\n Fval", s0.Fval(), "numOfCall", fcn.NumOfCalls(),
248 "\n Internal Parameter values", s0.Vec(), "\n Newton step", step);
249
250 double gdel = inner_product(step, s0.Gradient().Grad());
251 //not sure if this is needed ??
252 if (gdel > 0.) {
253 print.Warn("Matrix not pos.def, gdel =", gdel, " > 0");
254
256 s0 = psdf(s0, prec);
257 step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
258 gdel = inner_product(step, s0.Gradient().Grad());
259
260 print.Warn("After correction, gdel =", gdel);
261
262 if (gdel > 0.) {
263 result.push_back(s0);
264 return FunctionMinimum(seed, result, fcn.Up());
265 }
266 }
267
268 // take a full step
269
270 //evaluate function only if doing a line search
271 double fval2 = (doLineSearch) ? fcnCaller(s0.Vec() + step) : 0;
272 MinimumParameters p(s0.Vec() + step, fval2);
273
274 // check that taking the full step does not deteriorate minimum
275 // in that case do a line search
276 if (doLineSearch && p.Fval() >= s0.Fval()) {
277 print.Debug("Do a line search", fcn.NumOfCalls());
279 auto pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
280
281 if (std::fabs(pp.Y() - s0.Fval()) < prec.Eps()) {
282 // std::cout<<"FumiliBuilder: no improvement"<<std::endl;
283 break; // no improvement
284 }
285 p = MinimumParameters(s0.Vec() + pp.X() * step, pp.Y());
286 print.Debug("New point after Line Search :", "\n FVAL ", p.Fval(), "\n Parameter", p.Vec());
287 }
288 // use as scaling matrix diagonal of Hessian
289 auto & H = s0.Error().Hessian();
290 unsigned int n = (scaleTR) ? H.Nrow() : 0;
291
296 // set the scaling matrix to the sqrt(diagoinal hessian)
297 for (unsigned int i = 0; i < n; i++){
298 double d = std::sqrt(H(i,i));
299 D(i,i) = d;
300 Dinv(i,i) = 1./d;
301 Dinv2(i,i) = 1./(d*d);
302 }
303
304 if (doTrustRegion) {
305 // do simple trust region using some control delta and eta
306 // compute norm of Newton step
307 double norm = 0;
308 if (scaleTR) {
309 print.Debug("scaling Trust region with diagonal matrix D ",D);
310 stepScaled = D * step;
312 }
313 else
314 norm = sqrt(inner_product(step,step));
315 // some conditions
316
317 // double tr_radius = norm;
318 // if (norm < 0.1 * delta) tr_radius = 0.1*delta;
319 // else if (norm > delta) tr_radius = delta;
320
321 // update new point using reduced step
322 //step = (tr_radius/norm) * step;
323
324 // if the proposed point (newton step) is inside the trust region radius accept it
325 if (norm <= delta) {
326 p = MinimumParameters(s0.Vec() + step, fcnCaller(s0.Vec() + step));
327 print.Debug("Accept full Newton step - it is inside TR ",delta);
328 } else {
329 //step = - (delta/norm) * step;
330
331 // if Newton step is outside try to use the Cauchy point ?
332 double normGrad2 = 0;
333 double gHg = 0;
334 if (scaleTR) {
335 auto gScaled = Dinv * s0.Gradient().Grad();
337 // compute D-1 H D-1 = H(i,j) * D(i,i) * D(j,j)
339 for (unsigned int i = 0; i < n; i++) {
340 for (unsigned int j = 0; j <=i; j++) {
341 Hscaled(i,j) = H(i,j) * Dinv(i,i) * Dinv(j,j);
342 }
343 }
345 } else {
346 normGrad2 = inner_product(s0.Gradient().Grad(),s0.Gradient().Grad());
347 gHg = similarity(s0.Gradient().Grad(), s0.Error().Hessian());
348 }
349 double normGrad = sqrt(normGrad2);
350 // need to compute gTHg :
351
352
353 print.Debug("computed gdel gHg and normGN",gdel,gHg,norm*norm, normGrad2);
354 if (gHg <= 0.) {
355 // Cauchy point is at the trust region boundary)
356 step = - (delta/ normGrad) * s0.Gradient().Grad();
357 if (scaleTR) step = Dinv2 * step;
358 print.Debug("Use as new point the Cauchy point - along gradient with norm=delta ", delta);
359 } else {
360 // Cauchy point can be inside the trust region
361 double tau = std::min( (normGrad2*normGrad)/(gHg*delta), 1.0);
362
364 stepC = -tau * (delta / normGrad) * s0.Gradient().Grad();
365 if (scaleTR)
366 stepC = Dinv2 * stepC;
367
368 print.Debug("Use as new point the Cauchy point - along gradient with tau ", tau, "delta = ", delta);
369
370 // compute dog -leg step solving quadratic equation a * t^2 + b * t + c = 0
371 // where a = ||p_n - P_c ||^2
372 // b = 2 * p_c * (P_n - P_c)
373 // c = || pc ||^2 - ||Delta||^2
374 auto diffP = step - stepC;
375 double a = inner_product(diffP, diffP);
376 double b = 2. * inner_product(stepC, diffP);
377 // norm cauchy point is tau*delta
378 double c = (scaleTR) ? inner_product(stepC, stepC) - delta * delta : delta * delta * (tau * tau - 1.);
379
380 print.Debug(" dogleg equation", a, b, c);
381 // solution is
382 double t = 0;
383 // a cannot be negative or zero, otherwise point was accepted
384 if (a <= 0) {
385 // case a is zero
386 print.Warn("a is equal to zero! a = ", a);
387 print.Info(" delta ", delta, " tau ", tau, " gHg ", gHg, " normgrad2 ", normGrad2);
388 t = -b / c;
389 } else {
390 double t1 = (-b + sqrt(b * b - 4. * a * c)) / (2.0 * a);
391 double t2 = (-b - sqrt(b * b - 4. * a * c)) / (2.0 * a);
392 // if b is positive solution is
393 print.Debug(" solution dogleg equation", t1, t2);
394 if (t1 >= 0 && t1 <= 1.)
395 t = t1;
396 else
397 t = t2;
398 }
399 step = stepC + t * diffP;
400 // need to rescale point per D^-1 >
401 print.Debug("New dogleg point is t = ", t);
402 }
403 print.Debug("New accepted step is ",step);
404
405 p = MinimumParameters(s0.Vec() + step, fcnCaller(s0.Vec() + step));
406 norm = delta;
407 gdel = inner_product(step, s0.Gradient().Grad());
408 }
409
410 // compute real changes (require an evaluation)
411
412 // expected change is gdel (can correct for Hessian )
413 //double fcexp = (-gdel - 0.5 * dot(step, hess(x) * step)
414
415 double svs = 0.5 * similarity(step, s0.Error().Hessian());
416 double rho = (p.Fval()-s0.Fval())/(gdel+svs);
417 if (rho < 0.25) {
418 delta = tr_factor_down * delta;
419 } else {
420 if (rho > 0.75 && norm == delta) {
421 delta = std::min(tr_factor_up * delta, 100.*norm); // 1000. is the delta max
422 }
423 }
424 print.Debug("New point after Trust region :", "norm tr ",norm," rho ", rho," delta ", delta,
425 " FVAL ", p.Fval(), "\n Parameter", p.Vec());
426
427 // check if we need to accept the point ?
428 acceptNewPoint = (rho > eta);
429 if (acceptNewPoint) {
430 // we accept the point
431 // we have already p = s0 + step
432 print.Debug("Trust region: accept new point p = x + step since rho is larger than eta");
433 }
434 else {
435 // we keep old point
436 print.Debug("Trust region reject new point and repeat since rho is smaller than eta");
437 p = MinimumParameters(s0.Vec(), s0.Fval());
438 }
439
440 }
441
442 FunctionGradient g = s0.Gradient();
443
444 if (acceptNewPoint || result.size() == 1) {
445
446 print.Debug("Before Gradient - NCalls = ", fcn.NumOfCalls());
447
448 g = gc(p, s0.Gradient());
449
450 print.Debug("After Gradient - NCalls = ", fcn.NumOfCalls());
451
452 }
453
454 // move Error updator after Gradient since the Value is cached inside
455
456 MinimumError e = ErrorUpdator().Update(s0, p, gc, lambda);
457
458 // should I use here new error instead of old one (so.Error) ?
459 edm = Estimator().Estimate(g, s0.Error());
460
461 print.Debug("Updated new point:", "\n FVAL ", p.Fval(), "\n Parameter", p.Vec(), "\n Gradient", g.Vec(),
462 "\n InvHessian", e.InvHessian(), "\n Hessian", e.Hessian(), "\n Edm", edm);
463
464 if (edm < 0.) {
465 print.Warn("Matrix not pos.def., Edm < 0");
466
468 s0 = psdf(s0, prec);
469 edm = Estimator().Estimate(g, s0.Error());
470 if (edm < 0.) {
471 result.push_back(s0);
472 if (TraceIter())
473 TraceIteration(result.size() - 1, result.back());
474 return FunctionMinimum(seed, result, fcn.Up());
475 }
476 }
477
478 // check lambda according to step
479 if (doLineSearch) {
480 if (p.Fval() < s0.Fval())
481 // fcn is decreasing along the step
482 lambda *= 0.1;
483 else {
484 lambda *= 10;
485 // if close to minimum stop to avoid oscillations around minimum value
486 if (edm < 0.1)
487 break;
488 }
489 }
490
491 print.Debug("finish iteration -", result.size(), "lambda =", lambda, "f1 =", p.Fval(), "f0 =", s0.Fval(),
492 "num of calls =", fcn.NumOfCalls(), "edm =", edm);
493
494 result.push_back(MinimumState(p, e, g, edm, fcn.NumOfCalls()));
495 if (TraceIter())
496 TraceIteration(result.size() - 1, result.back());
497
498 print.Info(MnPrint::Oneline(result.back(), result.size()));
499
500 edm *= (1. + 3. * e.Dcovar());
501
502 //} // endif on acceptNewPoint
503
504
505 } while (edm > edmval && fcn.NumOfCalls() < maxfcn);
506
507 if (fcn.NumOfCalls() >= maxfcn) {
508 print.Warn("Call limit exceeded",fcn.NumOfCalls(), maxfcn);
509
511 }
512
513 if (edm > edmval) {
514 if (edm < std::fabs(prec.Eps2() * result.back().Fval())) {
515 print.Warn("Machine accuracy limits further improvement");
516
517 return FunctionMinimum(seed, result, fcn.Up());
518 } else if (edm < 10 * edmval) {
519 return FunctionMinimum(seed, result, fcn.Up());
520 } else {
521
522 print.Warn("No convergence; Edm", edm, "is above tolerance", 10 * edmval);
523
525 }
526 }
527
528 print.Debug("Exiting successfully", "Ncalls", fcn.NumOfCalls(), "FCN", result.back().Fval(), "Edm", edm, "Requested",
529 edmval);
530
531 return FunctionMinimum(seed, result, fcn.Up());
532}
533
534} // namespace Minuit2
535
536} // namespace ROOT
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define s0(x)
Definition RSha256.hxx:90
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
FunctionMinimum Minimum(const MnFcn &fMnFcn, const GradientCalculator &fGradienCalculator, const MinimumSeed &fMinimumSeed, const MnStrategy &fMnStrategy, unsigned int maxfcn, double edmval) const override
Class the member function calculating the Minimum and verifies the result depending on the strategy.
const FumiliErrorUpdator & ErrorUpdator() const
Accessor to the Error updator of the builder.
const VariableMetricEDMEstimator & Estimator() const
Accessor to the EDM (expected vertical distance to the Minimum) estimator.
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
interface class for gradient calculators
Class describing a symmetric matrix of size n.
Definition MnMatrix.h:438
unsigned int size() const
Definition MnMatrix.h:1032
void TraceIteration(int iter, const MinimumState &state) const
MinimumError keeps the inv.
const FunctionGradient & Gradient() const
Definition MinimumSeed.h:32
const MinimumError & Error() const
Definition MinimumSeed.h:31
const MinimumParameters & Parameters() const
Definition MinimumSeed.h:30
const MnMachinePrecision & Precision() const
Definition MinimumSeed.h:34
const MinimumState & State() const
Definition MinimumSeed.h:29
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Wrapper class to FCNBase interface used internally by Minuit.
Definition MnFcn.h:34
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition MnHesse.h:41
Implements a 1-dimensional minimization along a given direction (i.e.
Sets the relative floating point (double) arithmetic precision.
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:135
void Error(const Ts &... args)
Definition MnPrint.h:117
void Info(const Ts &... args)
Definition MnPrint.h:129
void Warn(const Ts &... args)
Definition MnPrint.h:123
API class for defining four levels of strategies: low (0), medium (1), high (2), very high (>=3); act...
Definition MnStrategy.h:27
const Int_t n
Definition legend1.C:16
#define H(x, y, z)
double similarity(const LAVector &, const LASymMatrix &)
Definition MnMatrix.cxx:629
double inner_product(const LAVector &, const LAVector &)
Definition MnMatrix.cxx:316
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
auto * t1
Definition textangle.C:20