Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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#include "Minuit2/MnPrint.h"
27
28#include <cmath>
29#include <cassert>
30
31namespace ROOT {
32
33namespace Minuit2 {
34
35double inner_product(const LAVector &, const LAVector &);
36
37void VariableMetricBuilder::AddResult(std::vector<MinimumState> &result, const MinimumState &state) const
38{
39 // // if (!store) store = StorageLevel();
40 // // store |= (result.size() == 0);
41 // if (store)
42 result.push_back(state);
43 // else {
44 // result.back() = state;
45 // }
46 if (TraceIter())
47 TraceIteration(result.size() - 1, result.back());
48 else {
49 MnPrint print("VariableMetricBuilder", PrintLevel());
50 print.Info(MnPrint::Oneline(result.back(), result.size() - 1));
51 }
52}
53
55 const MnStrategy &strategy, unsigned int maxfcn, double edmval) const
56{
57 MnPrint print("VariableMetricBuilder", PrintLevel());
58
59 // top level function to find minimum from a given initial seed
60 // iterate on a minimum search in case of first attempt is not successful
61
62 // to be consistent with F77 Minuit
63 // in Minuit2 edm is correct and is ~ a factor of 2 smaller than F77Minuit
64 // There are also a check for convergence if (edm < 0.1 edmval for exiting the loop)
65 // LM: change factor to 2E-3 to be consistent with F77Minuit
66 edmval *= 0.002;
67
68 // set global printlevel to the local one so all calls to MN_INFO_MSG can be controlled in the same way
69 // at exit of this function the BuilderPrintLevelConf object is destructed and automatically the
70 // previous level will be restored
71
72 // double edm = Estimator().Estimate(seed.Gradient(), seed.Error());
73 double edm = seed.State().Edm();
74
75 FunctionMinimum min(seed, fcn.Up());
76
77 if (seed.Parameters().Vec().size() == 0) {
78 print.Warn("No free parameters.");
79 return min;
80 }
81
82 if (!seed.IsValid()) {
83 print.Error("Minimum seed invalid.");
84 return min;
85 }
86
87 if (edm < 0.) {
88 print.Error("Initial matrix not pos.def.");
89
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 // do actual iterations
101 print.Info("Start iterating until Edm is <", edmval, "with call limit =", maxfcn);
102
103 AddResult(result, seed.State());
104
105 // try first with a maxfxn = 80% of maxfcn
106 int maxfcn_eff = maxfcn;
107 int ipass = 0;
108 bool iterate = false;
109
110 do {
111
112 iterate = false;
113
114 print.Debug(ipass > 0 ? "Continue" : "Start", "iterating...");
115
116 min = Minimum(fcn, gc, seed, result, maxfcn_eff, edmval);
117
118 // if max function call reached exits
119 if (min.HasReachedCallLimit()) {
120 print.Warn("FunctionMinimum is invalid, reached function call limit");
121 return min;
122 }
123
124 // second time check for validity of function Minimum
125 if (ipass > 0) {
126 if (!min.IsValid()) {
127 print.Warn("FunctionMinimum is invalid after second try");
128 return min;
129 }
130 }
131
132 // resulting edm of minimization
133 edm = result.back().Edm();
134 // need to correct again for Dcovar: edm *= (1. + 3. * e.Dcovar()) ???
135
136 if ((strategy.Strategy() >= 2) || (strategy.Strategy() == 1 && min.Error().Dcovar() > 0.05)) {
137
138 print.Debug("MnMigrad will verify convergence and Error matrix; dcov =", min.Error().Dcovar());
139
140 MnStrategy strat(strategy);
141 strat.SetHessianForcePosDef(1); // ensure no matter what strategy is used, we force the result positive-definite if required
142 MinimumState st = MnHesse(strat)(fcn, min.State(), min.Seed().Trafo(), maxfcn);
143
144 print.Info("After Hessian");
145
146 AddResult(result, st);
147
148 if (!st.IsValid()) {
149 print.Warn("Invalid Hessian - exit the minimization");
150 break;
151 }
152
153 // check new edm
154 edm = st.Edm();
155
156 print.Debug("New Edm", edm, "Requested", edmval);
157
158 if (edm > edmval) {
159 // be careful with machine precision and avoid too small edm
160 double machineLimit = std::fabs(seed.Precision().Eps2() * result.back().Fval());
161 if (edm >= machineLimit) {
162 iterate = true;
163
164 print.Info("Tolerance not sufficient, continue minimization; "
165 "Edm",
166 edm, "Required", edmval);
167 } else {
168 print.Warn("Reached machine accuracy limit; Edm", edm, "is smaller than machine limit", machineLimit,
169 "while", edmval, "was requested");
170 }
171 }
172 }
173
174 // end loop on iterations
175 // ? need a maximum here (or max of function calls is enough ? )
176 // continnue iteration (re-calculate function Minimum if edm IS NOT sufficient)
177 // no need to check that hesse calculation is done (if isnot done edm is OK anyway)
178 // count the pass to exit second time when function Minimum is invalid
179 // increase by 20% maxfcn for doing some more tests
180 if (ipass == 0)
181 maxfcn_eff = int(maxfcn * 1.3);
182 ipass++;
183 } while (iterate);
184
185 // Add latest state (Hessian calculation)
186 const MinimumState &latest = result.back();
187
188 // check edm (add a factor of 10 in tolerance )
189 if (edm > 10 * edmval) {
190 min.Add(latest, FunctionMinimum::MnAboveMaxEdm);
191 print.Warn("No convergence; Edm", edm, "is above tolerance", 10 * edmval);
192 } else if (latest.Error().HasReachedCallLimit()) {
193 // communicate to user that call limit was reached in MnHesse
195 } else if (latest.Error().IsAvailable()) {
196 // check if minimum had edm above max before
197 if (min.IsAboveMaxEdm())
198 print.Info("Edm has been re-computed after Hesse; Edm", edm, "is now within tolerance");
199 min.Add(latest);
200 }
201
202 print.Debug("Minimum found", min);
203
204 return min;
205}
206
208 std::vector<MinimumState> &result, unsigned int maxfcn,
209 double edmval) const
210{
211 // function performing the minimum searches using the Variable Metric algorithm (MIGRAD)
212 // perform first a line search in the - Vg direction and then update using the Davidon formula (Davidon Error
213 // updator) stop when edm reached is less than required (edmval)
214
215 // after the modification when I iterate on this functions, so it can be called many times,
216 // the seed is used here only to get precision and construct the returned FunctionMinimum object
217
218 MnPrint print("VariableMetricBuilder", PrintLevel());
219
220 const MnMachinePrecision &prec = seed.Precision();
221
222 // result.push_back(MinimumState(seed.Parameters(), seed.Error(), seed.Gradient(), edm, fcn.NumOfCalls()));
223 const MinimumState &initialState = result.back();
224
225 double edm = initialState.Edm();
226
227 print.Debug("Initial State:", "\n Parameter:", initialState.Vec(), "\n Gradient:", initialState.Gradient().Vec(),
228 "\n InvHessian:", initialState.Error().InvHessian(), "\n Edm:", initialState.Edm());
229
230 // iterate until edm is small enough or max # of iterations reached
231 edm *= (1. + 3. * initialState.Error().Dcovar());
232 MnLineSearch lsearch;
233 MnAlgebraicVector step(initialState.Gradient().Vec().size());
234 // keep also prevStep
235 MnAlgebraicVector prevStep(initialState.Gradient().Vec().size());
236
237 MinimumState s0 = result.back();
238
239 do {
240
241 // MinimumState s0 = result.back();
242
243 step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
244
245 print.Debug("Iteration", result.size(), "Fval", s0.Fval(), "numOfCall", fcn.NumOfCalls(),
246 "\n Internal parameters", s0.Vec(), "\n Newton step", step);
247
248 // check if derivatives are not zero
249 if (inner_product(s0.Gradient().Vec(), s0.Gradient().Vec()) <= 0) {
250 print.Debug("all derivatives are zero - return current status");
251 break;
252 }
253
254 // gdel = s^T * g = -g^T H g (since s = - Hg) so it must be negative
255 double gdel = inner_product(step, s0.Gradient().Grad());
256
257 if (gdel > 0.) {
258 print.Warn("Matrix not pos.def, gdel =", gdel, "> 0");
259
260 MnPosDef psdf;
261 s0 = psdf(s0, prec);
262 step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
263 // #ifdef DEBUG
264 // std::cout << "After MnPosdef - Error " << s0.Error().InvHessian() << " Gradient " <<
265 // s0.Gradient().Vec() << " step " << step << std::endl;
266 // #endif
267 gdel = inner_product(step, s0.Gradient().Grad());
268
269 print.Warn("gdel =", gdel);
270
271 if (gdel > 0.) {
273
274 return FunctionMinimum(seed, result, fcn.Up());
275 }
276 }
277
278 MnParabolaPoint pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
279
280 // <= needed for case 0 <= 0
281 if (std::fabs(pp.Y() - s0.Fval()) <= std::fabs(s0.Fval()) * prec.Eps()) {
282
283 print.Warn("No improvement in line search");
284
285 // no improvement exit (is it really needed LM ? in vers. 1.22 tried alternative )
286 // add new state when only fcn changes
287 if (result.size() <= 1)
288 AddResult(result, MinimumState(s0.Parameters(), s0.Error(), s0.Gradient(), s0.Edm(), fcn.NumOfCalls()));
289 else
290 // no need to re-store the state
291 AddResult(result, MinimumState(pp.Y(), s0.Edm(), fcn.NumOfCalls()));
292
293 break;
294 }
295
296 print.Debug("Result after line search :", "\n x =", pp.X(), "\n Old Fval =", s0.Fval(),
297 "\n New Fval =", pp.Y(), "\n NFcalls =", fcn.NumOfCalls());
298
299 MinimumParameters p(s0.Vec() + pp.X() * step, pp.Y());
300
301 FunctionGradient g = gc(p, s0.Gradient());
302
303 edm = Estimator().Estimate(g, s0.Error());
304
305 if (std::isnan(edm)) {
306 print.Warn("Edm is NaN; stop iterations");
308 return FunctionMinimum(seed, result, fcn.Up());
309 }
310
311 if (edm < 0.) {
312 print.Warn("Matrix not pos.def., try to make pos.def.");
313
314 MnPosDef psdf;
315 s0 = psdf(s0, prec);
316 edm = Estimator().Estimate(g, s0.Error());
317 if (edm < 0.) {
318 print.Warn("Matrix still not pos.def.; stop iterations");
319
321
322 return FunctionMinimum(seed, result, fcn.Up());
323 }
324 }
326
327 // avoid print Hessian that will invert the matrix
328 print.Debug("Updated new point:", "\n Parameter:", p.Vec(), "\n Gradient:", g.Vec(),
329 "\n InvHessian:", e.Matrix(), "\n Edm:", edm);
330
331 // update the state
332 s0 = MinimumState(p, e, g, edm, fcn.NumOfCalls());
333 if (StorageLevel() || result.size() <= 1)
335 else
336 // use a reduced state for not-final iterations
337 AddResult(result, MinimumState(p.Fval(), edm, fcn.NumOfCalls()));
338
339 // correct edm
340 edm *= (1. + 3. * e.Dcovar());
341
342 print.Debug("Dcovar =", e.Dcovar(), "\tCorrected edm =", edm);
343
344 } while (edm > edmval && fcn.NumOfCalls() < maxfcn); // end of iteration loop
345
346 // save last result in case of no complete final states
347 // when the result is filled above (reduced storage) the resulting state will not be valid
348 // since they will not have parameter values and error
349 // the line above will fill as last element a valid state
350 if (!result.back().IsValid())
351 result.back() = s0;
352
353 if (fcn.NumOfCalls() >= maxfcn) {
354 print.Warn("Call limit exceeded");
356 }
357
358 if (edm > edmval) {
359 if (edm < 10 * edmval) {
360 print.Info("Edm is close to limit - return current minimum");
361 return FunctionMinimum(seed, result, fcn.Up());
362 } else if (edm < std::fabs(prec.Eps2() * result.back().Fval())) {
363 print.Warn("Edm is limited by Machine accuracy - return current minimum");
364 return FunctionMinimum(seed, result, fcn.Up());
365 } else {
366 print.Warn("Iterations finish without convergence; Edm", edm, "Requested", edmval);
367
369 }
370 }
371 // std::cout<<"result.back().Error().Dcovar()= "<<result.back().Error().Dcovar()<<std::endl;
372
373 print.Debug("Exiting successfully;", "Ncalls", fcn.NumOfCalls(), "FCN", result.back().Fval(), "Edm", edm,
374 "Requested", edmval);
375
376 return FunctionMinimum(seed, result, fcn.Up());
377}
378
379} // namespace Minuit2
380
381} // namespace ROOT
#define g(i)
Definition RSha256.hxx:105
#define s0(x)
Definition RSha256.hxx:90
#define e(i)
Definition RSha256.hxx:103
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
const MnAlgebraicVector & Vec() const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
interface class for gradient calculators
unsigned int size() const
Definition LAVector.h:231
void TraceIteration(int iter, const MinimumState &state) const
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const =0
MinimumError keeps the inv.
const MnAlgebraicSymMatrix & InvHessian() const
const MnAlgebraicVector & Vec() const
const MinimumParameters & Parameters() const
Definition MinimumSeed.h:29
const MnMachinePrecision & Precision() const
Definition MinimumSeed.h:33
const MinimumState & State() const
Definition MinimumSeed.h:28
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
const MinimumError & Error() const
const MnAlgebraicVector & Vec() const
const FunctionGradient & Gradient() const
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
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.
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)
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:25
void Debug(const Ts &... args)
Definition MnPrint.h:147
void Error(const Ts &... args)
Definition MnPrint.h:129
void Info(const Ts &... args)
Definition MnPrint.h:141
void Warn(const Ts &... args)
Definition MnPrint.h:135
API class for defining four levels of strategies: low (0), medium (1), high (2), very high (>=3); act...
Definition MnStrategy.h:27
unsigned int Strategy() const
Definition MnStrategy.h:38
void SetHessianForcePosDef(unsigned int flag)
Definition MnStrategy.h:78
void AddResult(std::vector< MinimumState > &result, const MinimumState &state) const
FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const override
const VariableMetricEDMEstimator & Estimator() const
const MinimumErrorUpdator & ErrorUpdator() const
double Estimate(const FunctionGradient &, const MinimumError &) const
int iterate(rng_state_t *X)
Definition mixmax.icc:34
double inner_product(const LAVector &, const LAVector &)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.