Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.1.1
\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
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/lac/solver.h>
20 
21 #include <deal.II/numerics/history.h>
22 
23 #include <deal.II/optimization/line_minimization.h>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
55 template <typename VectorType>
56 class SolverBFGS : public SolverBase<VectorType>
57 {
58 public:
62  typedef typename VectorType::value_type Number;
63 
64 
69  {
73  explicit AdditionalData(const unsigned int max_history_size = 5,
74  const bool debug_output = false);
75 
79  unsigned int max_history_size;
80 
85  };
86 
87 
91  explicit SolverBFGS(SolverControl & residual_control,
92  const AdditionalData &data = AdditionalData());
93 
107  void
108  solve(
109  const std::function<Number(const VectorType &x, VectorType &g)> &compute,
110  VectorType & x);
111 
120  boost::signals2::connection
122  const std::function<
123  Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
124  &slot);
125 
152  boost::signals2::connection
154  const std::function<void(VectorType & g,
156  const FiniteSizeHistory<VectorType> &y)> &slot);
157 
158 
159 protected:
164 
168  boost::signals2::signal<
169  Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
171 
175  boost::signals2::signal<void(VectorType & g,
179 };
180 
181 
182 // ------------------- inline and template functions ----------------
183 #ifndef DOXYGEN
184 
185 template <typename VectorType>
187  const unsigned int max_history_size_,
188  const bool debug_output_)
189  : max_history_size(max_history_size_)
190  , debug_output(debug_output_)
191 {}
192 
193 
194 
195 template <typename VectorType>
197  const AdditionalData &data)
198  : SolverBase<VectorType>(solver_control)
199  , additional_data(data)
200 {}
201 
202 
203 
204 template <class VectorType>
205 boost::signals2::connection
207  const std::function<
208  Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> &slot)
209 {
210  Assert(line_search_signal.empty(),
211  ExcMessage("One should not attach more than one line search signal."));
212  return line_search_signal.connect(slot);
213 }
214 
215 
216 
217 template <class VectorType>
218 boost::signals2::connection
220  const std::function<void(VectorType & g,
222  const FiniteSizeHistory<VectorType> &y)> &slot)
223 {
224  Assert(preconditioner_signal.empty(),
225  ExcMessage(
226  "One should not attach more than one preconditioner signal."));
227  return preconditioner_signal.connect(slot);
228 }
229 
230 
231 
232 template <typename VectorType>
233 void
235  const std::function<typename VectorType::value_type(const VectorType &x,
236  VectorType &f)> &compute,
237  VectorType & x)
238 {
239  // Also see scipy Fortran implementation
240  // https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f
241  // and Octave-optim implementation:
242  // https://sourceforge.net/p/octave/optim/ci/default/tree/src/__bfgsmin.cc
243  LogStream::Prefix prefix("BFGS");
244 
245  // default line search:
246  bool first_step = true;
247  Number f_prev = 0.;
248  // provide default line search if no signal was attached
249  VectorType x0;
250  if (line_search_signal.empty())
251  {
252  x0.reinit(x);
253  const auto default_line_min =
254  [&](Number &f, VectorType &x, VectorType &g, const VectorType &p) {
255  const Number f0 = f;
256  const Number g0 = g * p;
257  Assert(g0 < 0,
258  ExcMessage(
259  "Function does not decrease along the current direction"));
260 
261  // save current solution value (to be used in line_search):
262  x0 = x;
263 
264  // see scipy implementation
265  // https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.line_search.html#scipy.optimize.line_search
266  // and Eq. 2.6.8 in Fletcher 2013, Practical methods of optimization
267  Number df = f_prev - f;
268  Assert(first_step || df >= 0.,
269  ExcMessage("Function value is not decreasing"));
270  df = std::max(df, 100. * std::numeric_limits<Number>::epsilon());
271  // guess a reasonable first step:
272  const Number a1 =
273  (first_step ? 1. : std::min(1., -1.01 * 2. * df / g0));
274  Assert(a1 > 0., ExcInternalError());
275  f_prev = f;
276 
277  // 1D line-search function
278  const auto line_func =
279  [&](const Number &x_line) -> std::pair<Number, Number> {
280  x = x0;
281  x.add(x_line, p);
282  f = compute(x, g);
283  const Number g_line = g * p;
284  return std::make_pair(f, g_line);
285  };
286 
287  // loose line search:
288  const auto res = LineMinimization::line_search<Number>(
289  line_func,
290  f0,
291  g0,
292  LineMinimization::poly_fit<Number>,
293  a1,
294  0.9,
295  0.001);
296 
297  if (first_step)
298  first_step = false;
299 
300  return res.first;
301  };
302  this->connect_line_search_slot(default_line_min);
303  }
304 
305  // FIXME: Octave has convergence in terms of:
306  // function change tolerance, default 1e-12
307  // parameter change tolerance, default 1e-6
308  // gradient tolerance, default 1e-5
309  // SolverBase and/or SolverControl need extension
310 
311  VectorType g(x), p(x), y_k(x), s_k(x);
312 
313  std::vector<Number> c1;
314  c1.reserve(additional_data.max_history_size);
315 
316  // limited history
317  FiniteSizeHistory<VectorType> y(additional_data.max_history_size);
318  FiniteSizeHistory<VectorType> s(additional_data.max_history_size);
319  FiniteSizeHistory<Number> rho(additional_data.max_history_size);
320 
321  unsigned int m = 0;
322  Number f;
323 
325  unsigned int k = 0;
326 
327  f = compute(x, g);
328 
329  conv = this->iteration_status(k, g.l2_norm(), x);
330  if (conv != SolverControl::iterate)
331  return;
332 
333  while (conv == SolverControl::iterate)
334  {
335  if (additional_data.debug_output)
336  deallog << "Iteration " << k << " history " << m << std::endl
337  << "f=" << f << std::endl;
338 
339  // 1. Two loop recursion to calculate p = - H*g
340  c1.resize(m);
341  p = g;
342  // first loop:
343  for (unsigned int i = 0; i < m; ++i)
344  {
345  c1[i] = rho[i] * (s[i] * p);
346  p.add(-c1[i], y[i]);
347  }
348  // H0
349  if (!preconditioner_signal.empty())
350  preconditioner_signal(p, s, y);
351 
352  // second loop:
353  for (int i = m - 1; i >= 0; --i)
354  {
355  Assert(i >= 0, ExcInternalError());
356  const Number c2 = rho[i] * (y[i] * p);
357  p.add(c1[i] - c2, s[i]);
358  }
359  p *= -1.;
360 
361  // 2. Line search
362  s_k = x;
363  y_k = g;
364  const Number alpha = line_search_signal(f, x, g, p)
365  .get(); // <-- signals return boost::optional
366  s_k.sadd(-1, 1, x);
367  y_k.sadd(-1, 1, g);
368 
369  if (additional_data.debug_output)
370  deallog << "Line search a=" << alpha << " f=" << f << std::endl;
371 
372  // 3. Check convergence
373  k++;
374  const Number g_l2 = g.l2_norm();
375  conv = this->iteration_status(k, g_l2, x);
376  if (conv != SolverControl::iterate)
377  break;
378 
379  // 4. Store s, y, rho
380  const Number curvature = s_k * y_k;
381  if (additional_data.debug_output)
382  deallog << "Curvature " << curvature << std::endl;
383 
384  if (curvature > 0. && additional_data.max_history_size > 0)
385  {
386  s.add(s_k);
387  y.add(y_k);
388  rho.add(1. / curvature);
389  m = s.size();
390 
391  Assert(y.size() == m, ExcInternalError());
392  Assert(rho.size() == m, ExcInternalError());
393  }
394 
395  Assert(m <= additional_data.max_history_size, ExcInternalError());
396  }
397 
398  // In the case of failure: throw exception.
400  SolverControl::NoConvergence(k, g.l2_norm()));
401 }
402 
403 #endif
404 
405 DEAL_II_NAMESPACE_CLOSE
406 
407 #endif
boost::signals2::signal< Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> line_search_signal
Definition: solver_bfgs.h:170
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
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
unsigned int max_history_size
Definition: solver_bfgs.h:79
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:163
#define Assert(cond, exc)
Definition: exceptions.h:1407
VectorType::value_type Number
Definition: solver_bfgs.h:62
boost::signals2::signal< void(VectorType &g, const FiniteSizeHistory< VectorType > &s, const FiniteSizeHistory< VectorType > &y)> preconditioner_signal
Definition: solver_bfgs.h:178
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()