Skip to content

Commit 00e5312

Browse files
AlePaluAlePalu
committed
is_gradient_based_opt_v and is_gradient_free_opt_v concepts, nelder_mead callbacks support, gradient-free optimizers expose x_curr and obj_curr
1 parent 0cfddd2 commit 00e5312

File tree

10 files changed

+102
-213
lines changed

10 files changed

+102
-213
lines changed

‎fdaPDE/optimization.h‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,39 @@
2424
#include "utility.h"
2525
#include "fields.h"
2626

27+
namespace fdapde {
28+
29+
template <typename Opt> struct is_gradient_based_opt {
30+
private:
31+
using Opt_ = std::decay_t<Opt>;
32+
static constexpr int N = Opt_::static_input_size;
33+
using vector_t = std::conditional_t<N == Dynamic, Eigen::Matrix<double, Dynamic, 1>, Eigen::Matrix<double, N, 1>>;
34+
public:
35+
static constexpr bool value = !Opt_::gradient_free && requires(Opt_ opt) {
36+
{ opt.x_old } -> std::convertible_to<vector_t>;
37+
{ opt.x_new } -> std::convertible_to<vector_t>;
38+
{ opt.update } -> std::convertible_to<vector_t>;
39+
{ opt.grad_old } -> std::convertible_to<vector_t>;
40+
{ opt.grad_new } -> std::convertible_to<vector_t>;
41+
};
42+
};
43+
template <typename Opt> static constexpr bool is_gradient_based_opt_v = is_gradient_based_opt<Opt>::value;
44+
45+
template <typename Opt> struct is_gradient_free_opt {
46+
private:
47+
using Opt_ = std::decay_t<Opt>;
48+
static constexpr int N = Opt_::static_input_size;
49+
using vector_t = std::conditional_t<N == Dynamic, Eigen::Matrix<double, Dynamic, 1>, Eigen::Matrix<double, N, 1>>;
50+
public:
51+
static constexpr bool value = Opt_::gradient_free && requires(Opt_ opt) {
52+
{ opt.x_curr } -> std::convertible_to<vector_t>;
53+
{ opt.obj_curr } -> std::convertible_to<double>;
54+
};
55+
};
56+
template <typename Opt> static constexpr bool is_gradient_free_opt_v = is_gradient_free_opt<Opt>::value;
57+
58+
} // namespace fdapde
59+
2760
// callbacks
2861
#include "src/optimization/callbacks.h"
2962
#include "src/optimization/backtracking.h"

‎fdaPDE/src/optimization/backtracking.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121

2222
namespace fdapde {
2323

24-
// implementation of the backatracking line search method for step selection
24+
// backtracking line search method for adaptive step size
2525
class BacktrackingLineSearch {
2626
private:
2727
double alpha_, beta_, gamma_;
@@ -32,6 +32,7 @@ class BacktrackingLineSearch {
3232

3333
// backtracking based step search
3434
template <typename Opt, typename Obj> bool adapt_hook(Opt& opt, Obj& obj) {
35+
fdapde_static_assert(is_gradient_based_opt_v<Opt>, THIS_METHOD_IS_FOR_GRADIENT_BASED_OPTIMIZATION_ONLY);
3536
double alpha = alpha_; // restore to user defined settings
3637
double m = opt.grad_old.dot(opt.update);
3738
if (m < 0) { // descent direction

‎fdaPDE/src/optimization/bfgs.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,8 @@ template <int N> class BFGS {
3737
double tol_; // tolerance on error before forced stop
3838
double step_; // update step
3939
public:
40+
static constexpr bool gradient_free = false;
41+
static constexpr int static_input_size = N;
4042
vector_t x_old, x_new, update, grad_old, grad_new;
4143
double h;
4244
// constructor

‎fdaPDE/src/optimization/conjugate_gradient.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,8 @@ template <int N, typename DirectionUpdate> class conjugate_gradient_impl {
3737
double tol_; // tolerance on error before forced stop
3838
double step_; // update step
3939
public:
40+
static constexpr bool gradient_free = false;
41+
static constexpr int static_input_size = N;
4042
vector_t x_old, x_new, update, grad_old, grad_new;
4143
double h;
4244
// constructor

‎fdaPDE/src/optimization/gradient_descent.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,8 @@ template <int N> class GradientDescent {
3636
double tol_; // tolerance on error before forced stop
3737
double step_; // update step
3838
public:
39+
static constexpr bool gradient_free = false;
40+
static constexpr int static_input_size = N;
3941
vector_t x_old, x_new, update, grad_old, grad_new;
4042
double h;
4143
// constructor

‎fdaPDE/src/optimization/grid_search.h‎

Lines changed: 19 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -31,13 +31,13 @@ template <int N> class GridSearch {
3131
std::vector<double> values_; // explored objective values during optimization
3232
int size_;
3333
public:
34-
vector_t x_old, x_new;
35-
double obj_old, obj_new;
34+
static constexpr bool gradient_free = true;
35+
static constexpr int static_input_size = N;
36+
vector_t x_curr;
37+
double obj_curr;
3638
// constructor
37-
GridSearch()
38-
requires(N != Dynamic)
39-
: size_(N) { }
40-
GridSearch(int size) : size_(N == Dynamic ? size : N) { }
39+
GridSearch() : size_(N) { fdapde_static_assert(N != Dynamic, THIS_METHOD_IS_FOR_STATIC_SIZED_GRID_SEARCH_ONLY); }
40+
GridSearch(int size) : size_(N == Dynamic ? size : N) { fdapde_assert(N == Dynamic || size == N); }
4141
GridSearch(const GridSearch& other) : size_(other.size_) { }
4242
GridSearch& operator=(const GridSearch& other) {
4343
size_ = other.size_;
@@ -61,29 +61,26 @@ template <int N> class GridSearch {
6161
grid_ = grid_t(grid.data(), grid.rows(), size_);
6262
}
6363
bool stop = false; // asserted true in case of forced stop
64-
grid_.row(0).assign_to(x_old);
65-
x_new = vector_t::Constant(size_, NaN);
66-
obj_old = objective(x_old);
67-
obj_new = NaN;
64+
grid_.row(0).assign_to(x_curr);
65+
obj_curr = objective(x_curr);
6866
stop |= internals::exec_eval_hooks(*this, objective, callbacks_);
69-
values_.push_back(obj_old);
70-
if (obj_old < value_) {
71-
value_ = obj_old;
72-
optimum_ = x_old;
67+
values_.push_back(obj_curr);
68+
if (obj_curr < value_) {
69+
value_ = obj_curr;
70+
optimum_ = x_curr;
7371
}
7472
// optimize field over supplied grid
7573
for (std::size_t i = 1; i < grid_.rows() && !stop; ++i) {
76-
grid_.row(i).assign_to(x_new);
77-
obj_new = objective(x_new);
74+
grid_.row(i).assign_to(x_curr);
75+
obj_curr = objective(x_curr);
7876
stop |= internals::exec_eval_hooks(*this, objective, callbacks_);
79-
values_.push_back(obj_new);
77+
values_.push_back(obj_curr);
8078
// update minimum if better optimum found
81-
if (obj_new < value_) {
82-
value_ = obj_new;
83-
optimum_ = x_new;
79+
if (obj_curr < value_) {
80+
value_ = obj_curr;
81+
optimum_ = x_curr;
8482
}
85-
obj_old = obj_new;
86-
x_old = x_new;
83+
stop |= internals::exec_stop_if(*this, objective);
8784
}
8885
return optimum_;
8986
}

‎fdaPDE/src/optimization/lbfgs.h‎

Lines changed: 2 additions & 181 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,8 @@ template <int N> class LBFGS {
3838
int mem_size_ = 10; // number of vector used for approximating the objective hessian
3939
Eigen::Matrix<double, Dynamic, Dynamic> grad_mem_, x_mem_;
4040
public:
41+
static constexpr bool gradient_free = false;
42+
static constexpr int static_input_size = N;
4143
vector_t x_old, x_new, update, grad_old, grad_new;
4244
double h;
4345
// constructors
@@ -127,184 +129,3 @@ template <int N> class LBFGS {
127129
} // namespace fdapde
128130

129131
#endif // __FDAPDE_LBFGS_H__
130-
131-
132-
// This file is part of fdaPDE, a C++ library for physics-informed
133-
// spatial and functional data analysis.
134-
//
135-
// This program is free software: you can redistribute it and/or modify
136-
// it under the terms of the GNU General Public License as published by
137-
// the Free Software Foundation, either version 3 of the License, or
138-
// (at your option) any later version.
139-
//
140-
// This program is distributed in the hope that it will be useful,
141-
// but WITHOUT ANY WARRANTY; without even the implied warranty of
142-
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
143-
// GNU General Public License for more details.
144-
//
145-
// You should have received a copy of the GNU General Public License
146-
// along with this program. If not, see <http://www.gnu.org/licenses/>.
147-
148-
// #ifndef __FDAPDE_LBFGS_H__
149-
// #define __FDAPDE_LBFGS_H__
150-
151-
// #include "header_check.h"
152-
153-
// namespace fdapde {
154-
155-
// // Implementation of Broyden–Fletcher–Goldfarb–Shanno algorithm for unconstrained nonlinear optimization with optimized
156-
// // memory usage
157-
// template <int N, typename... Args> class LBFGS {
158-
// private:
159-
// using vector_t = std::conditional_t<N == Dynamic, Eigen::Matrix<double, Dynamic, 1>, Eigen::Matrix<double, N, 1>>;
160-
// using matrix_t =
161-
// std::conditional_t<N == Dynamic, Eigen::Matrix<double, Dynamic, Dynamic>, Eigen::Matrix<double, N, N>>;
162-
163-
// std::tuple<Args...> callbacks_;
164-
// vector_t optimum_;
165-
// double value_; // objective value at optimum
166-
// int max_iter_; // maximum number of iterations before forced stop
167-
// int n_iter_ = 0; // current iteration number
168-
// double tol_; // tolerance on error before forced stop
169-
// double step_; // update step
170-
// int mem_size_ = 10; // number of vector used for approximating the objective hessian
171-
// Eigen::Matrix<double, Dynamic, Dynamic> grad_mem_, x_mem_;
172-
// public:
173-
// vector_t x_old, x_new, update, grad_old, grad_new;
174-
// double h;
175-
// // constructors
176-
// LBFGS() = default;
177-
// LBFGS(int max_iter, double tol, double step, int mem_size)
178-
// requires(sizeof...(Args) != 0) : max_iter_(max_iter), tol_(tol), step_(step), mem_size_(mem_size) {
179-
// fdapde_assert(mem_size_ >= 0);
180-
// }
181-
// LBFGS(int max_iter, double tol, double step, int mem_size, Args&&... callbacks) :
182-
// callbacks_(std::make_tuple(std::forward<Args>(callbacks)...)),
183-
// max_iter_(max_iter),
184-
// tol_(tol),
185-
// step_(step),
186-
// mem_size_(mem_size) {
187-
// assert(mem_size_ >= 0);
188-
// }
189-
// // copy semantic
190-
// LBFGS(const LBFGS& other) :
191-
// callbacks_(other.callbacks_),
192-
// max_iter_(other.max_iter_), tol_(other.tol_), step_(other.step_), mem_size_(other.mem_size_) { }
193-
// LBFGS& operator=(const LBFGS& other) {
194-
// callbacks_ = other.callbacks_;
195-
// max_iter_ = other.max_iter_;
196-
// tol_ = other.tol_;
197-
// step_ = other.step_;
198-
// mem_size_ = other.mem_size_;
199-
// return *this;
200-
// }
201-
202-
// template <typename ObjectiveT, typename... Functor>
203-
// requires(sizeof...(Functor) < 2) && ((requires(Functor f, double value) { f(value); }) && ...)
204-
// vector_t optimize(ObjectiveT&& objective, const vector_t& x0, Functor&&... func) {
205-
// fdapde_static_assert(
206-
// std::is_same<decltype(std::declval<ObjectiveT>().operator()(vector_t())) FDAPDE_COMMA double>::value,
207-
// INVALID_CALL_TO_OPTIMIZE__OBJECTIVE_FUNCTOR_NOT_ACCEPTING_VECTORTYPE);
208-
// bool stop = false; // asserted true in case of forced stop
209-
// double error = 0;
210-
// double gamma = 1.0;
211-
// auto grad = objective.gradient();
212-
// n_iter_ = 0;
213-
// h = step_;
214-
// x_old = x0, x_new = x0;
215-
// vector_t zero;
216-
// if constexpr (N == Dynamic) {
217-
// zero = vector_t::Zero(x0.rows());
218-
// } else {
219-
// zero = vector_t::Zero();
220-
// }
221-
// update = zero;
222-
// grad_old = grad(x_old);
223-
// if (grad_old.isApprox(zero)) { // already at stationary point
224-
// optimum_ = x_old;
225-
// value_ = objective(optimum_);
226-
// if constexpr (sizeof...(Functor) == 1) { (func(value_), ...); }
227-
// return optimum_;
228-
// }
229-
// error = grad_old.norm();
230-
// x_mem_.resize(x0.rows(), mem_size_);
231-
// grad_mem_.resize(x0.rows(), mem_size_);
232-
233-
// while (n_iter_ < max_iter_ && error > tol_ && !stop) {
234-
// // compute update direction
235-
// vector_t q = grad_old;
236-
// int current_mem = n_iter_ < mem_size_ ? n_iter_ : mem_size_;
237-
// std::vector<double> alpha(current_mem, 0);
238-
// for (int i = 0; std::cmp_less(i, current_mem); ++i) {
239-
// int k = (n_iter_ + mem_size_ - i - 1) % mem_size_;
240-
// alpha[i] = x_mem_.col(k).dot(q) / grad_mem_.col(k).dot(x_mem_.col(k));
241-
// q -= alpha[i] * grad_mem_.col(k);
242-
// std::cout << "aggiorno q" << std::endl;
243-
// }
244-
// // H_0^k = I (initial guess of the inverse hessian)
245-
// std::cout << "q: " << q.transpose() << std::endl;
246-
// std::cout << "gamma: " << gamma << std::endl;
247-
248-
// update = -gamma * q;
249-
250-
// for (int i = current_mem - 1; i >= 0; --i) {
251-
// int k = (n_iter_ + mem_size_ - i - 1) % mem_size_;
252-
// double beta = grad_mem_.col(k).dot(update) / grad_mem_.col(k).dot(x_mem_.col(k));
253-
// update -= x_mem_.col(k) * (alpha[i] + beta);
254-
// }
255-
256-
// std::cout << update << std::endl;
257-
// std::cout << "----" << std::endl;
258-
259-
260-
// stop |= internals::exec_adapt_hooks(*this, objective, callbacks_);
261-
// // update along descent direction
262-
// std::cout << "h: " << h << std::endl;
263-
// std::cout << "x_old: " << x_old.transpose() << std::endl;
264-
// x_new = x_old + h * update;
265-
// grad_new = grad(x_new);
266-
267-
268-
// std::cout << "grad_new: " << grad_new.transpose() << std::endl;
269-
// std::cout << "x_new: " << x_new.transpose() << std::endl;
270-
271-
272-
// if (grad_new.isApprox(zero)) { // already at stationary point
273-
// optimum_ = x_old;
274-
// value_ = objective(optimum_);
275-
// if constexpr (sizeof...(Functor) == 1) { (func(value_), ...); }
276-
// return optimum_;
277-
// }
278-
// // mem update
279-
// // update inverse Hessian approximation
280-
// int col_idx = n_iter_ % mem_size_;
281-
// grad_mem_.col(col_idx) = grad_new - grad_old;
282-
// x_mem_.col(col_idx) = x_new - x_old;
283-
// gamma = x_mem_.col(col_idx).dot(grad_mem_.col(col_idx)) / grad_mem_.col(col_idx).norm();
284-
285-
// std::cout << "q: " << q.transpose() << std::endl;
286-
// std::cout << "gamma: " << gamma << std::endl;
287-
288-
// // prepare next iteration
289-
// if constexpr (sizeof...(Functor) == 1) { (func(objective(x_old)), ...); }
290-
// error = grad_new.norm();
291-
// // stop |=
292-
// // (execute_post_update_step(*this, objective, callbacks_) || execute_stopping_criterion(*this, objective));
293-
// x_old = x_new;
294-
// grad_old = grad_new;
295-
// ++n_iter_;
296-
// }
297-
// optimum_ = x_old;
298-
// value_ = objective(optimum_);
299-
// if constexpr (sizeof...(Functor) == 1) { (func(value_), ...); }
300-
// return optimum_;
301-
// }
302-
// // observers
303-
// vector_t optimum() const { return optimum_; }
304-
// double value() const { return value_; }
305-
// int n_iter() const { return n_iter_; }
306-
// };
307-
308-
// } // namespace fdapde
309-
310-
// #endif // __FDAPDE_LBFGS_H__

0 commit comments

Comments
 (0)