Reference documentation for deal.II version GIT relicensing-756-gc37f912904 2024-05-25 08:20:01+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_bfgs_h
16#define dealii_solver_bfgs_h
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/lac/solver.h>
21
23
25
26#include <limits>
27
29
56template <typename VectorType>
57class SolverBFGS : public SolverBase<VectorType>
58{
59public:
63 using Number = typename VectorType::value_type;
64
65
70 {
74 explicit AdditionalData(const unsigned int max_history_size = 5,
75 const bool debug_output = false);
76
80 unsigned int max_history_size;
81
86 };
87
88
92 explicit SolverBFGS(SolverControl &residual_control,
93 const AdditionalData &data = AdditionalData());
94
108 void
110 const std::function<Number(const VectorType &x, VectorType &g)> &compute,
111 VectorType &x);
112
121 boost::signals2::connection
123 const std::function<
124 Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
125 &slot);
126
153 boost::signals2::connection
155 const std::function<void(VectorType &g,
157 const FiniteSizeHistory<VectorType> &y)> &slot);
158
159
160protected:
165
169 boost::signals2::signal<
170 Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
172
176 boost::signals2::signal<void(VectorType &g,
180};
181
182
183// ------------------- inline and template functions ----------------
184#ifndef DOXYGEN
185
186template <typename VectorType>
188 const unsigned int max_history_size_,
189 const bool debug_output_)
190 : max_history_size(max_history_size_)
191 , debug_output(debug_output_)
192{}
193
194
195
196template <typename VectorType>
198 const AdditionalData &data)
199 : SolverBase<VectorType>(solver_control)
200 , additional_data(data)
201{}
202
203
204
205template <typename VectorType>
206boost::signals2::connection
208 const std::function<
209 Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> &slot)
210{
211 Assert(line_search_signal.empty(),
212 ExcMessage("One should not attach more than one line search signal."));
213 return line_search_signal.connect(slot);
214}
215
216
217
218template <typename VectorType>
219boost::signals2::connection
221 const std::function<void(VectorType &g,
223 const FiniteSizeHistory<VectorType> &y)> &slot)
224{
225 Assert(preconditioner_signal.empty(),
227 "One should not attach more than one preconditioner signal."));
228 return preconditioner_signal.connect(slot);
229}
230
231
232
233template <typename VectorType>
234void
236 const std::function<typename VectorType::value_type(const VectorType &x,
237 VectorType &f)> &compute,
238 VectorType &x)
239{
240 // Also see scipy Fortran implementation
241 // https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f
242 // and Octave-optim implementation:
243 // https://sourceforge.net/p/octave/optim/ci/default/tree/src/__bfgsmin.cc
244 LogStream::Prefix prefix("BFGS");
245
246 // default line search:
247 bool first_step = true;
248 Number f_prev = 0.;
249 // provide default line search if no signal was attached
250 VectorType x0;
251 if (line_search_signal.empty())
252 {
253 x0.reinit(x);
254 const auto default_line_min =
255 [&](Number &f, VectorType &x, VectorType &g, const VectorType &p) {
256 const Number f0 = f;
257 const Number g0 = g * p;
258 Assert(g0 < 0,
260 "Function does not decrease along the current direction"));
261
262 // save current solution value (to be used in line_search):
263 x0 = x;
264
265 // see scipy implementation
266 // https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.line_search.html#scipy.optimize.line_search
267 // and Eq. 2.6.8 in Fletcher 2013, Practical methods of optimization
268 Number df = f_prev - f;
269 Assert(first_step || df >= 0.,
270 ExcMessage("Function value is not decreasing"));
271 df = std::max(df, 100. * std::numeric_limits<Number>::epsilon());
272 // guess a reasonable first step:
273 const Number a1 =
274 (first_step ? 1. : std::min(1., -1.01 * 2. * df / g0));
275 Assert(a1 > 0., ExcInternalError());
276 f_prev = f;
277
278 // 1d line-search function
279 const auto line_func =
280 [&](const Number &x_line) -> std::pair<Number, Number> {
281 x = x0;
282 x.add(x_line, p);
283 f = compute(x, g);
284 const Number g_line = g * p;
285 return std::make_pair(f, g_line);
286 };
287
288 // loose line search:
289 const auto res = LineMinimization::line_search<Number>(
290 line_func,
291 f0,
292 g0,
293 LineMinimization::poly_fit<Number>,
294 a1,
295 0.9,
296 0.001);
297
298 if (first_step)
299 first_step = false;
300
301 return res.first;
302 };
303 this->connect_line_search_slot(default_line_min);
304 }
305
306 // FIXME: Octave has convergence in terms of:
307 // function change tolerance, default 1e-12
308 // parameter change tolerance, default 1e-6
309 // gradient tolerance, default 1e-5
310 // SolverBase and/or SolverControl need extension
311
312 VectorType g(x), p(x), y_k(x), s_k(x);
313
314 std::vector<Number> c1;
315 c1.reserve(additional_data.max_history_size);
316
317 // limited history
318 FiniteSizeHistory<VectorType> y(additional_data.max_history_size);
319 FiniteSizeHistory<VectorType> s(additional_data.max_history_size);
320 FiniteSizeHistory<Number> rho(additional_data.max_history_size);
321
322 unsigned int m = 0;
323 Number f;
324
326 unsigned int k = 0;
327
328 f = compute(x, g);
329
330 conv = this->iteration_status(k, g.l2_norm(), x);
331 if (conv != SolverControl::iterate)
332 return;
333
334 while (conv == SolverControl::iterate)
335 {
336 if (additional_data.debug_output)
337 deallog << "Iteration " << k << " history " << m << std::endl
338 << "f=" << f << std::endl;
339
340 // 1. Two loop recursion to calculate p = - H*g
341 c1.resize(m);
342 p = g;
343 // first loop:
344 for (unsigned int i = 0; i < m; ++i)
345 {
346 c1[i] = rho[i] * (s[i] * p);
347 p.add(-c1[i], y[i]);
348 }
349 // H0
350 if (!preconditioner_signal.empty())
351 preconditioner_signal(p, s, y);
352
353 // second loop:
354 for (int i = m - 1; i >= 0; --i)
355 {
356 Assert(i >= 0, ExcInternalError());
357 const Number c2 = rho[i] * (y[i] * p);
358 p.add(c1[i] - c2, s[i]);
359 }
360 p *= -1.;
361
362 // 2. Line search
363 s_k = x;
364 y_k = g;
365 const Number alpha = line_search_signal(f, x, g, p)
366 .get(); // <-- signals return boost::optional
367 s_k.sadd(-1, 1, x);
368 y_k.sadd(-1, 1, g);
369
370 if (additional_data.debug_output)
371 deallog << "Line search a=" << alpha << " f=" << f << std::endl;
372
373 // 3. Check convergence
374 ++k;
375 const Number g_l2 = g.l2_norm();
376 conv = this->iteration_status(k, g_l2, x);
377 if (conv != SolverControl::iterate)
378 break;
379
380 // 4. Store s, y, rho
381 const Number curvature = s_k * y_k;
382 if (additional_data.debug_output)
383 deallog << "Curvature " << curvature << std::endl;
384
385 if (curvature > 0. && additional_data.max_history_size > 0)
386 {
387 s.add(s_k);
388 y.add(y_k);
389 rho.add(1. / curvature);
390 m = s.size();
391
392 Assert(y.size() == m, ExcInternalError());
393 Assert(rho.size() == m, ExcInternalError());
394 }
395
396 Assert(m <= additional_data.max_history_size, ExcInternalError());
397 }
398
399 // In the case of failure: throw exception.
401 SolverControl::NoConvergence(k, g.l2_norm()));
402}
403
404#endif
405
407
408#endif
std::size_t size() const
void add(const T &element)
typename VectorType::value_type Number
Definition solver_bfgs.h:63
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:36
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
::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:80