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