broyden.hpp 1.34 KB
 Praetorius, Simon committed Oct 13, 2016 1 // Software License for MTL  Praetorius, Simon committed Oct 14, 2016 2 3 // // Copyright (c) 2007 The Trustees of Indiana University.  Praetorius, Simon committed Oct 13, 2016 4 // 2008 Dresden University of Technology and the Trustees of Indiana University.  Praetorius, Simon committed Oct 14, 2016 5 // 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.  Praetorius, Simon committed Oct 13, 2016 6 7 // All rights reserved. // Authors: Peter Gottschling and Andrew Lumsdaine  Praetorius, Simon committed Oct 14, 2016 8 //  Praetorius, Simon committed Oct 13, 2016 9 // This file is part of the Matrix Template Library  Praetorius, Simon committed Oct 14, 2016 10 //  Praetorius, Simon committed Oct 13, 2016 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 // See also license.mtl.txt in the distribution. #ifndef ITL_BROYDEN_INCLUDE #define ITL_BROYDEN_INCLUDE #include #include #include namespace itl { /// Update of Hessian matrix for e.g. Quasi-Newton by Broyden formula struct broyden { /// \f$H_{k+1}=B_{k+1}^{-1}=H_k+\frac{(s_k-H_k\cdot y_k)\cdot y_k^T\cdot H_k}{y_k^T\cdot H_k\cdot s_k} \f$ template void operator() (Matrix& H, const Vector& y, const Vector& s) { typedef typename mtl::Collection::value_type value_type; assert(num_rows(H) == num_cols(H)); Vector h(H * y), d(s - h); value_type gamma= 1 / dot(y, h); MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality()); Matrix A(gamma * d * trans(y)), H2(H + A * H); swap(H2, H); // faster than H= H2 }  Praetorius, Simon committed Oct 14, 2016 39 };  Praetorius, Simon committed Oct 13, 2016 40 41 42 43 44 45  } // namespace itl #endif // ITL_BROYDEN_INCLUDE