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#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 MinimumState st = MnHesse(strategy)(fcn, min.State(), min.Seed().Trafo(), maxfcn);
141
142 print.Info("After Hessian");
143
144 AddResult(result, st);
145
146 if (!st.IsValid()) {
147 print.Warn("Invalid Hessian - exit the minimization");
148 break;
149 }
150
151 // check new edm
152 edm = st.Edm();
153
154 print.Debug("New Edm", edm, "Requested", edmval);
155
156 if (edm > edmval) {
157 // be careful with machine precision and avoid too small edm
158 double machineLimit = std::fabs(seed.Precision().Eps2() * result.back().Fval());
159 if (edm >= machineLimit) {
160 iterate = true;
161
162 print.Info("Tolerance not sufficient, continue minimization; "
163 "Edm",
164 edm, "Required", edmval);
165 } else {
166 print.Warn("Reached machine accuracy limit; Edm", edm, "is smaller than machine limit", machineLimit,
167 "while", edmval, "was requested");
168 }
169 }
170 }
171
172 // end loop on iterations
173 // ? need a maximum here (or max of function calls is enough ? )
174 // continnue iteration (re-calculate function Minimum if edm IS NOT sufficient)
175 // no need to check that hesse calculation is done (if isnot done edm is OK anyway)
176 // count the pass to exit second time when function Minimum is invalid
177 // increase by 20% maxfcn for doing some more tests
178 if (ipass == 0)
179 maxfcn_eff = int(maxfcn * 1.3);
180 ipass++;
181 } while (iterate);
182
183 // Add latest state (Hessian calculation)
184 const MinimumState &latest = result.back();
185
186 // check edm (add a factor of 10 in tolerance )
187 if (edm > 10 * edmval) {
189 print.Warn("No convergence; Edm", edm, "is above tolerance", 10 * edmval);
190 } else if (latest.Error().HasReachedCallLimit()) {
191 // communicate to user that call limit was reached in MnHesse
193 } else if (latest.Error().IsAvailable()) {
194 // check if minimum had edm above max before
195 if (min.IsAboveMaxEdm())
196 print.Info("Edm has been re-computed after Hesse; Edm", edm, "is now within tolerance");
197 min.Add(latest);
198 }
199
200 print.Debug("Minimum found", min);
201
202 return min;
203}
204
206 std::vector<MinimumState> &result, unsigned int maxfcn,
207 double edmval) const
208{
209 // function performing the minimum searches using the Variable Metric algorithm (MIGRAD)
210 // perform first a line search in the - Vg direction and then update using the Davidon formula (Davidon Error
211 // updator) stop when edm reached is less than required (edmval)
212
213 // after the modification when I iterate on this functions, so it can be called many times,
214 // the seed is used here only to get precision and construct the returned FunctionMinimum object
215
216 MnPrint print("VariableMetricBuilder", PrintLevel());
217
218 const MnMachinePrecision &prec = seed.Precision();
219
220 // result.push_back(MinimumState(seed.Parameters(), seed.Error(), seed.Gradient(), edm, fcn.NumOfCalls()));
221 const MinimumState &initialState = result.back();
222
223 double edm = initialState.Edm();
224
225 print.Debug("Initial State:", "\n Parameter:", initialState.Vec(), "\n Gradient:", initialState.Gradient().Vec(),
226 "\n InvHessian:", initialState.Error().InvHessian(), "\n Edm:", initialState.Edm());
227
228 // iterate until edm is small enough or max # of iterations reached
229 edm *= (1. + 3. * initialState.Error().Dcovar());
230 MnLineSearch lsearch;
231 MnAlgebraicVector step(initialState.Gradient().Vec().size());
232 // keep also prevStep
233 MnAlgebraicVector prevStep(initialState.Gradient().Vec().size());
234
235 MinimumState s0 = result.back();
236
237 do {
238
239 // MinimumState s0 = result.back();
240
241 step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
242
243 print.Debug("Iteration", result.size(), "Fval", s0.Fval(), "numOfCall", fcn.NumOfCalls(),
244 "\n Internal parameters", s0.Vec(), "\n Newton step", step);
245
246 // check if derivatives are not zero
247 if (inner_product(s0.Gradient().Vec(), s0.Gradient().Vec()) <= 0) {
248 print.Debug("all derivatives are zero - return current status");
249 break;
250 }
251
252 double gdel = inner_product(step, s0.Gradient().Grad());
253
254 if (gdel > 0.) {
255 print.Warn("Matrix not pos.def, gdel =", gdel, "> 0");
256
257 MnPosDef psdf;
258 s0 = psdf(s0, prec);
259 step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
260 // #ifdef DEBUG
261 // std::cout << "After MnPosdef - Error " << s0.Error().InvHessian() << " Gradient " <<
262 // s0.Gradient().Vec() << " step " << step << std::endl;
263 // #endif
264 gdel = inner_product(step, s0.Gradient().Grad());
265
266 print.Warn("gdel =", gdel);
267
268 if (gdel > 0.) {
270
271 return FunctionMinimum(seed, result, fcn.Up());
272 }
273 }
274
275 MnParabolaPoint pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
276
277 // <= needed for case 0 <= 0
278 if (std::fabs(pp.Y() - s0.Fval()) <= std::fabs(s0.Fval()) * prec.Eps()) {
279
280 print.Warn("No improvement in line search");
281
282 // no improvement exit (is it really needed LM ? in vers. 1.22 tried alternative )
283 // add new state when only fcn changes
284 if (result.size() <= 1)
285 AddResult(result, MinimumState(s0.Parameters(), s0.Error(), s0.Gradient(), s0.Edm(), fcn.NumOfCalls()));
286 else
287 // no need to re-store the state
288 AddResult(result, MinimumState(pp.Y(), s0.Edm(), fcn.NumOfCalls()));
289
290 break;
291 }
292
293 print.Debug("Result after line search :", "\n x =", pp.X(), "\n Old Fval =", s0.Fval(),
294 "\n New Fval =", pp.Y(), "\n NFcalls =", fcn.NumOfCalls());
295
296 MinimumParameters p(s0.Vec() + pp.X() * step, pp.Y());
297
298 FunctionGradient g = gc(p, s0.Gradient());
299
300 edm = Estimator().Estimate(g, s0.Error());
301
302 if (std::isnan(edm)) {
303 print.Warn("Edm is NaN; stop iterations");
305 return FunctionMinimum(seed, result, fcn.Up());
306 }
307
308 if (edm < 0.) {
309 print.Warn("Matrix not pos.def., try to make pos.def.");
310
311 MnPosDef psdf;
312 s0 = psdf(s0, prec);
313 edm = Estimator().Estimate(g, s0.Error());
314 if (edm < 0.) {
315 print.Warn("Matrix still not pos.def.; stop iterations");
316
318
319 return FunctionMinimum(seed, result, fcn.Up());
320 }
321 }
323
324 // avoid print Hessian that will invert the matrix
325 print.Debug("Updated new point:", "\n Parameter:", p.Vec(), "\n Gradient:", g.Vec(),
326 "\n InvHessian:", e.Matrix(), "\n Edm:", edm);
327
328 // update the state
329 s0 = MinimumState(p, e, g, edm, fcn.NumOfCalls());
330 if (StorageLevel() || result.size() <= 1)
332 else
333 // use a reduced state for not-final iterations
334 AddResult(result, MinimumState(p.Fval(), edm, fcn.NumOfCalls()));
335
336 // correct edm
337 edm *= (1. + 3. * e.Dcovar());
338
339 print.Debug("Dcovar =", e.Dcovar(), "\tCorrected edm =", edm);
340
341 } while (edm > edmval && fcn.NumOfCalls() < maxfcn); // end of iteration loop
342
343 // save last result in case of no complete final states
344 // when the result is filled above (reduced storage) the resulting state will not be valid
345 // since they will not have parameter values and error
346 // the line above will fill as last element a valid state
347 if (!result.back().IsValid())
348 result.back() = s0;
349
350 if (fcn.NumOfCalls() >= maxfcn) {
351 print.Warn("Call limit exceeded");
353 }
354
355 if (edm > edmval) {
356 if (edm < 10 * edmval) {
357 print.Info("Edm is close to limit - return current minimum");
358 return FunctionMinimum(seed, result, fcn.Up());
359 } else if (edm < std::fabs(prec.Eps2() * result.back().Fval())) {
360 print.Warn("Edm is limited by Machine accuracy - return current minimum");
361 return FunctionMinimum(seed, result, fcn.Up());
362 } else {
363 print.Warn("Iterations finish without convergence; Edm", edm, "Requested", edmval);
364
366 }
367 }
368 // std::cout<<"result.back().Error().Dcovar()= "<<result.back().Error().Dcovar()<<std::endl;
369
370 print.Debug("Exiting successfully;", "Ncalls", fcn.NumOfCalls(), "FCN", result.back().Fval(), "Edm", edm,
371 "Requested", edmval);
372
373 return FunctionMinimum(seed, result, fcn.Up());
374}
375
376} // namespace Minuit2
377
378} // namespace ROOT
#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 char Point_t Rectangle_t WindowAttributes_t Float_t Float_t g
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) ...
void Add(const MinimumState &state, Status status=MnValid)
add latest minimization state (for example add Hesse result after Migrad)
const MinimumError & Error() const
const MinimumState & State() const
const MinimumSeed & Seed() const
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.
Definition: MinimumError.h:28
bool HasReachedCallLimit() const
Definition: MinimumError.h:86
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:50
const MnAlgebraicVector & Vec() const
const MnUserTransformation & Trafo() const
Definition: MinimumSeed.h:32
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 (...
Definition: MinimumState.h:27
const MinimumError & Error() const
Definition: MinimumState.h:61
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:58
const FunctionGradient & Gradient() const
Definition: MinimumState.h:62
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.
Definition: MnLineSearch.h:40
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: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 three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
unsigned int Strategy() const
Definition: MnStrategy.h:38
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
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 &)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.