Agreed! You should be less defensive when someone finds a higher-performance solution to a cherry-picked function.Give this discussion a chance to run its course. Going all defensive is a bit silly.
And less cherry picking, por favor
Agreed! You should be less defensive when someone finds a higher-performance solution to a cherry-picked function.Give this discussion a chance to run its course. Going all defensive is a bit silly.
And less cherry picking, por favor
Are you that 'someone'? As I said, using Mathematica is cheating in this case. So tell us what's the best solution? It's also OK to say you don't know.Agreed! You should be less defensive when someone finds a higher-performance solution to a cherry-picked function.Give this discussion a chance to run its course. Going all defensive is a bit silly.
And less cherry picking, por favor
So now accessing Wolfram Alpha is out? (But accessing old math papers is still permitted?) Moving the goal posts once again?Are you that 'someone'? As I said, using Mathematica is cheating in this case. So tell us what's the best solution? It's also OK to say you don't know.Agreed! You should be less defensive when someone finds a higher-performance solution to a cherry-picked function.Give this discussion a chance to run its course. Going all defensive is a bit silly.
And less cherry picking, por favor
Cool.The use of the complex-step in calculating the derivative of an analytic function was introduced by Lyness & Moler.I complexify [$]f(z)[$] where [$]z = x + ih[$] and h = 0.001, e.g.
Then compute the imaginary part of [$]f(z)/h[$] and you are done.(*)
For x = 1.87 I get 1.3684..... e+289 for both exact and this method.(12 digits accuracy)
No subtraction cancellation as with FD classic.
(*) Squire and Trapp 1998.
Numerical Differentiation of Analytics Functions, SIAM Vol 4, N2, 1967 available here: link to the pdf paper
PS: In a previous life, I used the complex-step gradients to compute the sensittivity derivatives of the aerodynamic cost function.
+ DiffSharp?For the AD method for greeks, have to look to Mike Giles work.
There's a lot to be said for this approach for this computationally intensive problem.I stay with universal languages, e.g C++/Fortran
// TestComplexStep.cpp
//
// Complex-step method to compute approximate derivatives.
// Example is scalar-valued function of a scalar argument.
//
// https://pdfs.semanticscholar.org/3de7/e8ae217a4214507b9abdac66503f057aaae9.pdf
//
// http://mdolab.engin.umich.edu/sites/default/files/Martins2003CSD.pdf
//
// (C) Datasim Education BV 2018
//
#include <functional>
#include <complex>
#include <iostream>
#include <iomanip>
#include <cmath>
// Notation and function spaces
using value_type = double;
template <typename T>
using FunctionType = std::function < T(const T& c)>;
using CFunctionType = FunctionType<std::complex<value_type>>;
// Test case from Squire&Trapp 1998
template <typename T> T func(const T& t)
{
T n1 = std::exp(t);
T d1 = std::sin(t);
T d2 = std::cos(t);
return n1 / (d1*d1*d1 + d2*d2*d2);
}
template <typename T> T func2(const T& t)
{ // Derivative of e^t, sanity check
return std::exp(std::pow(t,1));
// return std::exp(std::pow(t, 5));
}
value_type Derivative(const CFunctionType& f, value_type x, value_type h)
{ // df/dx at x using tbe Complex step method
std::complex<value_type> z(x, h); // x + ih, i = sqrt(-1)
return std::imag(f(z)) / h;
}
int main()
{
// Squire Trapp
double x = 1.5; double h = 0.1;
do
{
std::cout << std::setprecision(12) << Derivative(func<std::complex<value_type>>, x, h) << '\n';
h *= 0.1;
} while (h > 1.0e-300);
// Exponential function (101 sanity check)
x = 5.0;
h = 1.0e-10;
std::cout << "Exponential 1: " << std::setprecision(12) << Derivative(func2<std::complex<value_type>>, x, h) << '\n';
return 0;
}
ublas::matrix<W> d = ...
std::size_t n = ...
for (std::size_t v = 0; v < n; ++v)
{
d(v, v) = 0.0;
}
for (std::size_t k = 0; k < n; ++k)
{
for (std::size_t i = 0; i < n; ++i)
{
for (std::size_t j = 0; j < n; ++j)
{
d(i, j) = std::min<W>(d(i, j), d(i, k) + d(k,j));
}
}
}
elised?
It is much more "complex" to use it with a computer code. Suppose you have a C++ program which output some "cost" function and you want compute the derivative of this cost function to do some optimisation problem with respect to the cost function.Computing derivatives without the pain of 1) differentiation, 2) catastrophic round-off errors. Sample code.
The scalar (almost a 1-liner). I leave the vector case as an exercise.
Code: Select all// TestComplexStep.cpp // // Complex-step method to compute approximate derivatives. // Example is scalar-valued function of a scalar argument. // // https://pdfs.semanticscholar.org/3de7/e8ae217a4214507b9abdac66503f057aaae9.pdf // // http://mdolab.engin.umich.edu/sites/default/files/Martins2003CSD.pdf // // (C) Datasim Education BV 2018 // #include <functional> #include <complex> #include <iostream> #include <iomanip> #include <cmath> // Notation and function spaces using value_type = double; template <typename T> using FunctionType = std::function < T(const T& c)>; using CFunctionType = FunctionType<std::complex<value_type>>; // Test case from Squire&Trapp 1998 template <typename T> T func(const T& t) { T n1 = std::exp(t); T d1 = std::sin(t); T d2 = std::cos(t); return n1 / (d1*d1*d1 + d2*d2*d2); } template <typename T> T func2(const T& t) { // Derivative of e^t, sanity check return std::exp(std::pow(t,1)); // return std::exp(std::pow(t, 5)); } value_type Derivative(const CFunctionType& f, value_type x, value_type h) { // df/dx at x using tbe Complex step method std::complex<value_type> z(x, h); // x + ih, i = sqrt(-1) return std::imag(f(z)) / h; } int main() { // Squire Trapp double x = 1.5; double h = 0.1; do { std::cout << std::setprecision(12) << Derivative(func<std::complex<value_type>>, x, h) << '\n'; h *= 0.1; } while (h > 1.0e-300); // Exponential function (101 sanity check) x = 5.0; h = 1.0e-10; std::cout << "Exponential 1: " << std::setprecision(12) << Derivative(func2<std::complex<value_type>>, x, h) << '\n'; return 0; }
/***********START***************/
void quad(double a, double b, double c, double & x1, double & x2 ){
double r = sqrt(b*b - 4 * a*c);
x1 = (-b + r) / (2 * a);
x2 = (-b - r) / (2 * a);
}
double quad2(double a, double b, double c, double x ) {
return a*x*x +b*x+c;
}
void QuadraticLossPrecision(){
double x1(0), x2(0);
quad(1, 1e8, 1, x1, x2);
fout << "\n\tx1:" << x1 << "\n\tx2:" << x2;
fout << "\n\tx1: " << x1 << "\n\tquad2(x1):" << quad2(1, 1e8, 1, x1);
}
/***********END QuadraticLossPrecision*****************/
°ExSan++ High Performance C++ Computing V.6.00.L°
NOT VALID FOR DISTRIBUTION ---BE---careful
Thu Nov 29 08:36:52 2018
[email protected] https://twitter.com/#!/ExSan_com
JOB: quad_V.6.00.L__3652
x1:-0.000000007450580596923828125000
x2:-100000000.000000000000000000000000000000
x1:-0.000000007450580596923828125000
quad2(x1): 0.254941940307617187500000000000
ENDS quad_V.6.00.L__3652 Elapsed Time: 0.008000000000000000166533453694 sec
NOT VALID FOR DISTRIBUTION ---BE---careful
EXIT FROM EXSAN