Reference documentation for deal.II version Git f40be01994 2020-04-09 07:13:12 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
solver_bfgs.h
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 //---------------------------------------------------------------
15 
16 #ifndef dealii_solver_bfgs_h
17 #define dealii_solver_bfgs_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/lac/solver.h>
22 
23 #include <deal.II/numerics/history.h>
24 
25 #include <deal.II/optimization/line_minimization.h>
26 
27 DEAL_II_NAMESPACE_OPEN
28 
57 template <typename VectorType>
58 class SolverBFGS : public SolverBase<VectorType>
59 {
60 public:
64  typedef typename VectorType::value_type Number;
65 
66 
71  {
75  explicit AdditionalData(const unsigned int max_history_size = 5,
76  const bool debug_output = false);
77 
81  unsigned int max_history_size;
82 
87  };
88 
89 
93  explicit SolverBFGS(SolverControl & residual_control,
94  const AdditionalData &data = AdditionalData());
95 
109  void
110  solve(
111  const std::function<Number(const VectorType &x, VectorType &g)> &compute,
112  VectorType & x);
113 
122  boost::signals2::connection
124  const std::function<
125  Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
126  &slot);
127 
154  boost::signals2::connection
156  const std::function<void(VectorType & g,
158  const FiniteSizeHistory<VectorType> &y)> &slot);
159 
160 
161 protected:
166 
170  boost::signals2::signal<
171  Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
173 
177  boost::signals2::signal<void(VectorType & g,
181 };
182 
183 
184 // ------------------- inline and template functions ----------------
185 #ifndef DOXYGEN
186 
187 template <typename VectorType>
189  const unsigned int max_history_size_,
190  const bool debug_output_)
191  : max_history_size(max_history_size_)
192  , debug_output(debug_output_)
193 {}
194 
195 
196 
197 template <typename VectorType>
199  const AdditionalData &data)
200  : SolverBase<VectorType>(solver_control)
201  , additional_data(data)
202 {}
203 
204 
205 
206 template <class VectorType>
207 boost::signals2::connection
209  const std::function<
210  Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> &slot)
211 {
212  Assert(line_search_signal.empty(),
213  ExcMessage("One should not attach more than one line search signal."));
214  return line_search_signal.connect(slot);
215 }
216 
217 
218 
219 template <class VectorType>
220 boost::signals2::connection
222  const std::function<void(VectorType & g,
224  const FiniteSizeHistory<VectorType> &y)> &slot)
225 {
227  ExcMessage(
228  "One should not attach more than one preconditioner signal."));
229  return preconditioner_signal.connect(slot);
230 }
231 
232 
233 
234 template <typename VectorType>
235 void
237  const std::function<typename VectorType::value_type(const VectorType &x,
238  VectorType &f)> &compute,
239  VectorType & x)
240 {
241  // Also see scipy Fortran implementation
242  // https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f
243  // and Octave-optim implementation:
244  // https://sourceforge.net/p/octave/optim/ci/default/tree/src/__bfgsmin.cc
245  LogStream::Prefix prefix("BFGS");
246 
247  // default line search:
248  bool first_step = true;
249  Number f_prev = 0.;
250  // provide default line search if no signal was attached
251  VectorType x0;
252  if (line_search_signal.empty())
253  {
254  x0.reinit(x);
255  const auto default_line_min =
256  [&](Number &f, VectorType &x, VectorType &g, const VectorType &p) {
257  const Number f0 = f;
258  const Number g0 = g * p;
259  Assert(g0 < 0,
260  ExcMessage(
261  "Function does not decrease along the current direction"));
262 
263  // save current solution value (to be used in line_search):
264  x0 = x;
265 
266  // see scipy implementation
267  // https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.line_search.html#scipy.optimize.line_search
268  // and Eq. 2.6.8 in Fletcher 2013, Practical methods of optimization
269  Number df = f_prev - f;
270  Assert(first_step || df >= 0.,
271  ExcMessage("Function value is not decreasing"));
272  df = std::max(df, 100. * std::numeric_limits<Number>::epsilon());
273  // guess a reasonable first step:
274  const Number a1 =
275  (first_step ? 1. : std::min(1., -1.01 * 2. * df / g0));
276  Assert(a1 > 0., ExcInternalError());
277  f_prev = f;
278 
279  // 1D line-search function
280  const auto line_func =
281  [&](const Number &x_line) -> std::pair<Number, Number> {
282  x = x0;
283  x.add(x_line, p);
284  f = compute(x, g);
285  const Number g_line = g * p;
286  return std::make_pair(f, g_line);
287  };
288 
289  // loose line search:
290  const auto res = LineMinimization::line_search<Number>(
291  line_func,
292  f0,
293  g0,
294  LineMinimization::poly_fit<Number>,
295  a1,
296  0.9,
297  0.001);
298 
299  if (first_step)
300  first_step = false;
301 
302  return res.first;
303  };
304  this->connect_line_search_slot(default_line_min);
305  }
306 
307  // FIXME: Octave has convergence in terms of:
308  // function change tolerance, default 1e-12
309  // parameter change tolerance, default 1e-6
310  // gradient tolerance, default 1e-5
311  // SolverBase and/or SolverControl need extension
312 
313  VectorType g(x), p(x), y_k(x), s_k(x);
314 
315  std::vector<Number> c1;
316  c1.reserve(additional_data.max_history_size);
317 
318  // limited history
322 
323  unsigned int m = 0;
324  Number f;
325 
327  unsigned int k = 0;
328 
329  f = compute(x, g);
330 
331  conv = this->iteration_status(k, g.l2_norm(), x);
332  if (conv != SolverControl::iterate)
333  return;
334 
335  while (conv == SolverControl::iterate)
336  {
338  deallog << "Iteration " << k << " history " << m << std::endl
339  << "f=" << f << std::endl;
340 
341  // 1. Two loop recursion to calculate p = - H*g
342  c1.resize(m);
343  p = g;
344  // first loop:
345  for (unsigned int i = 0; i < m; ++i)
346  {
347  c1[i] = rho[i] * (s[i] * p);
348  p.add(-c1[i], y[i]);
349  }
350  // H0
351  if (!preconditioner_signal.empty())
352  preconditioner_signal(p, s, y);
353 
354  // second loop:
355  for (int i = m - 1; i >= 0; --i)
356  {
357  Assert(i >= 0, ExcInternalError());
358  const Number c2 = rho[i] * (y[i] * p);
359  p.add(c1[i] - c2, s[i]);
360  }
361  p *= -1.;
362 
363  // 2. Line search
364  s_k = x;
365  y_k = g;
366  const Number alpha = line_search_signal(f, x, g, p)
367  .get(); // <-- signals return boost::optional
368  s_k.sadd(-1, 1, x);
369  y_k.sadd(-1, 1, g);
370 
372  deallog << "Line search a=" << alpha << " f=" << f << std::endl;
373 
374  // 3. Check convergence
375  k++;
376  const Number g_l2 = g.l2_norm();
377  conv = this->iteration_status(k, g_l2, x);
378  if (conv != SolverControl::iterate)
379  break;
380 
381  // 4. Store s, y, rho
382  const Number curvature = s_k * y_k;
384  deallog << "Curvature " << curvature << std::endl;
385 
386  if (curvature > 0. && additional_data.max_history_size > 0)
387  {
388  s.add(s_k);
389  y.add(y_k);
390  rho.add(1. / curvature);
391  m = s.size();
392 
393  Assert(y.size() == m, ExcInternalError());
394  Assert(rho.size() == m, ExcInternalError());
395  }
396 
398  }
399 
400  // In the case of failure: throw exception.
402  SolverControl::NoConvergence(k, g.l2_norm()));
403 }
404 
405 #endif
406 
407 DEAL_II_NAMESPACE_CLOSE
408 
409 #endif
boost::signals2::signal< Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> line_search_signal
Definition: solver_bfgs.h:172
boost::signals2::connection connect_preconditioner_slot(const std::function< void(VectorType &g, const FiniteSizeHistory< VectorType > &s, const FiniteSizeHistory< VectorType > &y)> &slot)
Continue iteration.
std::size_t size() const
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:463
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
unsigned int max_history_size
Definition: solver_bfgs.h:81
AdditionalData(const unsigned int max_history_size=5, const bool debug_output=false)
static ::ExceptionBase & ExcMessage(std::string arg1)
Stop iteration, goal reached.
const AdditionalData additional_data
Definition: solver_bfgs.h:165
#define Assert(cond, exc)
Definition: exceptions.h:1419
VectorType::value_type Number
Definition: solver_bfgs.h:64
boost::signals2::signal< void(VectorType &g, const FiniteSizeHistory< VectorType > &s, const FiniteSizeHistory< VectorType > &y)> preconditioner_signal
Definition: solver_bfgs.h:180
SolverBFGS(SolverControl &residual_control, const AdditionalData &data=AdditionalData())
void add(const T &element)
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)
static ::ExceptionBase & ExcInternalError()