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