ROOT   Reference Guide
DavidonErrorUpdator.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
12#include "Minuit2/LaSum.h"
13#include "Minuit2/LaProd.h"
14
15//#define DEBUG
16
17#if defined(DEBUG) || defined(WARNINGMSG)
18#include "Minuit2/MnPrint.h"
19#endif
20
21
22
23namespace ROOT {
24
25 namespace Minuit2 {
26
27
28double inner_product(const LAVector&, const LAVector&);
29double similarity(const LAVector&, const LASymMatrix&);
30double sum_of_elements(const LASymMatrix&);
31
33 const MinimumParameters& p1,
34 const FunctionGradient& g1) const {
35
36 // update of the covarianze matrix (Davidon formula, see Tutorial, par. 4.8 pag 26)
37 // in case of delgam > gvg (PHI > 1) use rank one formula
38 // see par 4.10 pag 30
39
40 const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
41 MnAlgebraicVector dx = p1.Vec() - s0.Vec();
42 MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
43
44 double delgam = inner_product(dx, dg);
45 double gvg = similarity(dg, v0);
46
47
48#ifdef DEBUG
49 std::cout << "dx = " << dx << std::endl;
50 std::cout << "dg = " << dg << std::endl;
51 std::cout<<"delgam= "<<delgam<<" gvg= "<<gvg<<std::endl;
52#endif
53
54 if (delgam == 0 ) {
55#ifdef WARNINGMSG
56 MN_INFO_MSG("DavidonErrorUpdator: delgam = 0 : cannot update - return same matrix ");
57#endif
58 return s0.Error();
59 }
60#ifdef WARNINGMSG
61 if (delgam < 0) MN_INFO_MSG("DavidonErrorUpdator: delgam < 0 : first derivatives increasing along search line");
62#endif
63
64 if (gvg <= 0 ) {
65 // since v0 is pos def this gvg can be only = 0 if dg = 0 - should never be here
66#ifdef WARNINGMSG
67 MN_INFO_MSG("DavidonErrorUpdator: gvg <= 0 : cannot update - return same matrix ");
68#endif
69 return s0.Error();
70 }
71
72
73 MnAlgebraicVector vg = v0*dg;
74
75 MnAlgebraicSymMatrix vUpd = Outer_product(dx)/delgam - Outer_product(vg)/gvg;
76
77 if(delgam > gvg) {
78 // use rank 1 formula
79 vUpd += gvg*Outer_product(MnAlgebraicVector(dx/delgam - vg/gvg));
80 }
81
82 double sum_upd = sum_of_elements(vUpd);
83 vUpd += v0;
84
85 double dcov = 0.5*(s0.Error().Dcovar() + sum_upd/sum_of_elements(vUpd));
86
87 return MinimumError(vUpd, dcov);
88}
89
90/*
91MinimumError DavidonErrorUpdator::Update(const MinimumState& s0,
92 const MinimumParameters& p1,
93 const FunctionGradient& g1) const {
94
95 const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
96 MnAlgebraicVector dx = p1.Vec() - s0.Vec();
97 MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
98
99 double delgam = inner_product(dx, dg);
100 double gvg = similarity(dg, v0);
101
102// std::cout<<"delgam= "<<delgam<<" gvg= "<<gvg<<std::endl;
103 MnAlgebraicVector vg = v0*dg;
104// MnAlgebraicSymMatrix vUpd(v0.Nrow());
105
106// MnAlgebraicSymMatrix dd = ( 1./delgam )*outer_product(dx);
107// dd *= ( 1./delgam );
108// MnAlgebraicSymMatrix VggV = ( 1./gvg )*outer_product(vg);
109// VggV *= ( 1./gvg );
110// vUpd = dd - VggV;
111// MnAlgebraicSymMatrix vUpd = ( 1./delgam )*outer_product(dx) - ( 1./gvg )*outer_product(vg);
112 MnAlgebraicSymMatrix vUpd = Outer_product(dx)/delgam - Outer_product(vg)/gvg;
113
114 if(delgam > gvg) {
115// dx *= ( 1./delgam );
116// vg *= ( 1./gvg );
117// MnAlgebraicVector flnu = dx - vg;
118// MnAlgebraicSymMatrix tmp = Outer_product(flnu);
119// tmp *= gvg;
120// vUpd = vUpd + tmp;
121 vUpd += gvg*outer_product(dx/delgam - vg/gvg);
122 }
123
124//
125// MnAlgebraicSymMatrix dd = Outer_product(dx);
126// dd *= ( 1./delgam );
127// MnAlgebraicSymMatrix VggV = Outer_product(vg);
128// VggV *= ( 1./gvg );
129// vUpd = dd - VggV;
130//
131//
132// double phi = delgam/(delgam - gvg);
133
134// MnAlgebraicSymMatrix vUpd(v0.Nrow());
135// if(phi < 0) {
136// // rank-2 Update
137// MnAlgebraicSymMatrix dd = Outer_product(dx);
138// dd *= ( 1./delgam );
139// MnAlgebraicSymMatrix VggV = Outer_product(vg);
140// VggV *= ( 1./gvg );
141// vUpd = dd - VggV;
142// }
143// if(phi > 1) {
144// // rank-1 Update
145// MnAlgebraicVector tmp = dx - vg;
146// vUpd = Outer_product(tmp);
147// vUpd *= ( 1./(delgam - gvg) );
148// }
149//
150
151//
152// if(delgam > gvg) {
153// // rank-1 Update
154// MnAlgebraicVector tmp = dx - vg;
155// vUpd = Outer_product(tmp);
156// vUpd *= ( 1./(delgam - gvg) );
157// } else {
158// // rank-2 Update
159// MnAlgebraicSymMatrix dd = Outer_product(dx);
160// dd *= ( 1./delgam );
161// MnAlgebraicSymMatrix VggV = Outer_productn(vg);
162// VggV *= ( 1./gvg );
163// vUpd = dd - VggV;
164// }
165//
166
167 double sum_upd = sum_of_elements(vUpd);
168 vUpd += v0;
169
170// MnAlgebraicSymMatrix V1 = v0 + vUpd;
171
172 double dcov =
173 0.5*(s0.Error().Dcovar() + sum_upd/sum_of_elements(vUpd));
174
175 return MinimumError(vUpd, dcov);
176}
177*/
178
179 } // namespace Minuit2
180
181} // namespace ROOT
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define s0(x)
Definition: RSha256.hxx:90
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const
const MnAlgebraicVector & Vec() const
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
MinimumError keeps the inv.
Definition: MinimumError.h:26
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, double >, double >, double > Outer_product(const ABObj< vec, LAVector, double > &obj)
LAPACK Algebra function specialize the Outer_product function for LAVector;.
LAVector MnAlgebraicVector
Definition: MnMatrix.h:42
double sum_of_elements(const LASymMatrix &)
double similarity(const LAVector &, const LASymMatrix &)
double inner_product(const LAVector &, const LAVector &)
VSD Structures.
Definition: StringConv.hxx:21