Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
solver_bfgs.h
Go to the documentation of this file.
1//-----------------------------------------------------------
2//
3// Copyright (C) 2018 - 2023 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
24
26
27#include <limits>
28
30
57template <typename VectorType>
58class SolverBFGS : public SolverBase<VectorType>
59{
60public:
64 using Number = typename VectorType::value_type;
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
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
161protected:
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
187template <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
197template <typename VectorType>
199 const AdditionalData &data)
200 : SolverBase<VectorType>(solver_control)
201 , additional_data(data)
202{}
203
204
205
206template <class VectorType>
207boost::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
219template <class VectorType>
220boost::signals2::connection
222 const std::function<void(VectorType & g,
224 const FiniteSizeHistory<VectorType> &y)> &slot)
225{
226 Assert(preconditioner_signal.empty(),
228 "One should not attach more than one preconditioner signal."));
229 return preconditioner_signal.connect(slot);
230}
231
232
233
234template <typename VectorType>
235void
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,
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
319 FiniteSizeHistory<VectorType> y(additional_data.max_history_size);
320 FiniteSizeHistory<VectorType> s(additional_data.max_history_size);
321 FiniteSizeHistory<Number> rho(additional_data.max_history_size);
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 {
337 if (additional_data.debug_output)
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
371 if (additional_data.debug_output)
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;
383 if (additional_data.debug_output)
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
397 Assert(m <= additional_data.max_history_size, ExcInternalError());
398 }
399
400 // In the case of failure: throw exception.
402 SolverControl::NoConvergence(k, g.l2_norm()));
403}
404
405#endif
406
408
409#endif
std::size_t size() const
void add(const T &element)
typename VectorType::value_type Number
Definition solver_bfgs.h:64
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)
void solve(const std::function< Number(const VectorType &x, VectorType &g)> &compute, VectorType &x)
boost::signals2::signal< void(VectorType &g, const FiniteSizeHistory< VectorType > &s, const FiniteSizeHistory< VectorType > &y)> preconditioner_signal
const AdditionalData additional_data
boost::signals2::signal< Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> line_search_signal
SolverBFGS(SolverControl &residual_control, const AdditionalData &data=AdditionalData())
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:37
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_history_size=5, const bool debug_output=false)
unsigned int max_history_size
Definition solver_bfgs.h:81