ROOT   Reference Guide
Searching...
No Matches
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#include "Minuit2/MnPrint.h"
15
16#include <vector>
17
18namespace ROOT {
19
20namespace Minuit2 {
21
22double inner_product(const LAVector &, const LAVector &);
23double similarity(const LAVector &, const LASymMatrix &);
24double 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
30class LASquareMatrix {
31public:
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
48private:
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
54LASquareMatrix 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
67LASquareMatrix 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
83MinimumError
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
double
#define b(i)
Definition RSha256.hxx:100
#define s0(x)
Definition RSha256.hxx:90
#define a(i)
Definition RSha256.hxx:99
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:45
unsigned int Nrow() const
MinimumError keeps the inv.
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
void Debug(const Ts &... args)
Definition MnPrint.h:138
void Warn(const Ts &... args)
Definition MnPrint.h:126
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)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...