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
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
22namespace ROOT {
23
24 namespace Minuit2 {
25
26
27
28
29double inner_product(const LAVector&, const LAVector&);
30double similarity(const LAVector&, const LASymMatrix&);
31double 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
37class LASquareMatrix {
38public:
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; }
55private:
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
61LASquareMatrix 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
73LASquareMatrix 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();
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
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define b(i)
Definition: RSha256.hxx:100
#define s0(x)
Definition: RSha256.hxx:90
TRObject operator()(const T1 &t1) const
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
unsigned int Nrow() const
Definition: LASymMatrix.h:239
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
const Int_t n
Definition: legend1.C:16
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;.
double sum_of_elements(const LASymMatrix &)
double similarity(const LAVector &, const LASymMatrix &)
LASquareMatrix OuterProduct(const LAVector &v1, const LAVector &v2)
double inner_product(const LAVector &, const LAVector &)
LASquareMatrix MatrixProduct(const LASymMatrix &m1, const LASquareMatrix &m2)
VSD Structures.
Definition: StringConv.hxx:21
static constexpr double m2
auto * a
Definition: textangle.C:12