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#include "Minuit2/MnPrint.h"
19
20#include <cmath>
21
22namespace ROOT {
23
24namespace Minuit2 {
25
27 const MnMachinePrecision &prec) const
28{
29
30 // when the second derivatives are negative perform a line search along Parameter which gives
31 // negative second derivative and magnitude equal to the Gradient step size.
32 // Recalculate the gradients for all the Parameter after the correction and
33 // continue iteration in case the second derivatives are still negative
34 //
35 MnPrint print("NegativeG2LineSearch");
36
37 bool negG2 = HasNegativeG2(st.Gradient(), prec);
38 if (!negG2)
39 return st;
40
41 print.Info("Doing a NegativeG2LineSearch since one of the G2 component is negative");
42
43 unsigned int n = st.Parameters().Vec().size();
44 FunctionGradient dgrad = st.Gradient();
46 bool iterate = false;
47 unsigned int iter = 0;
48 // in case of analytical gradients we don't have step sizes
49 bool hasGStep = !dgrad.IsAnalytical();
50 print.Debug("Gradient ", dgrad.Vec(), "G2",dgrad.G2());
52 // gradient present in the state must have G2 otherwise something is wrong
53 //assert(dgrad.HasG2());
54 if (!dgrad.HasG2()) {
55 print.Error("Input gradient to NG2LS must have G2 already computed");
56 return st;
57 }
58 bool computeHessian = false;
59
60 do {
61 iterate = false;
62 for (unsigned int i = 0; i < n; i++) {
63
64 if (dgrad.G2()(i) <= 0) {
65
66 // check also the gradient (if it is zero ) I can skip the param)
67
68 if (std::fabs(dgrad.Vec()(i)) < prec.Eps() && std::fabs(dgrad.G2()(i)) < prec.Eps())
69 continue;
70 // if(dgrad.G2()(i) < prec.Eps()) {
71 // do line search if second derivative negative
72 MnAlgebraicVector step(n);
73 MnLineSearch lsearch;
74
75 if (dgrad.Vec()(i) < 0)
76 step(i) = (hasGStep) ? dgrad.Gstep()(i) : 1;
77 else
78 step(i) = (hasGStep) ? -dgrad.Gstep()(i) : -1;
79
80 double gdel = step(i) * dgrad.Vec()(i);
81
82 // if using sec derivative information
83 // double g2del = step(i)*step(i) * dgrad.G2()(i);
84
85 print.Debug("Iter", iter, "param", i, pa.Vec()(i), "grad2", dgrad.G2()(i), "grad",
86 dgrad.Vec()(i), "grad step", step(i), " gdel ", gdel);
87
88 MnParabolaPoint pp = lsearch(fcn, pa, step, gdel, prec);
89
90 print.Debug("Line search result", pp.X(), "f(0)", pa.Fval(), "f(1)", pp.Y());
91
92 step *= pp.X();
93 pa = MinimumParameters(pa.Vec() + step, pp.Y());
94
95 dgrad = gc(pa, dgrad);
96 // re-compute also G2 if needed
97 if (!dgrad.HasG2()) {
98 computeHessian = true;
99 print.Debug("Compute Hessian at the new point", pa.Vec());
100 bool ret = gc.Hessian(pa,mat);
101 if (!ret) {
102 print.Error("Cannot compute Hessian");
103 assert(false);
104 return st;
105 }
107 for (unsigned int j = 0; j < n; j++)
108 g2(j) = mat(j,j);
109 dgrad = FunctionGradient(dgrad.Grad(), g2);
110 }
111
112 print.Debug("New result after Line search - iter", iter, "param", i, pa.Vec()(i), "step", step(i), "new grad2",
113 dgrad.G2()(i), "new grad", dgrad.Vec()(i));
114
115 iterate = true;
116 break;
117 }
118 }
119 } while (iter++ < 2 * n && iterate);
120
121 // for the case where we did not compute Hessian
122 if (!computeHessian) {
123 print.Info("Approximate covariance using only G2");
124 for (unsigned int i = 0; i < n; i++) {
125 mat(i, i) = (std::fabs(dgrad.G2()(i)) > prec.Eps2() ? 1. / dgrad.G2()(i) :
126 (dgrad.G2()(i) >= 0) ? 1./prec.Eps2() : -1./prec.Eps2());
127 }
128 }
129
130 MinimumError err(mat, 1.);
131 double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
132
133 if (edm < 0) {
135 }
136
137 return MinimumState(pa, err, dgrad, edm, fcn.NumOfCalls());
138}
139
141{
142 // check if function gradient has any component which is negative
143
144 for (unsigned int i = 0; i < grad.Vec().size(); i++)
145
146 if (grad.G2()(i) <= 0) {
147
148 return true;
149 }
150
151 return false;
152}
153
154} // namespace Minuit2
155
156} // namespace ROOT
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
const MnAlgebraicVector & Gstep() const
const MnAlgebraicVector & Grad() 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:45
unsigned int size() const
Definition: LAVector.h:231
MinimumError keeps the inv.
Definition: MinimumError.h:28
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:27
const MinimumParameters & Parameters() const
Definition: MinimumState.h:57
const FunctionGradient & Gradient() const
Definition: MinimumState.h:62
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:30
unsigned int NumOfCalls() const
Definition: MnFcn.h:39
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.
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
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)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.