Optimist  0.0.0
A C++ library for optimization
Loading...
Searching...
No Matches
FiniteDifferences.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2 * Copyright (c) 2025, Davide Stocco, Mattia Piazza and Enrico Bertolazzi. *
3 * *
4 * The Optimist project is distributed under the BSD 2-Clause License. *
5 * *
6 * Davide Stocco Mattia Piazza Enrico Bertolazzi *
7 * University of Trento University of Trento University of Trento *
8 * davide.stocco@unitn.it mattia.piazza@unitn.it enrico.bertolazzi@unitn.it *
9\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10
11#pragma once
12
13#ifndef OPTIMIST_FINITE_DIFFERENCES_HH
14#define OPTIMIST_FINITE_DIFFERENCES_HH
15
16#include "Optimist.hh"
17
18namespace Optimist {
19
20 namespace FiniteDifferences {
21
26 template<typename Real>
27 class Epsilon {
28 OPTIMIST_BASIC_CONSTANTS(Real)
29 Real m_epsilon_1{std::sqrt(EPSILON)};
30 Real m_epsilon_2{std::pow(EPSILON, 0.75)};
31 Real m_epsilon_3{std::pow(EPSILON, 0.25)};
32 public:
33
35 Real epsilon_1(Real const v) const {return (std::abs(v) + 1.0)*this->m_epsilon_1;}
36
38 Real epsilon_2(Real const v) const {return (std::abs(v) + 1.0)*this->m_epsilon_2;}
39
41 Real epsilon_3(Real const v) const {return (std::abs(v) + 1.0)*this->m_epsilon_3;}
42 };
43
44
55 template<typename Real>
56 inline Real SideFiniteDifferences(Real const fun_0, Real const fun_1, Real const fun_2, Real h_1,
57 Real h_2)
58 {
59 Real d_fun_1{fun_1 - fun_0};
60 Real d_fun_2{fun_2 - fun_0};
61 Real d_h{h_1 - h_2};
62 return ((h_1/h_2)*d_fun_2 - (h_2/h_1)*d_fun_1) / d_h;
63 }
64
74 template<typename Real>
75 inline Real CenteredFiniteDifferences(Real const fun_l, Real const fun_c, Real const fun_r, Real h)
76 {
77 constexpr Real EPSILON{std::numeric_limits<Real>::epsilon()};
78 Real diff_r{(fun_r-fun_c)/h};
79 Real diff_l{(fun_c-fun_l)/h};
80 Real weight_r{std::abs(diff_r) + EPSILON};
81 Real weight_l{std::abs(diff_l) + EPSILON};
82 Real weight_max{std::max(weight_r, weight_l)};
83 weight_r = std::sqrt(weight_r/weight_max);
84 weight_l = std::sqrt(weight_l/weight_max);
85 return (diff_r*weight_l + diff_l*weight_r) / (weight_r+weight_l);
86 }
87
88 /*\
89 | ____ _
90 | | _ \ _ _ _ __ __ _ _ __ ___ (_) ___
91 | | | | | | | | '_ \ / _` | '_ ` _ \| |/ __|
92 | | |_| | |_| | | | | (_| | | | | | | | (__
93 | |____/ \__, |_| |_|\__,_|_| |_| |_|_|\___|
94 | |___/
95 \*/
96
98 template<typename Real>
99 using VectorXd = Eigen::Vector<Real, Eigen::Dynamic>;
100
102 template<typename Real>
103 using MatrixXd = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
104
115 template<typename Real>
116 inline VectorXd<Real> SideFiniteDifferences(VectorXd<Real> const & fun_0, VectorXd<Real> const & fun_1,
117 VectorXd<Real> const & fun_2, Real h_1, Real h_2)
118 {
119 Real d_h{h_1 - h_2};
120 Real t_1{-(h_2/h_1) / d_h};
121 Real t_2{(h_1/h_2) / d_h};
122 return t_1*(fun_1 - fun_0) + t_2*(fun_2 - fun_0);
123 }
124
134 template<typename Real>
135 inline VectorXd<Real> CenteredFiniteDifferences(VectorXd<Real> const & fun_l, VectorXd<Real> const & fun_c,
136 VectorXd<Real> const & fun_r, Real h)
137 {
138 constexpr Real EPSILON{std::numeric_limits<Real>::epsilon()};
139 VectorXd<Real> diff_r((fun_r - fun_c) / h);
140 VectorXd<Real> diff_l((fun_c - fun_l) / h);
141 VectorXd<Real> weight_r(diff_r.array().abs() + EPSILON);
142 VectorXd<Real> weight_l(diff_l.array().abs() + EPSILON);
143 VectorXd<Real> weight_max(weight_r.array().max(weight_l.array()));
144 weight_r = (weight_r.array()/weight_max.array()).sqrt();
145 weight_l = (weight_l.array()/weight_max.array()).sqrt();
146 return (diff_r.array()*weight_l.array() + diff_l.array()*weight_r.array()) / (weight_r.array()+weight_l.array());
147 }
148
158 template <typename Function, typename Real>
159 bool Gradient(VectorXd<Real> const & x, Function fun, VectorXd<Real> & grad)
160 {
161 Epsilon<Real> eps;
162 Real v_c{0.0};
163 if (!fun(x, v_c) || !std::isfinite(v_c)) {return false;}
164 grad.resize(x.size());
165 for (Integer i{0}; i < x.size(); ++i) {
166 VectorXd<Real> v_x(x);
167 Real tmp{x[i]};
168 Real h_1{eps.epsilon_1(tmp)};
169 Real h_2{eps.epsilon_2(tmp)};
170 v_x[i] = tmp + h_1; Real v_r{0.0}; bool is_finite_r{fun(v_x, v_r) && std::isfinite(v_r)};
171 v_x[i] = tmp - h_1; Real v_l{0.0}; bool is_finite_l{fun(v_x, v_l) && std::isfinite(v_l)};
172 Integer ic{(is_finite_r && is_finite_l) ? 0 : (is_finite_r ? 1 : (is_finite_l ? -1 : -2))};
173 switch (ic) {
174 case 0:
175 grad[i] = CenteredFiniteDifferences<Real>(v_l, v_c, v_r, h_1);
176 break;
177 case 1: {
178 v_x[i] = tmp + h_2; Real v_rr{0.0}; bool is_finite_rr{fun(v_x, v_rr) && std::isfinite(v_rr)};
179 grad[i] = is_finite_rr ? SideFiniteDifferences<Real>(v_c, v_r, v_rr, h_1, h_2) : (v_r-v_c)/h_1;
180 break;
181 }
182 case -1: {
183 v_x[i] = tmp - h_2; Real v_ll{0.0}; bool is_finite_ll{fun(v_x, v_ll) && std::isfinite(v_ll)};
184 grad[i] = is_finite_ll ? SideFiniteDifferences<Real>(v_c, v_l, v_ll, -h_1, -h_2) : (v_c-v_l)/h_1;
185 break;
186 }
187 case -2: {
188 grad[i] = 0.0;
189 return false;
190 }
191 }
192 }
193 return grad.allFinite();
194 }
195
205 template <typename Function, typename Real>
206 bool Jacobian(VectorXd<Real> const & x, Function fun, MatrixXd<Real> & jac)
207 {
208 Epsilon<Real> eps;
209 VectorXd<Real> v_c;
210 if (!fun(x, v_c) || !v_c.allFinite()) {return false;}
211 Integer dim_x{static_cast<Integer>(x.size())};
212 Integer dim_f{static_cast<Integer>(v_c.size())};
213 jac.resize(dim_f, dim_x);
214 for (Integer j{0}; j < dim_x; ++j) {
215 VectorXd<Real> v_x(x);
216 Real tmp{x[j]};
217 Real h_1{eps.epsilon_1(tmp)};
218 Real h_2{eps.epsilon_2(tmp)};
219 v_x[j] = tmp + h_1; VectorXd<Real> v_r; bool is_finite_r{fun(v_x, v_r) && v_r.allFinite()};
220 v_x[j] = tmp - h_1; VectorXd<Real> v_l; bool is_finite_l{fun(v_x, v_l) && v_l.allFinite()};
221 Integer ic{(is_finite_r && is_finite_l) ? 0 : (is_finite_r ? 1 : (is_finite_l ? -1 : -2))};
222 switch (ic) {
223 case 0:
224 jac.col(j) = CenteredFiniteDifferences<Real>(v_l, v_c, v_r, h_1);
225 break;
226 case 1: {
227 v_x[j] = tmp + h_2; VectorXd<Real> v_rr; bool is_finite_rr{fun(v_x, v_rr) && v_rr.allFinite()};
228 jac.col(j) = is_finite_rr ? SideFiniteDifferences<Real>(v_c, v_r, v_rr, h_1, h_2) : (v_r-v_c)/h_1;
229 break;
230 }
231 case -1: {
232 v_x[j] = tmp - h_2; VectorXd<Real> v_ll; bool is_finite_ll{fun(v_x, v_ll) && v_ll.allFinite()};
233 jac.col(j) = is_finite_ll ? SideFiniteDifferences<Real>(v_c, v_l, v_ll, -h_1, -h_2) : (v_c-v_l)/h_1;
234 break;
235 }
236 case -2: {
237 jac.col(j).setZero();
238 return false;
239 }
240 }
241 }
242 return jac.allFinite();
243 }
244
254 template <typename Function, typename Real>
255 bool Hessian(VectorXd<Real> const & x, Function fun, MatrixXd<Real> & hes)
256 {
257 Epsilon<Real> eps;
258 Integer dim_x{x.size()};
259 hes.resize(dim_x, dim_x);
260 Real fc{0.0};
261 if (!fun(x, fc) || !std::isfinite(fc)) {return false;}
262 for (Integer j{0}; j < dim_x; ++j) {
263 Real tmp_j{x[j]};
264 Real h_j{eps.epsilon_3(tmp_j)};
265 VectorXd<Real> v_x(x);
266 v_x[j] = tmp_j + h_j; Real fp{0.0}; if (!fun(v_x, fp) || !std::isfinite(fp)) {return false;}
267 v_x[j] = tmp_j - h_j; Real fm{0.0}; if (!fun(v_x, fm) || !std::isfinite(fm)) {return false;}
268 hes(j, j) = ((fp + fm) - 2.0 * fc) / (h_j * h_j);
269 for (Integer i{j + 1}; i < dim_x; ++i) {
270 Real tmp_i{x[i]};
271 Real h_i{eps.epsilon_3(tmp_i)};
272 v_x[i] = tmp_i + h_i; v_x[j] = tmp_j + h_j; Real fpp{0.0}; if (!fun(v_x, fpp) || !std::isfinite(fpp)) {return false;}
273 v_x[i] = tmp_i - h_i; Real fmp{0.0}; if (!fun(v_x, fmp) || !std::isfinite(fmp)) {return false;}
274 v_x[j] = tmp_j - h_j; Real fmm{0.0}; if (!fun(v_x, fmm) || !std::isfinite(fmm)) {return false;}
275 v_x[i] = tmp_i + h_i; Real fpm{0.0}; if (!fun(v_x, fpm) || !std::isfinite(fpm)) {return false;}
276 Real h_ij{4.0 * h_i * h_j};
277 Real value{((fpp + fmm) - (fpm + fmp)) / h_ij};
278 hes(j, i) = hes(i, j) = value;
279 v_x[i] = tmp_i;
280 }
281 v_x[j] = tmp_j;
282 }
283 return hes.allFinite();
284 }
285
286 /*\
287 | _____ _ _
288 | | ___(_)_ _____ __| |
289 | | |_ | \ \/ / _ \/ _` |
290 | | _| | |> < __/ (_| |
291 | |_| |_/_/\_\___|\__,_|
292 |
293 \*/
294
305 template<typename Vector>
306 inline Vector SideFiniteDifferences(Vector const & fun_0, Vector const & fun_1,
307 Vector const & fun_2, typename Vector::Scalar h_1, typename Vector::Scalar h_2)
308 {
309 using Real = typename Vector::Scalar;
310 Real d_h{h_1 - h_2};
311 Real t_1{-(h_2/h_1) / d_h};
312 Real t_2{(h_1/h_2) / d_h};
313 return t_1*(fun_1 - fun_0) + t_2*(fun_2 - fun_0);
314 }
315
325 template<typename Vector>
326 inline Vector CenteredFiniteDifferences(Vector const & fun_l, Vector const & fun_c, Vector const & fun_r,
327 typename Vector::Scalar h)
328 {
329 using Real = typename Vector::Scalar;
330 constexpr Real EPSILON{std::numeric_limits<Real>::epsilon()};
331 Vector diff_r((fun_r - fun_c) / h);
332 Vector diff_l((fun_c - fun_l) / h);
333 Vector weight_r(diff_r.array().abs() + EPSILON);
334 Vector weight_l(diff_l.array().abs() + EPSILON);
335 Vector weight_max(weight_r.array().max(weight_l.array()));
336 weight_r = (weight_r.array()/weight_max.array()).sqrt();
337 weight_l = (weight_l.array()/weight_max.array()).sqrt();
338 return (diff_r.array()*weight_l.array() + diff_l.array()*weight_r.array()) / (weight_r.array()+weight_l.array());
339 }
340
350 template <typename Function, typename Vector>
351 bool Gradient(Vector const & x, Function fun, Vector & grad)
352 {
353 using Real = typename Vector::Scalar;
354 Epsilon<Real> eps;
355 Real v_c{0.0};
356 if (!fun(x, v_c) || !std::isfinite(v_c)) {return false;}
357 constexpr Integer dim_x{Vector::RowsAtCompileTime};
358 grad.setZero();
359 for (Integer i = 0; i < dim_x; ++i) {
360 Vector v_x(x);
361 Real tmp{x(i)};
362 Real h_1{eps.epsilon_1(tmp)};
363 Real h_2{eps.epsilon_2(tmp)};
364 v_x(i) = tmp + h_1; Real v_r{0.0}; bool is_finite_r = fun(v_x, v_r) && std::isfinite(v_r);
365 v_x(i) = tmp - h_1; Real v_l{0.0}; bool is_finite_l = fun(v_x, v_l) && std::isfinite(v_l);
366 Integer ic{(is_finite_r && is_finite_l) ? 0 : (is_finite_r ? 1 : (is_finite_l ? -1 : -2))};
367 switch (ic) {
368 case 0:
369 grad(i) = CenteredFiniteDifferences<Real>(v_l, v_c, v_r, h_1);
370 break;
371 case 1: {
372 v_x(i) = tmp + h_2; Real v_rr{0.0}; bool is_finite_rr = fun(v_x, v_rr) && std::isfinite(v_rr);
373 grad(i) = is_finite_rr ? SideFiniteDifferences<Real>(v_c, v_r, v_rr, h_1, h_2) : (v_r-v_c)/h_1;
374 break;
375 }
376 case -1: {
377 v_x(i) = tmp - h_2; Real v_ll{0.0}; bool is_finite_ll = fun(v_x, v_ll) && std::isfinite(v_ll);
378 grad(i) = is_finite_ll ? SideFiniteDifferences<Real>(v_c, v_l, v_ll, -h_1, -h_2) : (v_c-v_l)/h_1;
379 break;
380 }
381 case -2: {
382 grad(i) = 0.0;
383 return false;
384 }
385 }
386 }
387 return grad.allFinite();
388 }
389
400 template <typename Function, typename Vector, typename Matrix>
401 bool Jacobian(Vector const & x, Function fun, Matrix & jac)
402 {
403 using Real = typename Vector::Scalar;
404 Epsilon<Real> eps;
405 Vector v_c;
406 if (!fun(x, v_c) || !v_c.allFinite()) {return false;}
407 constexpr Integer dim_x{Vector::RowsAtCompileTime};
408 jac.setZero();
409 for (Integer j{0}; j < dim_x; ++j) {
410 Vector v_x(x);
411 Real tmp{x(j)};
412 Real h_1{eps.epsilon_1(tmp)};
413 Real h_2{eps.epsilon_2(tmp)};
414 v_x(j) = tmp + h_1; Vector v_r; bool is_finite_r = fun(v_x, v_r) && v_r.allFinite();
415 v_x(j) = tmp - h_1; Vector v_l; bool is_finite_l = fun(v_x, v_l) && v_l.allFinite();
416 Integer ic{(is_finite_r && is_finite_l) ? 0 : (is_finite_r ? 1 : (is_finite_l ? -1 : -2))};
417 switch (ic) {
418 case 0:
419 jac.col(j) = CenteredFiniteDifferences<Real>(v_l, v_c, v_r, h_1);
420 break;
421 case 1: {
422 v_x(j) = tmp + h_2; Vector v_rr; bool is_finite_rr = fun(v_x, v_rr) && v_rr.allFinite();
423 jac.col(j) = is_finite_rr ? SideFiniteDifferences<Real>(v_c, v_r, v_rr, h_1, h_2) : (v_r-v_c)/h_1;
424 break;
425 }
426 case -1: {
427 v_x(j) = tmp - h_2; Vector v_ll; bool is_finite_ll = fun(v_x, v_ll) && v_ll.allFinite();
428 jac.col(j) = is_finite_ll ? SideFiniteDifferences<Real>(v_c, v_l, v_ll, -h_1, -h_2) : (v_c-v_l)/h_1;
429 break;
430 }
431 case -2: {
432 jac.col(j).setZero();
433 return false;
434 }
435 }
436 }
437 return jac.allFinite();
438 }
439
450 template <typename Function, typename Vector, typename Matrix>
451 bool Hessian(Vector const & x, Function fun, Matrix & hes)
452 {
453 using Real = typename Vector::Scalar;
454 Epsilon<Real> eps;
455 constexpr Integer dim_x = Vector::RowsAtCompileTime;
456 hes.setZero();
457 Real fc{0.0};
458 if (!fun(x, fc) || !std::isfinite(fc)) {return false;}
459 for (Integer j{0}; j < dim_x; ++j) {
460 Real tmp_j{x(j)};
461 Real h_j{eps.epsilon_3(tmp_j)};
462 Vector v_x = x;
463 v_x(j) = tmp_j + h_j; Real fp{0.0}; if (!fun(v_x, fp) || !std::isfinite(fp)) {return false;}
464 v_x(j) = tmp_j - h_j; Real fm{0.0}; if (!fun(v_x, fm) || !std::isfinite(fm)) {return false;}
465 hes(j, j) = ((fp + fm) - 2.0 * fc) / (h_j * h_j);
466 for (Integer i{j + 1}; i < dim_x; ++i) {
467 Real tmp_i{x(i)};
468 Real h_i{eps.epsilon_3(tmp_i)};
469 v_x(i) = tmp_i + h_i; v_x(j) = tmp_j + h_j; Real fpp{0.0}; if (!fun(v_x, fpp) || !std::isfinite(fpp)) {return false;}
470 v_x(i) = tmp_i - h_i; Real fmp{0.0}; if (!fun(v_x, fmp) || !std::isfinite(fmp)) {return false;}
471 v_x(j) = tmp_j - h_j; Real fmm{0.0}; if (!fun(v_x, fmm) || !std::isfinite(fmm)) {return false;}
472 v_x(i) = tmp_i + h_i; Real fpm{0.0}; if (!fun(v_x, fpm) || !std::isfinite(fpm)) {return false;}
473 Real h_ij{4.0 * h_i * h_j};
474 Real value{((fpp + fmm) - (fpm + fmp)) / h_ij};
475 hes(j, i) = hes(i, j) = value;
476 v_x(i) = tmp_i;
477 }
478 v_x(j) = tmp_j;
479 }
480 return hes.allFinite();
481 }
482
483 } // namespace FiniteDifferences
484
485} // namespace Optimist
486
487#endif // OPTIMIST_FINITE_DIFFERENCES_HH
OPTIMIST_BASIC_CONSTANTS
#define OPTIMIST_BASIC_CONSTANTS(Real)
Definition Optimist.hh:71
Optimist::FiniteDifferences::Epsilon
Definition FiniteDifferences.hh:27
Optimist::FiniteDifferences::Epsilon::m_epsilon_3
Real m_epsilon_3
Definition FiniteDifferences.hh:31
Optimist::FiniteDifferences::Epsilon::m_epsilon_1
Real m_epsilon_1
Definition FiniteDifferences.hh:29
Optimist::FiniteDifferences::Epsilon::epsilon_3
Real epsilon_3(Real const v) const
Definition FiniteDifferences.hh:41
Optimist::FiniteDifferences::Epsilon::epsilon_2
Real epsilon_2(Real const v) const
Definition FiniteDifferences.hh:38
Optimist::FiniteDifferences::Epsilon::m_epsilon_2
Real m_epsilon_2
Definition FiniteDifferences.hh:30
Optimist::FiniteDifferences::Epsilon::epsilon_1
Real epsilon_1(Real const v) const
Definition FiniteDifferences.hh:35
Optimist::Function
Class container for the vector-valued function.
Definition Function.hh:184
Optimist::FiniteDifferences
Definition FiniteDifferences.hh:20
Optimist::FiniteDifferences::Gradient
bool Gradient(VectorXd< Real > const &x, Function fun, VectorXd< Real > &grad)
Definition FiniteDifferences.hh:159
Optimist::FiniteDifferences::VectorXd
Eigen::Vector< Real, Eigen::Dynamic > VectorXd
Definition FiniteDifferences.hh:99
Optimist::FiniteDifferences::Hessian
bool Hessian(VectorXd< Real > const &x, Function fun, MatrixXd< Real > &hes)
Definition FiniteDifferences.hh:255
Optimist::FiniteDifferences::SideFiniteDifferences
Real SideFiniteDifferences(Real const fun_0, Real const fun_1, Real const fun_2, Real h_1, Real h_2)
Definition FiniteDifferences.hh:56
Optimist::FiniteDifferences::MatrixXd
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > MatrixXd
Definition FiniteDifferences.hh:103
Optimist::FiniteDifferences::Jacobian
bool Jacobian(VectorXd< Real > const &x, Function fun, MatrixXd< Real > &jac)
Definition FiniteDifferences.hh:206
Optimist::FiniteDifferences::CenteredFiniteDifferences
Real CenteredFiniteDifferences(Real const fun_l, Real const fun_c, Real const fun_r, Real h)
Definition FiniteDifferences.hh:75
Optimist
Namespace for the Optimist library.
Definition Optimist.hh:88
Optimist::Integer
OPTIMIST_DEFAULT_INTEGER_TYPE Integer
The Integer type as used for the API.
Definition Optimist.hh:96