Logo ROOT   6.08/07
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"
12 #include "Minuit2/MinimumState.h"
15 #include "Minuit2/MnLineSearch.h"
18 
19 #include <cmath>
20 //#define DEBUG
21 #ifdef DEBUG
22 #include <iostream>
23 #endif
24 
25 namespace 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();
45  MinimumParameters pa = st.Parameters();
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 
105  MnAlgebraicSymMatrix mat(n);
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  return MinimumState(pa, err, dgrad, edm, fcn.NumOfCalls());
113 }
114 
115 bool NegativeG2LineSearch::HasNegativeG2(const FunctionGradient& grad, const MnMachinePrecision& /*prec */ ) const {
116  // check if function gradient has any component which is neegative
117 
118  for(unsigned int i = 0; i < grad.Vec().size(); i++)
119 
120  if(grad.G2()(i) <= 0 ) {
121 
122  return true;
123  }
124 
125  return false;
126 }
127 
128 
129 
130 
131 
132 
133 
134  } // namespace Minuit2
135 
136 } // namespace ROOT
double Y() const
Accessor to the y (second) coordinate.
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
double Estimate(const FunctionGradient &, const MinimumError &) const
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
bool HasNegativeG2(const FunctionGradient &, const MnMachinePrecision &) const
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
const MnAlgebraicVector & Vec() const
determines the relative floating point arithmetic precision.
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
A point of a parabola.
const MnAlgebraicVector & G2() const
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double X() const
Accessor to the x (first) coordinate.
unsigned int size() const
Definition: LAVector.h:198
const MnAlgebraicVector & Gstep() const
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
double Eps2() const
eps2 returns 2*sqrt(eps)
MinimumError keeps the inv.
Definition: MinimumError.h:26
const MnAlgebraicVector & Dirin() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
MinimumState operator()(const MnFcn &, const MinimumState &, const GradientCalculator &, const MnMachinePrecision &) const
Implements a 1-dimensional minimization along a given direction (i.e.
Definition: MnLineSearch.h:47
const Int_t n
Definition: legend1.C:16
interface class for gradient calculators