dfp.hpp 1.46 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 39 // See also license.mtl.txt in the distribution. #ifndef ITL_DFP_INCLUDE #define ITL_DFP_INCLUDE #include #include #include namespace itl { /// Update of Hessian matrix for e.g. Quasi-Newton by Davidon, Fletcher and Powell formula struct dfp { /// \f$H_{k+1}=B_{k+1}^{-1}=H_k+\frac{s_k\cdot s_k^T}{y_k^T\cdot s_k}- \frac{H_k\cdot y_k\cdot y_k^T\cdot H_k^T}{y_k^TH_k\cdot y_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); value_type gamma= 1 / dot(y,s), alpha= 1 / dot(y,h); MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality()); MTL_THROW_IF(alpha == 0.0, unexpected_orthogonality()); Matrix A(alpha * y * trans(y)), H2(H - H * A * H + gamma * s * trans(s)); swap(H2, H); // faster than H= H2 }  Praetorius, Simon committed Oct 14, 2016 40 };  Praetorius, Simon committed Oct 13, 2016 41 42 43 44 45 46  } // namespace itl #endif // ITL_DFP_INCLUDE