Logo ROOT  
Reference Guide
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
19#include <cmath>
20//#define DEBUG
21#ifdef DEBUG
22#include <iostream>
23#endif
24
25namespace ROOT {
26
27 namespace Minuit2 {
28
29
30
31
33
34// when the second derivatives are negative perform a line search along Parameter which gives
35// negative second derivative and magnitude equal to the Gradient step size.
36// Recalculate the gradients for all the Parameter after the correction and
37// continue iteration in case the second derivatives are still negative
38//
39
40 bool negG2 = HasNegativeG2(st.Gradient(), prec);
41 if(!negG2) return st;
42
43 unsigned int n = st.Parameters().Vec().size();
44 FunctionGradient dgrad = st.Gradient();
46 bool iterate = false;
47 unsigned int iter = 0;
48 do {
49 iterate = false;
50 for(unsigned int i = 0; i < n; i++) {
51
52#ifdef DEBUG
53 std::cout << "negative G2 - iter " << iter << " param " << i << " " << pa.Vec()(i) << " grad2 " << dgrad.G2()(i) << " grad " << dgrad.Vec()(i)
54 << " grad step " << dgrad.Gstep()(i) << " step size " << pa.Dirin()(i) << std::endl;
55#endif
56 if(dgrad.G2()(i) <= 0) {
57
58 // check also the gradient (if it is zero ) I can skip the param)
59
60 if ( std::fabs(dgrad.Vec()(i) ) < prec.Eps() && std::fabs(dgrad.G2()(i) ) < prec.Eps() ) continue;
61 // if(dgrad.G2()(i) < prec.Eps()) {
62 // do line search if second derivative negative
63 MnAlgebraicVector step(n);
64 MnLineSearch lsearch;
65
66 if ( dgrad.Vec()(i) < 0)
67 step(i) = dgrad.Gstep()(i); //*dgrad.Vec()(i);
68 else
69 step(i) = - dgrad.Gstep()(i); // *dgrad.Vec()(i);
70
71 double gdel = step(i)*dgrad.Vec()(i);
72
73 // if using sec derivative information
74 // double g2del = step(i)*step(i) * dgrad.G2()(i);
75 bool debugLS = false;
76
77#ifdef DEBUG
78 std::cout << "step(i) " << step(i) << " gdel " << gdel << std::endl;
79// std::cout << " g2del " << g2del << std::endl;
80 debugLS = true;
81#endif
82 MnParabolaPoint pp = lsearch(fcn, pa, step, gdel, prec,debugLS);
83
84
85
86#ifdef DEBUG
87 std::cout << "\nLine search result " << pp.X() << " f(0) " << pa.Fval() << " f(1) " << pp.Y() << std::endl;
88#endif
89
90 step *= pp.X();
91 pa = MinimumParameters(pa.Vec() + step, pp.Y());
92
93 dgrad = gc(pa, dgrad);
94
95#ifdef DEBUG
96 std::cout << "Line search - iter" << iter << " param " << i << " " << pa.Vec()(i) << " step " << step(i) << " new grad2 " << dgrad.G2()(i) << " new grad " << dgrad.Vec()(i) << " grad step " << dgrad.Gstep()(i) << std::endl;
97#endif
98
99 iterate = true;
100 break;
101 }
102 }
103 } while(iter++ < 2*n && iterate);
104
106 for(unsigned int i = 0; i < n; i++)
107 mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
108
109 MinimumError err(mat, 1.);
110 double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
111
112 if (edm < 0) {
114 }
115
116 return MinimumState(pa, err, dgrad, edm, fcn.NumOfCalls());
117}
118
120 // check if function gradient has any component which is neegative
121
122 for(unsigned int i = 0; i < grad.Vec().size(); i++)
123
124 if(grad.G2()(i) <= 0 ) {
125
126 return true;
127 }
128
129 return false;
130}
131
132
133
134
135
136
137
138 } // namespace Minuit2
139
140} // namespace ROOT
const MnAlgebraicVector & Gstep() const
const MnAlgebraicVector & Vec() const
const MnAlgebraicVector & G2() const
interface class for gradient calculators
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int size() const
Definition: LAVector.h:198
MinimumError keeps the inv.
Definition: MinimumError.h:26
const MnAlgebraicVector & Dirin() const
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
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.
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
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
VSD Structures.
Definition: StringConv.hxx:21