Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
NegativeG2LineSearch.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
11#include "Minuit2/MnFcn.h"
18#include "Minuit2/MnPrint.h"
19
20#include "Math/Util.h"
21
22#include <cmath>
23
24namespace ROOT {
25
26namespace Minuit2 {
27
29 const MnMachinePrecision &prec) const
30{
31
32 // when the second derivatives are negative perform a line search along Parameter which gives
33 // negative second derivative and magnitude equal to the Gradient step size.
34 // Recalculate the gradients for all the Parameter after the correction and
35 // continue iteration in case the second derivatives are still negative
36 //
37 MnPrint print("NegativeG2LineSearch");
38
39 // Print the runtime on returning from the function
40 ROOT::Math::Util::TimingScope timingScope([&print](std::string const& s){ print.Info(s); }, "Done after");
41
42 bool negG2 = HasNegativeG2(st.Gradient(), prec);
43 if (!negG2)
44 return st;
45
46 print.Info("Doing a NegativeG2LineSearch since one of the G2 component is negative");
47
48 unsigned int n = st.Parameters().Vec().size();
49 FunctionGradient dgrad = st.Gradient();
50 MinimumParameters pa = st.Parameters();
51 bool iterate = false;
52 unsigned int iter = 0;
53 // in case of analytical gradients we don't have step sizes
54 bool hasGStep = !dgrad.IsAnalytical();
55 //print.Trace("Gradient ", dgrad.Vec(), "G2",dgrad.G2());
56 // gradient present in the state must have G2 otherwise something is wrong
57 if (!dgrad.HasG2()) {
58 print.Error("Input gradient to NG2LS must have G2 already computed");
59 return st;
60 }
61
62 do {
63 iterate = false;
64 for (unsigned int i = 0; i < n; i++) {
65
66 if (dgrad.G2()(i) <= 0) {
67
68 // check also the gradient (if it is zero ) I can skip the param)
69
70 if (std::fabs(dgrad.Vec()(i)) < prec.Eps() && std::fabs(dgrad.G2()(i)) < prec.Eps())
71 continue;
72 // if(dgrad.G2()(i) < prec.Eps()) {
73 // do line search if second derivative negative
74 MnAlgebraicVector step(n);
76
77 // when using analytical gradient use as step size a dummy value of 1
78 // maybe could do better using user given parameter step sizes
79 // tested using inverse of G2() gives worse behaviour
80 if (dgrad.Vec()(i) < 0)
81 step(i) = (hasGStep) ? dgrad.Gstep()(i) : 1;
82 else
83 step(i) = (hasGStep) ? -dgrad.Gstep()(i) : -1;
84
85 double gdel = step(i) * dgrad.Vec()(i);
86
87 // if using sec derivative information
88 // double g2del = step(i)*step(i) * dgrad.G2()(i);
89
90 print.Debug("Iter", iter, "param", i, pa.Vec()(i), "grad2", dgrad.G2()(i), "grad",
91 dgrad.Vec()(i), "grad step", step(i), " gdel ", gdel);
92
93 MnParabolaPoint pp = lsearch(fcn, pa, step, gdel, prec);
94
95 print.Debug("Line search result", pp.X(), "f(0)", pa.Fval(), "f(1)", pp.Y());
96
97 step *= pp.X();
98 pa = MinimumParameters(pa.Vec() + step, pp.Y());
99
100 dgrad = gc(pa, dgrad);
101 // re-compute also G2 if needed
102 if (!dgrad.HasG2()) {
103 //no need to compute Hessian here but only G2
104 print.Debug("Compute G2 at the new point", pa.Vec());
106 bool ret = gc.G2(pa,g2);
107 if (!ret) {
108 print.Error("Cannot compute G2");
109 assert(false);
110 return st;
111 }
112
113 dgrad = FunctionGradient(dgrad.Grad(), g2);
114 }
115
116 print.Debug("New result after Line search - iter", iter, "param", i, pa.Vec()(i), "step", step(i), "new grad2",
117 dgrad.G2()(i), "new grad", dgrad.Vec()(i));
118
119 iterate = true;
120 break;
121 }
122 }
123 } while (iter++ < 2 * n && iterate);
124
125 // even if we computed the Hessian it is still better to use the diagonal part, the G2
126 print.Debug("Approximate new covariance after NegativeG2LS using only G2");
128 for (unsigned int i = 0; i < n; i++) {
129 mat(i, i) = std::fabs(dgrad.G2()(i)) > prec.Eps() ? 1. / dgrad.G2()(i) :
130 1; // use an arbitrary value (e.g. 1)
131 }
132
133 MinimumError err(mat, 1.);
134 double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
135
136 if (edm < 0) {
138 }
139
140 return MinimumState(pa, err, dgrad, edm, fcn.NumOfCalls());
141}
142
144{
145 // check if function gradient has any component which is negative
146
147 for (unsigned int i = 0; i < grad.Vec().size(); i++)
148
149 if (grad.G2()(i) <= 0) {
150
151 return true;
152 }
153
154 return false;
155}
156
157} // namespace Minuit2
158
159} // namespace ROOT
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
const MnAlgebraicVector & Vec() const
const MnAlgebraicVector & G2() const
interface class for gradient calculators
Class describing a symmetric matrix of size n.
Definition MnMatrix.h:438
MinimumError keeps the inv.
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
Implements a 1-dimensional minimization along a given direction (i.e.
Sets the relative floating point (double) arithmetic precision.
double Y() const
Accessor to the y (second) coordinate.
double X() const
Accessor to the x (first) coordinate.
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
bool HasNegativeG2(const FunctionGradient &, const MnMachinePrecision &) const
MinimumState operator()(const MnFcn &, const MinimumState &, const GradientCalculator &, const MnMachinePrecision &) const
double Estimate(const FunctionGradient &, const MinimumError &) const
const Int_t n
Definition legend1.C:16
int iterate(rng_state_t *X)
Definition mixmax.icc:34
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...