// TestExp5.cpp
//
// dx/dt = - grad(U(x)) , grad == gradient
// Transform f(x) = 0 to a least-squares problem, U = f*f
//
// Tip iceberg: minimise a function using an ODE solver for a gradient system.
//
// compute z = exp(5); f(z) = z - exp(5); argmin f(z)*f(z)
//
// It's a sledgehammer but the method can be applied to constrained optimisation for ODE systems.
// To be discussed elsewhere. It can be applied in all sorts of places.
// Current ML wisdom uses Gradient Descent (GD) method, which is Euler's method, and much weaker in many
// ways to the approach taken here. It's much more robust and general than GD.
//
// Exercise:
// 1. generalise to ODE systems.
// 2. generalise to constrained optimisation and penalty method.
// 3. methods to compute gradient in n dimensions.
// 4. Use in an ANN instead of flaky GD.
//
// The world is continuous, but the mind is discrete. David Mumford
// Natura non facit saltum (Nature does nothing in jumps)
//
// My take is that Gradient Descent is a fabrication; ODE gradient system are closer
// to the physical world.
//
// DJD
//
// Using ODEs for optimisation.
//
// http://www.unige.ch/~hairer/preprints/gradientflow.pdf
// https://authors.library.caltech.edu/26703/2/postscript.pdf
// https://blogs.mathworks.com/cleve/2013/10/14/complex-step-differentiation/
// https://www.dam.brown.edu/people/geman/Homepage/Image%20processing,%20image%20analysis,%20Markov%20random%20fields,%20and%20MCMC/Diffusions%20and%20optimization.pdf
// https://www.cs.ubc.ca/~ascher/papers/adhs.pdf
//
// Nesterov's method
// https://epubs.siam.org/doi/epdf/10.1137/21M1390037
//
//
// For training courses, see
//
// https://www.datasim.nl/onlinecourses/97/distance-learning-ordinary-and-partial-differential-equations
//
// https://www.datasim.nl/application/files/9015/4809/1157/DL_Ordinary_and_Partial_Differential_Equations.pdf
//
// https://www.datasim.nl/onlinecourses
//
// See also my book on ODE/PDE/FDM in computational finance
//
// https://www.wiley.com/en-us/Numerical+Methods+in+Computational+Finance:+A+Partial+Differential+Equation+(PDE+FDM)+Approach-p-9781119719670
//
// (C) Datasim Education BV 2018-2023
//
//
#include <iostream>
#include <vector>
#include <cmath>
#include <complex>
// https://www.boost.org/doc/libs/1_82_0/libs/numeric/odeint/doc/html/index.html
#include <boost/numeric/odeint.hpp>
// Preliminary notation
using value_type = double;
using state_type = std::vector<value_type>;
// Using C++11 functions
template <typename T>
using FunctionType = std::function<T(const T& arg)>;
using CFunctionType = FunctionType<std::complex<value_type>>;
template <typename T>
FunctionType<T> operator * (FunctionType<T>& f, FunctionType<T>& g)
{ // Multiplication (higher-order function)
return [=](T x)
{
return f(x)*g(x);
};
}
// Complex Step Method for gradients (others: divided difference, AAD, analytic,..)
value_type CSM(const CFunctionType& f, value_type x, value_type h)
{ // df/dx at x using the Complex step method (Trapp/Squire)
std::complex<value_type> z(x, h); // x + ih, i = sqrt(-1)
return std::imag(f(z)) / h;
}
std::complex<double> ObjFunc(const std::complex<value_type>& z)
{ // My hard-coded specific function
// compute z = exp(5); f(z) = z - exp(5); argmin f(z)*f(z)
const std::complex<double> a(std::exp(5.0), 0.0); // => 148.413
// Original function f(x) = 0 (Newton)
CFunctionType f = [&](const std::complex<double>& z) { return z - a; };
CFunctionType f2 = f * f;
// Least squares function
return f2(z);
}
// Free function to model RHS in dy/dt = RHS(t,y)
void Ode(const state_type &x, state_type &dxdt, const value_type t)
{
//dxdt[0] = 2.0 * (x[0] - std::exp(5.0)); // Exact derivative
const double h = 1.0e-156; // Small h but no overflow!
dxdt[0] = CSM(ObjFunc, x[0], h);
// Transform semi-infinite time interval (0, infinity) to (0,1)
// Then we look at asymptotic behaviour..
double tau = (1.0 - t) * (1.0 - t);
dxdt[0] = -dxdt[0]/tau;
}
void print(std::string name, std::size_t steps, value_type v)
{
std::cout << "Number of steps " << name << std::setprecision(16) << steps << "approximate: " << v << '\n';
}
int main()
{
namespace Bode = boost::numeric::odeint;
// Initial condition
value_type L = 0.0;
value_type T = 0.99; // HEURISTIC simulates T = infinity
value_type initialCondition = 1110.548; // Convergence independent of IC?
value_type dt = 1.0e-5;
{
// Cash Karp, middle-of-road
state_type x{ initialCondition };
std::size_t steps = Bode::integrate(Ode, x, L, T, dt);
print("Cash Karp ", steps, x[0]);
}
{
// Bulirsch-Stoer method, high-class solver
state_type x{ initialCondition};
Bode::bulirsch_stoer<state_type, value_type> bsStepper; // O(variable), controlled stepper
std::size_t steps = Bode::integrate_const(bsStepper, Ode, x, L, T, dt);
print("Bulirsch-Stoer ", steps, x[0]);
}
return 0;
}
