solver_bfgs.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2018 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_solver_bfgs_h
#define dealii_solver_bfgs_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/lac/solver.h>
 
#include <deal.II/numerics/history.h>
 
#include <deal.II/optimization/line_minimization.h>
 
#include <limits>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename VectorType>
class SolverBFGS : public SolverBase<VectorType>
{
public:
  using Number = typename VectorType::value_type;
 
 
  struct AdditionalData
  {
    explicit AdditionalData(const unsigned int max_history_size = 5,
                            const bool         debug_output     = false);
 
    unsigned int max_history_size;
 
    bool debug_output;
  };
 
 
  explicit SolverBFGS(SolverControl        &residual_control,
                      const AdditionalData &data = AdditionalData());
 
  void
  solve(
    const std::function<Number(const VectorType &x, VectorType &g)> &compute,
    VectorType                                                      &x);
 
  boost::signals2::connection
  connect_line_search_slot(
    const std::function<
      Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
      &slot);
 
  boost::signals2::connection
  connect_preconditioner_slot(
    const std::function<void(VectorType                          &g,
                             const FiniteSizeHistory<VectorType> &s,
                             const FiniteSizeHistory<VectorType> &y)> &slot);
 
 
protected:
  const AdditionalData additional_data;
 
  boost::signals2::signal<
    Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
    line_search_signal;
 
  boost::signals2::signal<void(VectorType                          &g,
                               const FiniteSizeHistory<VectorType> &s,
                               const FiniteSizeHistory<VectorType> &y)>
    preconditioner_signal;
};
 
 
// -------------------  inline and template functions ----------------
#ifndef DOXYGEN
 
template <typename VectorType>
SolverBFGS<VectorType>::AdditionalData::AdditionalData(
  const unsigned int max_history_size_,
  const bool         debug_output_)
  : max_history_size(max_history_size_)
  , debug_output(debug_output_)
{}
 
 
 
template <typename VectorType>
SolverBFGS<VectorType>::SolverBFGS(SolverControl        &solver_control,
                                   const AdditionalData &data)
  : SolverBase<VectorType>(solver_control)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
boost::signals2::connection
SolverBFGS<VectorType>::connect_line_search_slot(
  const std::function<
    Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> &slot)
{
  Assert(line_search_signal.empty(),
         ExcMessage("One should not attach more than one line search signal."));
  return line_search_signal.connect(slot);
}
 
 
 
template <typename VectorType>
boost::signals2::connection
SolverBFGS<VectorType>::connect_preconditioner_slot(
  const std::function<void(VectorType                          &g,
                           const FiniteSizeHistory<VectorType> &s,
                           const FiniteSizeHistory<VectorType> &y)> &slot)
{
  Assert(preconditioner_signal.empty(),
         ExcMessage(
           "One should not attach more than one preconditioner signal."));
  return preconditioner_signal.connect(slot);
}
 
 
 
template <typename VectorType>
void
SolverBFGS<VectorType>::solve(
  const std::function<typename VectorType::value_type(const VectorType &x,
                                                      VectorType &f)> &compute,
  VectorType                                                          &x)
{
  // Also see scipy Fortran implementation
  // https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f
  // and Octave-optim implementation:
  // https://sourceforge.net/p/octave/optim/ci/default/tree/src/__bfgsmin.cc
  LogStream::Prefix prefix("BFGS");
 
  // default line search:
  bool   first_step = true;
  Number f_prev     = 0.;
  // provide default line search if no signal was attached
  VectorType x0;
  if (line_search_signal.empty())
    {
      x0.reinit(x);
      const auto default_line_min =
        [&](Number &f, VectorType &x, VectorType &g, const VectorType &p) {
          const Number f0 = f;
          const Number g0 = g * p;
          Assert(g0 < 0,
                 ExcMessage(
                   "Function does not decrease along the current direction"));
 
          // save current solution value (to be used in line_search):
          x0 = x;
 
          // see scipy implementation
          // https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.line_search.html#scipy.optimize.line_search
          // and Eq. 2.6.8 in Fletcher 2013, Practical methods of optimization
          Number df = f_prev - f;
          Assert(first_step || df >= 0.,
                 ExcMessage("Function value is not decreasing"));
          df = std::max(df, 100. * std::numeric_limits<Number>::epsilon());
          // guess a reasonable first step:
          const Number a1 =
            (first_step ? 1. : std::min(1., -1.01 * 2. * df / g0));
          Assert(a1 > 0., ExcInternalError());
          f_prev = f;
 
          // 1d line-search function
          const auto line_func =
            [&](const Number &x_line) -> std::pair<Number, Number> {
            x = x0;
            x.add(x_line, p);
            f                   = compute(x, g);
            const Number g_line = g * p;
            return std::make_pair(f, g_line);
          };
 
          // loose line search:
          const auto res = LineMinimization::line_search<Number>(
            line_func,
            f0,
            g0,
            LineMinimization::poly_fit<Number>,
            a1,
            0.9,
            0.001);
 
          if (first_step)
            first_step = false;
 
          return res.first;
        };
      this->connect_line_search_slot(default_line_min);
    }
 
  // FIXME: Octave has convergence in terms of:
  // function change tolerance, default 1e-12
  // parameter change tolerance, default 1e-6
  // gradient tolerance, default 1e-5
  // SolverBase and/or SolverControl need extension
 
  VectorType g(x), p(x), y_k(x), s_k(x);
 
  std::vector<Number> c1;
  c1.reserve(additional_data.max_history_size);
 
  // limited history
  FiniteSizeHistory<VectorType> y(additional_data.max_history_size);
  FiniteSizeHistory<VectorType> s(additional_data.max_history_size);
  FiniteSizeHistory<Number>     rho(additional_data.max_history_size);
 
  unsigned int m = 0;
  Number       f;
 
  SolverControl::State conv = SolverControl::iterate;
  unsigned int         k    = 0;
 
  f = compute(x, g);
 
  conv = this->iteration_status(k, g.l2_norm(), x);
  if (conv != SolverControl::iterate)
    return;
 
  while (conv == SolverControl::iterate)
    {
      if (additional_data.debug_output)
        deallog << "Iteration " << k << " history " << m << std::endl
                << "f=" << f << std::endl;
 
      // 1. Two loop recursion to calculate p = - H*g
      c1.resize(m);
      p = g;
      // first loop:
      for (unsigned int i = 0; i < m; ++i)
        {
          c1[i] = rho[i] * (s[i] * p);
          p.add(-c1[i], y[i]);
        }
      // H0
      if (!preconditioner_signal.empty())
        preconditioner_signal(p, s, y);
 
      // second loop:
      for (int i = m - 1; i >= 0; --i)
        {
          Assert(i >= 0, ExcInternalError());
          const Number c2 = rho[i] * (y[i] * p);
          p.add(c1[i] - c2, s[i]);
        }
      p *= -1.;
 
      // 2. Line search
      s_k                = x;
      y_k                = g;
      const Number alpha = line_search_signal(f, x, g, p)
                             .get(); // <-- signals return boost::optional
      s_k.sadd(-1, 1, x);
      y_k.sadd(-1, 1, g);
 
      if (additional_data.debug_output)
        deallog << "Line search a=" << alpha << " f=" << f << std::endl;
 
      // 3. Check convergence
      ++k;
      const Number g_l2 = g.l2_norm();
      conv              = this->iteration_status(k, g_l2, x);
      if (conv != SolverControl::iterate)
        break;
 
      // 4. Store s, y, rho
      const Number curvature = s_k * y_k;
      if (additional_data.debug_output)
        deallog << "Curvature " << curvature << std::endl;
 
      if (curvature > 0. && additional_data.max_history_size > 0)
        {
          s.add(s_k);
          y.add(y_k);
          rho.add(1. / curvature);
          m = s.size();
 
          Assert(y.size() == m, ExcInternalError());
          Assert(rho.size() == m, ExcInternalError());
        }
 
      Assert(m <= additional_data.max_history_size, ExcInternalError());
    }
 
  // In the case of failure: throw exception.
  AssertThrow(conv == SolverControl::success,
              SolverControl::NoConvergence(k, g.l2_norm()));
}
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif