Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include <panoc-alm/inner/directions/decl/lbfgs.hpp>
       4             : 
       5             : namespace pa {
       6             : 
       7             : /// Limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm that can
       8             : /// handle updates of the γ parameter.
       9             : /// @ingroup    grp_PANOCDirectionProviders
      10           4 : class SpecializedLBFGS {
      11             :   public:
      12             :     using Params = LBFGSParams;
      13             : 
      14             :     SpecializedLBFGS(Params params) : params(params) {}
      15          12 :     SpecializedLBFGS(Params params, size_t n, size_t history) : params(params) {
      16           4 :         resize(n, history);
      17           4 :     }
      18             : 
      19             :     /// Standard L-BFGS update without changing the step size γ.
      20             :     bool standard_update(crvec xₖ, crvec xₖ₊₁, crvec pₖ,
      21             :                          crvec pₖ₊₁, crvec gradₖ₊₁);
      22             :     /// L-BFGS update when changing the step size γ, recomputing everything.
      23             :     bool full_update(crvec xₖ, crvec xₖ₊₁, crvec pₖ_old_γ,
      24             :                      crvec pₖ₊₁, crvec gradₖ₊₁, const Box &C,
      25             :                      real_t γ);
      26             :     /// Update the inverse Hessian approximation using the new vectors xₖ₊₁
      27             :     /// and pₖ₊₁.
      28             :     bool update(crvec xₖ, crvec xₖ₊₁, crvec pₖ, crvec pₖ₊₁,
      29             :                 crvec gradₖ₊₁, const Box &C, real_t γ);
      30             : 
      31             :     /// Apply the inverse Hessian approximation to the given vector q.
      32             :     template <class Vec>
      33             :     void apply(Vec &&q);
      34             : 
      35             :     /// Initialize with the starting point x₀ and the gradient in that point.
      36             :     void initialize(crvec x₀, crvec grad₀);
      37             : 
      38             :     /// Throw away the approximation and all previous vectors s and y.
      39             :     void reset();
      40             :     /// Re-allocate storage for a problem with a different size. 
      41             :     void resize(size_t n, size_t history);
      42             : 
      43             :     std::string get_name() const { return "SpecializedLBFGS"; }
      44             : 
      45             :     /// Get the size of the s, y, x and g vectors in the buffer.
      46         582 :     size_t n() const { return sto.rows() - 1; }
      47             :     /// Get the number of previous vectors s, y, x and g stored in the buffer.
      48         428 :     size_t history() const { return (sto.cols() - 2) / 4; }
      49             :     /// Get the next index in the circular buffer of previous s, y, x and g
      50             :     /// vectors.
      51          96 :     size_t succ(size_t i) const { return i + 1 < history() ? i + 1 : 0; }
      52             :     /// Get the previous index in the circular buffer of previous s, y, x and g
      53             :     /// vectors.
      54          20 :     size_t pred(size_t i) const { return i == 0 ? history() - 1 : i - 1; }
      55             : 
      56         104 :     auto s(size_t i) { return sto.col(2 * i).topRows(n()); }
      57             :     auto s(size_t i) const { return sto.col(2 * i).topRows(n()); }
      58         124 :     auto y(size_t i) { return sto.col(2 * i + 1).topRows(n()); }
      59             :     auto y(size_t i) const { return sto.col(2 * i + 1).topRows(n()); }
      60         158 :     auto x(size_t i) { return sto.col(2 * history() + 2 * i).topRows(n()); }
      61             :     auto x(size_t i) const {
      62             :         return sto.col(2 * history() + 2 * i).topRows(n());
      63             :     }
      64         152 :     auto g(size_t i) { return sto.col(2 * history() + 2 * i + 1).topRows(n()); }
      65             :     auto g(size_t i) const {
      66             :         return sto.col(2 * history() + 2 * i + 1).topRows(n());
      67             :     }
      68           6 :     auto p() { return sto.col(4 * history()).topRows(n()); }
      69             :     auto p() const { return sto.col(4 * history()).topRows(n()); }
      70           6 :     auto w() { return sto.col(4 * history() + 1).topRows(n()); }
      71             :     auto w() const { return sto.col(4 * history() + 1).topRows(n()); }
      72          32 :     real_t &ρ(size_t i) { return sto.coeffRef(n(), 2 * i); }
      73             :     const real_t &ρ(size_t i) const { return sto.coeff(n(), 2 * i); }
      74             :     real_t &α(size_t i) { return sto.coeffRef(n(), 2 * i + 1); }
      75             :     const real_t &α(size_t i) const { return sto.coeff(n(), 2 * i + 1); }
      76             : 
      77             :   private:
      78             :     using storage_t = Eigen::Matrix<real_t, Eigen::Dynamic, Eigen::Dynamic>;
      79             : 
      80             :     storage_t sto;
      81           4 :     size_t idx   = 0;
      82           4 :     bool full    = false;
      83           4 :     real_t old_γ = 0;
      84             :     Params params;
      85             : };
      86             : 
      87             : } // namespace pa