Skip to content

Commit 0cfddd2

Browse files
AlePaluAlePalu
committed
optimization module refactored, removed callbacks as class template arguments, all optimizers (both gradient and gradient-free) expose the same API
1 parent 055e8da commit 0cfddd2

File tree

12 files changed

+499
-378
lines changed

12 files changed

+499
-378
lines changed

‎fdaPDE/optimization.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -26,11 +26,11 @@
2626

2727
// callbacks
2828
#include "src/optimization/callbacks.h"
29-
#include "src/optimization/backtracking_line_search.h"
30-
#include "src/optimization/wolfe_line_search.h"
29+
#include "src/optimization/backtracking.h"
30+
#include "src/optimization/wolfe.h"
3131

3232
// algorithms
33-
#include "src/optimization/grid.h"
33+
#include "src/optimization/grid_search.h"
3434
#include "src/optimization/newton.h"
3535
#include "src/optimization/gradient_descent.h"
3636
#include "src/optimization/conjugate_gradient.h"

‎fdaPDE/src/fields/matrix_field.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -830,7 +830,6 @@ template <int StaticInputSize, typename Derived> struct MatrixFieldBase {
830830
if constexpr (
831831
Derived::Rows == Dynamic || Derived::Cols == Dynamic || InputType_::RowsAtCompileTime == Dynamic ||
832832
InputType_::ColsAtCompileTime == Dynamic) {
833-
fdapde_assert(p.rows() == derived().rows() && p.cols() == derived().cols());
834833
out.resize(derived().rows(), derived().cols());
835834
}
836835
eval_at(p, out);

‎fdaPDE/src/optimization/backtracking_line_search.h‎ renamed to ‎fdaPDE/src/optimization/backtracking.h‎

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -24,21 +24,20 @@ namespace fdapde {
2424
// implementation of the backatracking line search method for step selection
2525
class BacktrackingLineSearch {
2626
private:
27-
double alpha_ = 2.0;
28-
double beta_ = 0.5;
29-
double gamma_ = 0.5;
27+
double alpha_, beta_, gamma_;
3028
public:
3129
// constructors
32-
BacktrackingLineSearch() = default;
30+
BacktrackingLineSearch() : alpha_(2.0), beta_(0.5), gamma_(0.5) { }
3331
BacktrackingLineSearch(double alpha, double beta, double gamma) : alpha_(alpha), beta_(beta), gamma_(gamma) { }
3432

3533
// backtracking based step search
36-
template <typename Opt, typename Obj> bool pre_update_step(Opt& opt, Obj& obj) {
34+
template <typename Opt, typename Obj> bool adapt_hook(Opt& opt, Obj& obj) {
3735
double alpha = alpha_; // restore to user defined settings
3836
double m = opt.grad_old.dot(opt.update);
3937
if (m < 0) { // descent direction
4038
while (obj(opt.x_old) - obj(opt.x_old + alpha * opt.update) // Armijo–Goldstein condition
41-
+ gamma_ * alpha * m < 0) {
39+
+ gamma_ * alpha * m <
40+
0) {
4241
alpha *= beta_;
4342
}
4443
}

‎fdaPDE/src/optimization/bfgs.h‎

Lines changed: 36 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -21,110 +21,92 @@
2121

2222
namespace fdapde {
2323

24-
// implementation of Broyden–Fletcher–Goldfarb–Shanno algorithm for unconstrained nonlinear optimization
25-
template <int N, typename... Args> class BFGS {
24+
template <int N> class BFGS {
2625
private:
2726
using vector_t = std::conditional_t<N == Dynamic, Eigen::Matrix<double, Dynamic, 1>, Eigen::Matrix<double, N, 1>>;
2827
using matrix_t =
2928
std::conditional_t<N == Dynamic, Eigen::Matrix<double, Dynamic, Dynamic>, Eigen::Matrix<double, N, N>>;
3029

31-
std::tuple<Args...> callbacks_;
3230
vector_t optimum_;
33-
double value_; // objective value at optimum
31+
double value_; // objective value at optimum
32+
int n_iter_ = 0; // current iteration number
33+
std::vector<double> values_; // explored objective values during optimization
34+
matrix_t inv_hessian_;
35+
3436
int max_iter_; // maximum number of iterations before forced stop
35-
int n_iter_ = 0; // current iteration number
3637
double tol_; // tolerance on error before forced stop
3738
double step_; // update step
3839
public:
3940
vector_t x_old, x_new, update, grad_old, grad_new;
40-
matrix_t inv_hessian;
4141
double h;
4242
// constructor
43-
BFGS() = default;
44-
BFGS(int max_iter, double tol, double step)
45-
requires(sizeof...(Args) != 0)
46-
: max_iter_(max_iter), tol_(tol), step_(step) { }
47-
BFGS(int max_iter, double tol, double step, Args&&... callbacks) :
48-
callbacks_(std::make_tuple(std::forward<Args>(callbacks)...)), max_iter_(max_iter), tol_(tol), step_(step) { }
49-
// copy semantic
50-
BFGS(const BFGS& other) :
51-
callbacks_(other.callbacks_), max_iter_(other.max_iter_), tol_(other.tol_), step_(other.step_) { }
43+
BFGS() : max_iter_(500), tol_(1e-5), step_(1e-2) { }
44+
BFGS(int max_iter, double tol, double step) : max_iter_(max_iter), tol_(tol), step_(step) { }
45+
BFGS(const BFGS& other) : max_iter_(other.max_iter_), tol_(other.tol_), step_(other.step_) { }
5246
BFGS& operator=(const BFGS& other) {
5347
max_iter_ = other.max_iter_;
5448
tol_ = other.tol_;
5549
step_ = other.step_;
56-
callbacks_ = other.callbacks_;
5750
return *this;
5851
}
59-
template <typename ObjectiveT, typename... Functor>
60-
requires(sizeof...(Functor) < 2) && ((requires(Functor f, double value) { f(value); }) && ...)
61-
vector_t optimize(ObjectiveT&& objective, const vector_t& x0, Functor&&... func) {
52+
template <typename ObjectiveT, typename... Callbacks>
53+
vector_t optimize(ObjectiveT&& objective, const vector_t& x0, Callbacks&&... callbacks) {
6254
fdapde_static_assert(
6355
std::is_same<decltype(std::declval<ObjectiveT>().operator()(vector_t())) FDAPDE_COMMA double>::value,
64-
INVALID_CALL_TO_OPTIMIZE__OBJECTIVE_FUNCTOR_NOT_ACCEPTING_VECTORTYPE);
56+
INVALID_CALL_TO_OPTIMIZE__OBJECTIVE_FUNCTOR_NOT_CALLABLE_AT_VECTOR_TYPE);
57+
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
58+
std::tuple<Callbacks...> callbacks_ {callbacks...};
6559
bool stop = false; // asserted true in case of forced stop
66-
vector_t zero;
67-
double error = 0;
68-
auto grad = objective.gradient();
69-
n_iter_ = 0;
60+
double error = std::numeric_limits<double>::max();
61+
int size = N == Dynamic ? x0.rows() : N;
62+
auto grad = objective.gradient();
7063
h = step_;
71-
x_old = x0, x_new = x0;
64+
n_iter_ = 0;
65+
x_old = x0, x_new = vector_t::Constant(size, NaN);
66+
grad_old = grad(x_old), grad_new = vector_t::Constant(size, NaN);
7267
if constexpr (N == Dynamic) { // inv_hessian approximated with identity matrix
73-
inv_hessian = matrix_t::Identity(x0.rows(), x0.rows());
74-
zero = vector_t::Zero(x0.rows());
68+
inv_hessian_ = matrix_t::Identity(size, size);
7569
} else {
76-
inv_hessian = matrix_t::Identity();
77-
zero = vector_t::Zero();
70+
inv_hessian_ = matrix_t::Identity();
7871
}
79-
grad_old = grad(x_old);
80-
if (grad_old.isApprox(zero)) { // already at stationary point
81-
optimum_ = x_old;
82-
value_ = objective(optimum_);
83-
if constexpr (sizeof...(Functor) == 1) { (func(value_), ...); }
84-
return optimum_;
85-
}
72+
update = -inv_hessian_ * grad_old;
73+
stop |= internals::exec_grad_hooks(*this, objective, callbacks_);
8674
error = grad_old.norm();
75+
values_.push_back(objective(x_old));
8776

8877
while (n_iter_ < max_iter_ && error > tol_ && !stop) {
89-
// compute update direction
90-
update = -inv_hessian * grad_old;
91-
stop |= execute_pre_update_step(*this, objective, callbacks_);
78+
stop |= internals::exec_adapt_hooks(*this, objective, callbacks_);
9279
// update along descent direction
9380
x_new = x_old + h * update;
9481
grad_new = grad(x_new);
95-
if (grad_new.isApprox(zero)) { // already at stationary point
96-
optimum_ = x_old;
97-
value_ = objective(optimum_);
98-
if constexpr (sizeof...(Functor) == 1) { (func(value_), ...); }
99-
return optimum_;
100-
}
10182
// update inverse hessian approximation
10283
vector_t delta_x = x_new - x_old;
10384
vector_t delta_grad = grad_new - grad_old;
10485
double xg = delta_x.dot(delta_grad);
105-
vector_t hx = inv_hessian * delta_grad;
86+
vector_t hx = inv_hessian_ * delta_grad;
10687

10788
matrix_t U = (1 + (delta_grad.dot(hx)) / xg) * ((delta_x * delta_x.transpose()) / xg);
10889
matrix_t V = ((hx * delta_x.transpose() + delta_x * hx.transpose())) / xg;
109-
inv_hessian += (U - V);
90+
inv_hessian_ += (U - V);
11091
// prepare next iteration
111-
if constexpr (sizeof...(Functor) == 1) { (func(objective(x_old)), ...); }
112-
error = grad_new.norm();
92+
update = -inv_hessian_ * grad_new;
11393
stop |=
114-
(execute_post_update_step(*this, objective, callbacks_) || execute_stopping_criterion(*this, objective));
94+
(internals::exec_grad_hooks(*this, objective, callbacks_) || internals::exec_stop_if(*this, objective));
11595
x_old = x_new;
11696
grad_old = grad_new;
97+
error = grad_new.norm();
98+
values_.push_back(objective(x_old));
11799
n_iter_++;
118100
}
119101
optimum_ = x_old;
120-
value_ = objective(optimum_);
121-
if constexpr (sizeof...(Functor) == 1) { (func(value_), ...); }
102+
value_ = values_.back();
122103
return optimum_;
123104
}
124-
// getters
125-
vector_t optimum() const { return optimum_; }
105+
// observers
106+
const vector_t& optimum() const { return optimum_; }
126107
double value() const { return value_; }
127108
int n_iter() const { return n_iter_; }
109+
const std::vector<double>& values() const { return values_; }
128110
};
129111

130112
} // namespace fdapde

‎fdaPDE/src/optimization/callbacks.h‎

Lines changed: 45 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -20,36 +20,58 @@
2020
#include "header_check.h"
2121

2222
namespace fdapde {
23+
namespace internals {
2324

24-
template <typename Opt, typename Obj, typename... Args>
25-
bool execute_pre_update_step(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
25+
template <typename Hook, typename... Args>
26+
bool opt_hooks_loop(Hook hook, std::tuple<Args...>& callbacks) {
2627
bool b = false;
27-
auto exec_callback = [&](auto&& callback) {
28-
if constexpr (requires(std::decay_t<decltype(callback)> c, Opt opt, Obj obj) {
29-
{ c.pre_update_step(opt, obj) } -> std::same_as<bool>;
30-
}) {
31-
b |= callback.pre_update_step(optimizer, objective);
32-
}
33-
};
34-
std::apply([&](auto&&... callback) { (exec_callback(callback), ...); }, callbacks);
28+
std::apply([&](auto&&... callback) { ([&]() { b |= hook(callback); }(), ...); }, callbacks);
3529
return b;
3630
}
3731

3832
template <typename Opt, typename Obj, typename... Args>
39-
bool execute_post_update_step(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
40-
bool b = false;
41-
auto exec_callback = [&](auto&& callback) {
42-
if constexpr (requires(std::decay_t<decltype(callback)> c, Opt opt, Obj obj) {
43-
{ c.post_update_step(opt, obj) } -> std::same_as<bool>;
44-
}) {
45-
b |= callback.post_update_step(optimizer, objective);
46-
}
47-
};
48-
std::apply([&](auto&&... callback) { (exec_callback(callback), ...); }, callbacks);
49-
return b;
33+
bool exec_eval_hooks(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
34+
return opt_hooks_loop(
35+
[&](auto&& callback) {
36+
if constexpr (requires(std::decay_t<decltype(callback)> c, Opt opt, Obj obj) {
37+
{ c.eval_hook(opt, obj) } -> std::same_as<bool>;
38+
}) {
39+
return callback.eval_hook(optimizer, objective);
40+
}
41+
return false;
42+
},
43+
callbacks);
44+
}
45+
46+
template <typename Opt, typename Obj, typename... Args>
47+
bool exec_grad_hooks(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
48+
return opt_hooks_loop(
49+
[&](auto&& callback) {
50+
if constexpr (requires(std::decay_t<decltype(callback)> c, Opt opt, Obj obj) {
51+
{ c.grad_hook(opt, obj) } -> std::same_as<bool>;
52+
}) {
53+
return callback.grad_hook(optimizer, objective);
54+
}
55+
return false;
56+
},
57+
callbacks);
58+
}
59+
60+
template <typename Opt, typename Obj, typename... Args>
61+
bool exec_adapt_hooks(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
62+
return opt_hooks_loop(
63+
[&](auto&& callback) {
64+
if constexpr (requires(std::decay_t<decltype(callback)> c, Opt opt, Obj obj) {
65+
{ c.adapt_hook(opt, obj) } -> std::same_as<bool>;
66+
}) {
67+
return callback.adapt_hook(optimizer, objective);
68+
}
69+
return false;
70+
},
71+
callbacks);
5072
}
5173

52-
template <typename Opt, typename Obj> bool execute_stopping_criterion(Opt& optimizer, Obj& objective) {
74+
template <typename Opt, typename Obj> bool exec_stop_if(Opt& optimizer, Obj& objective) {
5375
bool b = false;
5476
if constexpr (requires(Opt opt, Obj obj) {
5577
{ obj.stop_if(opt) } -> std::same_as<bool>;
@@ -59,6 +81,7 @@ template <typename Opt, typename Obj> bool execute_stopping_criterion(Opt& optim
5981
return b;
6082
}
6183

84+
} // namespace internals
6285
} // namespace fdapde
6386

6487
#endif // __FDAPDE_OPTIMIZATION_CALLBACKS_H__

0 commit comments

Comments
 (0)