deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14: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 <boost/signals2.hpp>
27
28#include <limits>
29
31
58template <typename VectorType>
59class SolverBFGS : public SolverBase<VectorType>
60{
61public:
65 using Number = typename VectorType::value_type;
66
67
72 {
76 explicit AdditionalData(const unsigned int max_history_size = 5,
77 const bool debug_output = false);
78
82 unsigned int max_history_size;
83
88 };
89
90
96
110 void
112 const std::function<Number(const VectorType &x, VectorType &g)> &compute,
113 VectorType &x);
114
123 boost::signals2::connection
125 const std::function<
126 Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
127 &slot);
128
155 boost::signals2::connection
157 const std::function<void(VectorType &g,
160
161
162protected:
167
171 boost::signals2::signal<
172 Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)>
174
178 boost::signals2::signal<void(VectorType &g,
182};
183
184
185// ------------------- inline and template functions ----------------
186#ifndef DOXYGEN
187
188template <typename VectorType>
190 const unsigned int max_history_size_,
191 const bool debug_output_)
192 : max_history_size(max_history_size_)
193 , debug_output(debug_output_)
194{}
195
196
197
198template <typename VectorType>
200 const AdditionalData &data)
201 : SolverBase<VectorType>(solver_control)
202 , additional_data(data)
203{}
204
205
206
207template <typename VectorType>
208boost::signals2::connection
210 const std::function<
211 Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> &slot)
212{
213 Assert(line_search_signal.empty(),
214 ExcMessage("One should not attach more than one line search signal."));
215 return line_search_signal.connect(slot);
216}
217
218
219
220template <typename VectorType>
221boost::signals2::connection
223 const std::function<void(VectorType &g,
226{
227 Assert(preconditioner_signal.empty(),
229 "One should not attach more than one preconditioner signal."));
230 return preconditioner_signal.connect(slot);
231}
232
233
234
235template <typename VectorType>
236void
238 const std::function<typename VectorType::value_type(const VectorType &x,
239 VectorType &f)> &compute,
240 VectorType &x)
241{
242 // Also see scipy Fortran implementation
243 // https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f
244 // and Octave-optim implementation:
245 // https://sourceforge.net/p/octave/optim/ci/default/tree/src/__bfgsmin.cc
247
248 // default line search:
249 bool first_step = true;
250 Number f_prev = 0.;
251 // provide default line search if no signal was attached
253 if (line_search_signal.empty())
254 {
255 x0.reinit(x);
256 const auto default_line_min =
257 [&](Number &f, VectorType &x, VectorType &g, const VectorType &p) {
258 const Number f0 = f;
259 const Number g0 = g * p;
260 Assert(g0 < 0,
262 "Function does not decrease along the current direction"));
263
264 // save current solution value (to be used in line_search):
265 x0 = x;
266
267 // see scipy implementation
268 // https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.line_search.html#scipy.optimize.line_search
269 // and Eq. 2.6.8 in Fletcher 2013, Practical methods of optimization
270 Number df = f_prev - f;
271 Assert(first_step || df >= 0.,
272 ExcMessage("Function value is not decreasing"));
273 df = std::max(df, 100. * std::numeric_limits<Number>::epsilon());
274 // guess a reasonable first step:
275 const Number a1 =
276 (first_step ? 1. : std::min(1., -1.01 * 2. * df / g0));
277 Assert(a1 > 0., ExcInternalError());
278 f_prev = f;
279
280 // 1d line-search function
281 const auto line_func =
282 [&](const Number &x_line) -> std::pair<Number, Number> {
283 x = x0;
284 x.add(x_line, p);
285 f = compute(x, g);
286 const Number g_line = g * p;
287 return std::make_pair(f, g_line);
288 };
289
290 // loose line search:
291 const auto res = LineMinimization::line_search<Number>(
292 line_func,
293 f0,
294 g0,
295 LineMinimization::poly_fit<Number>,
296 a1,
297 0.9,
298 0.001);
299
300 if (first_step)
301 first_step = false;
302
303 return res.first;
304 };
305 this->connect_line_search_slot(default_line_min);
306 }
307
308 // FIXME: Octave has convergence in terms of:
309 // function change tolerance, default 1e-12
310 // parameter change tolerance, default 1e-6
311 // gradient tolerance, default 1e-5
312 // SolverBase and/or SolverControl need extension
313
314 VectorType g(x), p(x), y_k(x), s_k(x);
315
316 std::vector<Number> c1;
317 c1.reserve(additional_data.max_history_size);
318
319 // limited history
320 FiniteSizeHistory<VectorType> y(additional_data.max_history_size);
321 FiniteSizeHistory<VectorType> s(additional_data.max_history_size);
322 FiniteSizeHistory<Number> rho(additional_data.max_history_size);
323
324 unsigned int m = 0;
325 Number f;
326
328 unsigned int k = 0;
329
330 f = compute(x, g);
331
332 conv = this->iteration_status(k, g.l2_norm(), x);
334 return;
335
337 {
338 if (additional_data.debug_output)
339 deallog << "Iteration " << k << " history " << m << std::endl
340 << "f=" << f << std::endl;
341
342 // 1. Two loop recursion to calculate p = - H*g
343 c1.resize(m);
344 p = g;
345 // first loop:
346 for (unsigned int i = 0; i < m; ++i)
347 {
348 c1[i] = rho[i] * (s[i] * p);
349 p.add(-c1[i], y[i]);
350 }
351 // H0
352 if (!preconditioner_signal.empty())
353 preconditioner_signal(p, s, y);
354
355 // second loop:
356 for (int i = m - 1; i >= 0; --i)
357 {
358 Assert(i >= 0, ExcInternalError());
359 const Number c2 = rho[i] * (y[i] * p);
360 p.add(c1[i] - c2, s[i]);
361 }
362 p *= -1.;
363
364 // 2. Line search
365 s_k = x;
366 y_k = g;
367 const Number alpha = line_search_signal(f, x, g, p)
368 .get(); // <-- signals return boost::optional
369 s_k.sadd(-1, 1, x);
370 y_k.sadd(-1, 1, g);
371
372 if (additional_data.debug_output)
373 deallog << "Line search a=" << alpha << " f=" << f << std::endl;
374
375 // 3. Check convergence
376 ++k;
377 const Number g_l2 = g.l2_norm();
378 conv = this->iteration_status(k, g_l2, x);
380 break;
381
382 // 4. Store s, y, rho
383 const Number curvature = s_k * y_k;
384 if (additional_data.debug_output)
385 deallog << "Curvature " << curvature << std::endl;
386
387 if (curvature > 0. && additional_data.max_history_size > 0)
388 {
389 s.add(s_k);
390 y.add(y_k);
391 rho.add(1. / curvature);
392 m = s.size();
393
394 Assert(y.size() == m, ExcInternalError());
395 Assert(rho.size() == m, ExcInternalError());
396 }
397
398 Assert(m <= additional_data.max_history_size, ExcInternalError());
399 }
400
401 // In the case of failure: throw exception.
403 SolverControl::NoConvergence(k, g.l2_norm()));
404}
405
406#endif
407
409
410#endif
typename VectorType::value_type Number
Definition solver_bfgs.h:65
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.
friend class Tensor
Definition tensor.h:865
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:38
std::vector< index_type > data
Definition mpi.cc:740
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
::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:82