1 #ifndef STAN_OPTIMIZATION_LBFGS_UPDATE_HPP
2 #define STAN_OPTIMIZATION_LBFGS_UPDATE_HPP
4 #include <stan/math/prim/mat/fun/Eigen.hpp>
5 #include <boost/tuple/tuple.hpp>
6 #include <boost/circular_buffer.hpp>
10 namespace optimization {
16 template<
typename Scalar = double,
17 int DimAtCompile = Eigen::Dynamic>
20 typedef Eigen::Matrix<Scalar, DimAtCompile, 1>
VectorT;
21 typedef Eigen::Matrix<Scalar, DimAtCompile, DimAtCompile>
HessianT;
23 typedef boost::tuple<Scalar, VectorT, VectorT>
UpdateT;
33 _buf.rset_capacity(L);
47 inline Scalar
update(
const VectorT &yk,
const VectorT &sk,
49 Scalar skyk = yk.dot(sk);
53 B0fact = yk.squaredNorm()/skyk;
60 Scalar invskyk = 1.0/skyk;
61 _gammak = skyk/yk.squaredNorm();
63 _buf.back() = boost::tie(invskyk, yk, sk);
77 std::vector<Scalar> alphas(
_buf.size());
79 boost::circular_buffer<UpdateT>::const_reverse_iterator buf_rit;
80 typename boost::circular_buffer<UpdateT>::const_iterator buf_it;
81 typename std::vector<Scalar>::const_iterator alpha_it;
82 typename std::vector<Scalar>::reverse_iterator alpha_rit;
85 for (buf_rit =
_buf.rbegin(), alpha_rit = alphas.rbegin();
86 buf_rit !=
_buf.rend();
87 buf_rit++, alpha_rit++) {
89 const Scalar &rhoi(boost::get<0>(*buf_rit));
90 const VectorT &yi(boost::get<1>(*buf_rit));
91 const VectorT &si(boost::get<2>(*buf_rit));
93 alpha = rhoi * si.dot(pk);
98 for (buf_it =
_buf.begin(), alpha_it = alphas.begin();
100 buf_it++, alpha_it++) {
102 const Scalar &rhoi(boost::get<0>(*buf_it));
103 const VectorT &yi(boost::get<1>(*buf_it));
104 const VectorT &si(boost::get<2>(*buf_it));
106 beta = rhoi*yi.dot(pk);
107 pk += (*alpha_it - beta)*si;
112 boost::circular_buffer<UpdateT>
_buf;
void search_direction(VectorT &pk, const VectorT &gk) const
Compute the search direction based on the current (inverse) Hessian approximation and given gradient...
Eigen::Matrix< Scalar, DimAtCompile, 1 > VectorT
Probability, optimization and sampling library.
Scalar update(const VectorT &yk, const VectorT &sk, bool reset=false)
Add a new set of update vectors to the history.
Eigen::Matrix< Scalar, DimAtCompile, DimAtCompile > HessianT
boost::circular_buffer< UpdateT > _buf
void set_history_size(size_t L)
Set the number of inverse Hessian updates to keep.
boost::tuple< Scalar, VectorT, VectorT > UpdateT
Implement a limited memory version of the BFGS update.