Logo ROOT  
Reference Guide
BFGSErrorUpdator.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/MinimumState.h"
12 #include "Minuit2/LaSum.h"
13 #include "Minuit2/LaProd.h"
14 
15 #include <vector>
16 
17 //#define DEBUG
18 
19 #if defined(DEBUG) || defined(WARNINGMSG)
20 #include "Minuit2/MnPrint.h"
21 #endif
22 
23 
24 namespace ROOT {
25 
26  namespace Minuit2 {
27 
28 
29 
30 
31 double inner_product(const LAVector&, const LAVector&);
32 double similarity(const LAVector&, const LASymMatrix&);
33 double sum_of_elements(const LASymMatrix&);
34 
35 // define here a square matrix that it is needed for computingthe BFGS update
36 // define just the class, no need for defining operatipons as dane for the Symmetric matrices
37 // since the square matrix will be converted in a symmetric one afterwards
38 
39 class LASquareMatrix {
40 public:
41  LASquareMatrix(unsigned int n) :
42  fNRow(n),
43  fData(std::vector<double> (n*n) )
44  {}
45 
46  double operator()(unsigned int row, unsigned int col) const {
47  assert(row<fNRow && col < fNRow);
48  return fData[col+row*fNRow];
49  }
50 
51  double& operator()(unsigned int row, unsigned int col) {
52  assert(row<fNRow && col < fNRow);
53  return fData[col+row*fNRow];
54  }
55 
56  unsigned int Nrow() const { return fNRow; }
57 private:
58  unsigned int fNRow;
59  std::vector<double> fData;
60 };
61 
62 // compute outer product of two vector of same size to return a square matrix
63 LASquareMatrix OuterProduct(const LAVector& v1, const LAVector& v2) {
64  assert(v1.size() == v2.size() );
65  LASquareMatrix a(v1.size() );
66  for (unsigned int i = 0; i < v1.size() ; ++i) {
67  for (unsigned int j = 0; j < v2.size() ; ++j) {
68  a(i,j) = v1[i]*v2[j];
69  }
70  }
71  return a;
72 }
73 // compute product of symmetric matrix with square matrix
74 
75 LASquareMatrix MatrixProduct(const LASymMatrix& m1, const LASquareMatrix& m2) {
76  unsigned int n = m1.Nrow();
77  assert(n == m2.Nrow() );
78  LASquareMatrix a( n );
79  for (unsigned int i = 0; i < n ; ++i) {
80  for (unsigned int j = 0; j < n ; ++j) {
81  a(i,j) = 0;
82  for (unsigned int k = 0; k < n ; ++k) {
83  a(i,j) += m1(i,k)*m2(k,j);
84  }
85  }
86  }
87  return a;
88 }
89 
90 
91 
92 MinimumError BFGSErrorUpdator::Update(const MinimumState& s0,
93  const MinimumParameters& p1,
94  const FunctionGradient& g1) const {
95 
96  // update of the covarianze matrix (BFGS formula, see Tutorial, par. 4.8 pag 26)
97  // in case of delgam > gvg (PHI > 1) use rank one formula
98  // see par 4.10 pag 30
99 
100  const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
101  MnAlgebraicVector dx = p1.Vec() - s0.Vec();
102  MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
103 
104  double delgam = inner_product(dx, dg); // this is s^T y using wikipedia conventions
105  double gvg = similarity(dg, v0); // this is y^T B^-1 y
106 
107 
108 #ifdef DEBUG
109  std::cout << "dx = " << dx << std::endl;
110  std::cout << "dg = " << dg << std::endl;
111  std::cout<<"delgam= "<<delgam<<" gvg= "<<gvg<<std::endl;
112 #endif
113 
114  if (delgam == 0 ) {
115 #ifdef WARNINGMSG
116  MN_INFO_MSG("BFGSErrorUpdator: delgam = 0 : cannot update - return same matrix ");
117 #endif
118  return s0.Error();
119  }
120 #ifdef WARNINGMSG
121  if (delgam < 0) MN_INFO_MSG("BFGSErrorUpdator: delgam < 0 : first derivatives increasing along search line");
122 #endif
123 
124  if (gvg <= 0 ) {
125  // since v0 is pos def this gvg can be only = 0 if dg = 0 - should never be here
126 #ifdef WARNINGMSG
127  MN_INFO_MSG("BFGSErrorUpdator: gvg <= 0 ");
128 #endif
129  //return s0.Error();
130  }
131 
132 
133  // compute update formula for BFGS
134  // see wikipedia https://en.wikipedia.org/wiki/Broyden–Fletcher–Goldfarb–Shanno_algorithm
135  // need to compute outer product dg . dx and it is not symmetric anymore
136 
137  LASquareMatrix a = OuterProduct(dg,dx);
138  LASquareMatrix b = MatrixProduct(v0, a);
139 
140  unsigned int n = v0.Nrow();
142  for (unsigned int i = 0; i < n; ++i) {
143  for (unsigned int j = i; j < n; ++j) {
144  v2(i,j) = (b(i,j) + b(j,i))/(delgam);
145  }
146  }
147 
148  MnAlgebraicSymMatrix vUpd = ( delgam + gvg) * Outer_product(dx) / (delgam * delgam);
149  vUpd -= v2;
150 
151 
152  double sum_upd = sum_of_elements(vUpd);
153  vUpd += v0;
154 
155  double dcov = 0.5*(s0.Error().Dcovar() + sum_upd/sum_of_elements(vUpd));
156 
157 #ifdef DEBUG
158  std::cout << "BFGSErrorUpdator - dcov is " << dcov << std::endl;
159 #endif
160 
161  return MinimumError(vUpd, dcov);
162 }
163 
164 
165  } // namespace Minuit2
166 
167 } // namespace ROOT
n
const Int_t n
Definition: legend1.C:16
v0
@ v0
Definition: rootcling_impl.cxx:3636
ROOT::Minuit2::LAVector
Definition: LAVector.h:41
BFGSErrorUpdator.h
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:59
ROOT::Minuit2::MnAlgebraicSymMatrix
LASymMatrix MnAlgebraicSymMatrix
Definition: MnMatrix.h:36
operator()
TRObject operator()(const T1 &t1) const
Definition: TRFunctionImport__oprtr.h:14
b
#define b(i)
Definition: RSha256.hxx:118
LaProd.h
ROOT::Minuit2::MatrixProduct
LASquareMatrix MatrixProduct(const LASymMatrix &m1, const LASquareMatrix &m2)
Definition: BFGSErrorUpdator.cxx:83
s0
#define s0(x)
Definition: RSha256.hxx:108
TGeant4Unit::m2
static constexpr double m2
Definition: TGeant4SystemOfUnits.h:129
ROOT::Minuit2::inner_product
double inner_product(const LAVector &, const LAVector &)
Definition: LaInnerProduct.cxx:27
a
auto * a
Definition: textangle.C:12
ROOT::Minuit2::OuterProduct
LASquareMatrix OuterProduct(const LAVector &v1, const LAVector &v2)
Definition: BFGSErrorUpdator.cxx:71
double
double
Definition: Converters.cxx:921
ROOT::Minuit2::sum_of_elements
double sum_of_elements(const LASymMatrix &)
Definition: LaSumOfElements.cxx:34
MN_INFO_MSG
#define MN_INFO_MSG(str)
Definition: MnPrint.h:111
v1
@ v1
Definition: rootcling_impl.cxx:3637
ROOT::Minuit2::similarity
double similarity(const LAVector &, const LASymMatrix &)
Definition: LaVtMVSimilarity.cxx:29
v2
@ v2
Definition: rootcling_impl.cxx:3638
ROOT::Minuit2::BFGSErrorUpdator::Update
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const
Definition: BFGSErrorUpdator.cxx:100
ROOT::Minuit2::Outer_product
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;.
Definition: LaOuterProduct.h:36
LaSum.h
MnPrint.h
ROOT
VSD Structures.
Definition: StringConv.hxx:21
MinimumState.h