Stan  2.10.0
probability, sampling & optimization
lbfgs_update.hpp
Go to the documentation of this file.
1 #ifndef STAN_OPTIMIZATION_LBFGS_UPDATE_HPP
2 #define STAN_OPTIMIZATION_LBFGS_UPDATE_HPP
3 
4 #include <stan/math/prim/mat/fun/Eigen.hpp>
5 #include <boost/tuple/tuple.hpp>
6 #include <boost/circular_buffer.hpp>
7 #include <vector>
8 
9 namespace stan {
10  namespace optimization {
16  template<typename Scalar = double,
17  int DimAtCompile = Eigen::Dynamic>
18  class LBFGSUpdate {
19  public:
20  typedef Eigen::Matrix<Scalar, DimAtCompile, 1> VectorT;
21  typedef Eigen::Matrix<Scalar, DimAtCompile, DimAtCompile> HessianT;
22  // NOLINTNEXTLINE(build/include_what_you_use)
23  typedef boost::tuple<Scalar, VectorT, VectorT> UpdateT;
24 
25  explicit LBFGSUpdate(size_t L = 5) : _buf(L) {}
26 
32  void set_history_size(size_t L) {
33  _buf.rset_capacity(L);
34  }
35 
47  inline Scalar update(const VectorT &yk, const VectorT &sk,
48  bool reset = false) {
49  Scalar skyk = yk.dot(sk);
50 
51  Scalar B0fact;
52  if (reset) {
53  B0fact = yk.squaredNorm()/skyk;
54  _buf.clear();
55  } else {
56  B0fact = 1.0;
57  }
58 
59  // New updates are pushed to the "back" of the circular buffer
60  Scalar invskyk = 1.0/skyk;
61  _gammak = skyk/yk.squaredNorm();
62  _buf.push_back();
63  _buf.back() = boost::tie(invskyk, yk, sk);
64 
65  return B0fact;
66  }
67 
76  inline void search_direction(VectorT &pk, const VectorT &gk) const {
77  std::vector<Scalar> alphas(_buf.size());
78  typename
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;
83 
84  pk.noalias() = -gk;
85  for (buf_rit = _buf.rbegin(), alpha_rit = alphas.rbegin();
86  buf_rit != _buf.rend();
87  buf_rit++, alpha_rit++) {
88  Scalar alpha;
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));
92 
93  alpha = rhoi * si.dot(pk);
94  pk -= alpha * yi;
95  *alpha_rit = alpha;
96  }
97  pk *= _gammak;
98  for (buf_it = _buf.begin(), alpha_it = alphas.begin();
99  buf_it != _buf.end();
100  buf_it++, alpha_it++) {
101  Scalar beta;
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));
105 
106  beta = rhoi*yi.dot(pk);
107  pk += (*alpha_it - beta)*si;
108  }
109  }
110 
111  protected:
112  boost::circular_buffer<UpdateT> _buf;
113  Scalar _gammak;
114  };
115  }
116 }
117 
118 #endif
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.

     [ Stan Home Page ] © 2011–2016, Stan Development Team.