Tutorials¶
We present here some examples demonstrating the use of autodiff for computing different types of derivatives. We welcome any contribution towards improving and expanding this list of examples. We would also love to hear your suggestions on how to better demonstrate the capabilities of autodiff.
Forward mode¶
Derivatives of a single-variable function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
// The single-variable function for which derivatives are needed
dual f(dual x)
{
return 1 + x + x*x + 1/x + log(x);
}
int main()
{
dual x = 2.0; // the input variable x
dual u = f(x); // the output variable u
double dudx = derivative(f, wrt(x), at(x)); // evaluate the derivative du/dx
std::cout << "u = " << u << std::endl; // print the evaluated output u
std::cout << "du/dx = " << dudx << std::endl; // print the evaluated derivative du/dx
}
Derivatives of a single-variable function using a custom scalar (complex)¶
// C++ includes
#include <iostream>
#include <complex>
using namespace std;
// autodiff include
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
// Specialize isArithmetic for complex to make it compatible with dual
namespace autodiff::detail {
template<typename T>
struct ArithmeticTraits<complex<T>> : ArithmeticTraits<T> {};
} // autodiff::detail
using cxdual = Dual<complex<double>, complex<double>>;
// The single-variable function for which derivatives are needed
cxdual f(cxdual x)
{
return 1 + x + x*x + 1/x + log(x);
}
int main()
{
cxdual x = 2.0; // the input variable x
cxdual u = f(x); // the output variable u
cxdual dudx = derivative(f, wrt(x), at(x)); // evaluate the derivative du/dx
cout << "u = " << u << endl; // print the evaluated output u
cout << "du/dx = " << dudx << endl; // print the evaluated derivative du/dx
}
Derivatives of a multi-variable function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
// The multi-variable function for which derivatives are needed
dual f(dual x, dual y, dual z)
{
return 1 + x + y + z + x*y + y*z + x*z + x*y*z + exp(x/y + y/z);
}
int main()
{
dual x = 1.0;
dual y = 2.0;
dual z = 3.0;
dual u = f(x, y, z);
double dudx = derivative(f, wrt(x), at(x, y, z));
double dudy = derivative(f, wrt(y), at(x, y, z));
double dudz = derivative(f, wrt(z), at(x, y, z));
std::cout << "u = " << u << std::endl; // print the evaluated output u = f(x, y, z)
std::cout << "du/dx = " << dudx << std::endl; // print the evaluated derivative du/dx
std::cout << "du/dy = " << dudy << std::endl; // print the evaluated derivative du/dy
std::cout << "du/dz = " << dudz << std::endl; // print the evaluated derivative du/dz
}
Derivatives of a multi-variable function with parameters¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
// A type defining parameters for a function of interest
struct Params
{
dual a;
dual b;
dual c;
};
// The function that depends on parameters for which derivatives are needed
dual f(dual x, const Params& params)
{
return params.a * sin(x) + params.b * cos(x) + params.c * sin(x)*cos(x);
}
int main()
{
Params params; // initialize the parameter variables
params.a = 1.0; // the parameter a of type dual, not double!
params.b = 2.0; // the parameter b of type dual, not double!
params.c = 3.0; // the parameter c of type dual, not double!
dual x = 0.5; // the input variable x
dual u = f(x, params); // the output variable u
double dudx = derivative(f, wrt(x), at(x, params)); // evaluate the derivative du/dx
double duda = derivative(f, wrt(params.a), at(x, params)); // evaluate the derivative du/da
double dudb = derivative(f, wrt(params.b), at(x, params)); // evaluate the derivative du/db
double dudc = derivative(f, wrt(params.c), at(x, params)); // evaluate the derivative du/dc
std::cout << "u = " << u << std::endl; // print the evaluated output u
std::cout << "du/dx = " << dudx << std::endl; // print the evaluated derivative du/dx
std::cout << "du/da = " << duda << std::endl; // print the evaluated derivative du/da
std::cout << "du/db = " << dudb << std::endl; // print the evaluated derivative du/db
std::cout << "du/dc = " << dudc << std::endl; // print the evaluated derivative du/dc
}
/*-------------------------------------------------------------------------------------------------
=== Note ===
---------------------------------------------------------------------------------------------------
This example would also work if real was used instead of dual. Should you
need higher-order cross derivatives, however, e.g.,:
double d2udxda = derivative(f, wrt(x, params.a), at(x, params));
then higher-order dual types are the right choicesince real types are
optimally designed for higher-order directional derivatives.
-------------------------------------------------------------------------------------------------*/
Derivatives of a multi-variable function that also relies on analytical derivatives¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
// Define functions A, Ax, Ay using double; analytical derivatives are available.
double A(double x, double y) { return x*y; }
double Ax(double x, double y) { return y; }
double Ay(double x, double y) { return x; }
// Define functions B, Bx, By using double; analytical derivatives are available.
double B(double x, double y) { return x + y; }
double Bx(double x, double y) { return 1.0; }
double By(double x, double y) { return 1.0; }
// Wrap A into Adual function so that it can be used within autodiff-enabled codes.
dual Adual(dual const& x, dual const& y)
{
dual res = A(x.val, y.val);
if(x.grad != 0.0)
res.grad += x.grad * Ax(x.val, y.val);
if(y.grad != 0.0)
res.grad += y.grad * Ay(x.val, y.val);
return res;
}
// Wrap B into Bdual function so that it can be used within autodiff-enabled codes.
dual Bdual(dual const& x, dual const& y)
{
dual res = B(x.val, y.val);
if(x.grad != 0.0)
res.grad += x.grad * Bx(x.val, y.val);
if(y.grad != 0.0)
res.grad += y.grad * By(x.val, y.val);
return res;
}
// Define autodiff-enabled C function that relies on Adual and Bdual
dual C(dual const& x, dual const& y)
{
const auto A = Adual(x, y);
const auto B = Bdual(x, y);
return A*A + B;
}
int main()
{
dual x = 1.0;
dual y = 2.0;
auto C0 = C(x, y);
// Compute derivatives of C with respect to x and y using autodiff!
auto Cx = derivative(C, wrt(x), at(x, y));
auto Cy = derivative(C, wrt(y), at(x, y));
// Compute expected analytical derivatives of C with respect to x and y
auto x0 = x.val;
auto y0 = y.val;
auto expectedCx = 2.0*A(x0, y0)*Ax(x0, y0) + Bx(x0, y0);
auto expectedCy = 2.0*A(x0, y0)*Ay(x0, y0) + By(x0, y0);
std::cout << "C0 = " << C0 << "\n";
std::cout << "Cx(computed) = " << Cx << "\n";
std::cout << "Cx(expected) = " << expectedCx << "\n";
std::cout << "Cy(computed) = " << Cy << "\n";
std::cout << "Cy(expected) = " << expectedCy << "\n";
}
// Output:
// C0 = 7
// Cx(computed) = 5
// Cx(expected) = 5
// Cy(computed) = 9
// Cy(expected) = 9
Gradient vector of a scalar function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The scalar function for which the gradient is needed
real f(const ArrayXreal& x)
{
return (x * x.exp()).sum(); // sum([xi * exp(xi) for i = 1:5])
}
int main()
{
using Eigen::VectorXd;
ArrayXreal x(5); // the input array x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
real u; // the output scalar u = f(x) evaluated together with gradient below
VectorXd g = gradient(f, wrt(x), at(x), u); // evaluate the function value u and its gradient vector g = du/dx
std::cout << "u = " << u << std::endl; // print the evaluated output u
std::cout << "g = \n" << g << std::endl; // print the evaluated gradient vector g = du/dx
}
Gradient vector of a scalar function with parameters¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The scalar function for which the gradient is needed
real f(const ArrayXreal& x, const ArrayXreal& p, const real& q)
{
return (x * x).sum() * p.sum() * exp(q); // sum([xi * xi for i = 1:5]) * sum(p) * exp(q)
}
int main()
{
using Eigen::VectorXd;
ArrayXreal x(5); // the input vector x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
ArrayXreal p(3); // the input parameter vector p with 3 variables
p << 1, 2, 3; // p = [1, 2, 3]
real q = -2; // the input parameter q as a single variable
real u; // the output scalar u = f(x, p, q) evaluated together with gradient below
VectorXd gx = gradient(f, wrt(x), at(x, p, q), u); // evaluate the function value u and its gradient vector gx = du/dx
VectorXd gp = gradient(f, wrt(p), at(x, p, q), u); // evaluate the function value u and its gradient vector gp = du/dp
VectorXd gq = gradient(f, wrt(q), at(x, p, q), u); // evaluate the function value u and its gradient vector gq = du/dq
VectorXd gqpx = gradient(f, wrt(q, p, x), at(x, p, q), u); // evaluate the function value u and its gradient vector gqpx = [du/dq, du/dp, du/dx]
std::cout << "u = " << u << std::endl; // print the evaluated output u
std::cout << "gx = \n" << gx << std::endl; // print the evaluated gradient vector gx = du/dx
std::cout << "gp = \n" << gp << std::endl; // print the evaluated gradient vector gp = du/dp
std::cout << "gq = \n" << gq << std::endl; // print the evaluated gradient vector gq = du/dq
std::cout << "gqpx = \n" << gqpx << std::endl; // print the evaluated gradient vector gqpx = [du/dq, du/dp, du/dx]
}
/*-------------------------------------------------------------------------------------------------
=== Note ===
---------------------------------------------------------------------------------------------------
This example would also work if dual was used instead. However, if gradient,
Jacobian, and directional derivatives are all you need, then real types are your
best option. You want to use dual types when evaluating higher-order cross
derivatives, which is not supported for real types.
-------------------------------------------------------------------------------------------------*/
Jacobian matrix of a vector function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The vector function for which the Jacobian is needed
VectorXreal f(const VectorXreal& x)
{
return x * x.sum();
}
int main()
{
using Eigen::MatrixXd;
VectorXreal x(5); // the input vector x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
VectorXreal F; // the output vector F = f(x) evaluated together with Jacobian matrix below
MatrixXd J = jacobian(f, wrt(x), at(x), F); // evaluate the output vector F and the Jacobian matrix dF/dx
std::cout << "F = \n" << F << std::endl; // print the evaluated output vector F
std::cout << "J = \n" << J << std::endl; // print the evaluated Jacobian matrix dF/dx
}
Jacobian matrix of a vector function with parameters¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The vector function with parameters for which the Jacobian is needed
VectorXreal f(const VectorXreal& x, const VectorXreal& p, const real& q)
{
return x * p.sum() * exp(q);
}
int main()
{
using Eigen::MatrixXd;
VectorXreal x(5); // the input vector x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
VectorXreal p(3); // the input parameter vector p with 3 variables
p << 1, 2, 3; // p = [1, 2, 3]
real q = -2; // the input parameter q as a single variable
VectorXreal F; // the output vector F = f(x, p, q) evaluated together with Jacobian below
MatrixXd Jx = jacobian(f, wrt(x), at(x, p, q), F); // evaluate the function and the Jacobian matrix Jx = dF/dx
MatrixXd Jp = jacobian(f, wrt(p), at(x, p, q), F); // evaluate the function and the Jacobian matrix Jp = dF/dp
MatrixXd Jq = jacobian(f, wrt(q), at(x, p, q), F); // evaluate the function and the Jacobian matrix Jq = dF/dq
MatrixXd Jqpx = jacobian(f, wrt(q, p, x), at(x, p, q), F); // evaluate the function and the Jacobian matrix Jqpx = [dF/dq, dF/dp, dF/dx]
std::cout << "F = \n" << F << std::endl; // print the evaluated output vector F
std::cout << "Jx = \n" << Jx << std::endl; // print the evaluated Jacobian matrix dF/dx
std::cout << "Jp = \n" << Jp << std::endl; // print the evaluated Jacobian matrix dF/dp
std::cout << "Jq = \n" << Jq << std::endl; // print the evaluated Jacobian matrix dF/dq
std::cout << "Jqpx = \n" << Jqpx << std::endl; // print the evaluated Jacobian matrix [dF/dq, dF/dp, dF/dx]
}
Jacobian matrix of a vector function using memory maps¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The vector function for which the Jacobian is needed
VectorXreal f(const VectorXreal& x)
{
return x * x.sum();
}
int main()
{
using Eigen::Map;
using Eigen::MatrixXd;
VectorXreal x(5); // the input vector x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
double y[25]; // the output Jacobian as a flat array
VectorXreal F; // the output vector F = f(x) evaluated together with Jacobian matrix below
Map<MatrixXd> J(y, 5, 5); // the output Jacobian dF/dx mapped onto the flat array
jacobian(f, wrt(x), at(x), F, J); // evaluate the output vector F and the Jacobian matrix dF/dx
std::cout << "F = \n" << F << std::endl; // print the evaluated output vector F
std::cout << "J = \n" << J << std::endl; // print the evaluated Jacobian matrix dF/dx
std::cout << "y = "; // print the flat array
for(int i = 0 ; i < 25 ; ++i)
std::cout << y[i] << " ";
std::cout << std::endl;
}
Higher-order cross derivatives of a scalar function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
// The multi-variable function for which higher-order derivatives are needed (up to 4th order)
dual4th f(dual4th x, dual4th y, dual4th z)
{
return 1 + x + y + z + x*y + y*z + x*z + x*y*z + exp(x/y + y/z);
}
int main()
{
dual4th x = 1.0;
dual4th y = 2.0;
dual4th z = 3.0;
auto [u0, ux, uxy, uxyx, uxyxz] = derivatives(f, wrt(x, y, x, z), at(x, y, z));
std::cout << "u0 = " << u0 << std::endl; // print the evaluated value of u = f(x, y, z)
std::cout << "ux = " << ux << std::endl; // print the evaluated derivative du/dx
std::cout << "uxy = " << uxy << std::endl; // print the evaluated derivative d²u/dxdy
std::cout << "uxyx = " << uxyx << std::endl; // print the evaluated derivative d³u/dx²dy
std::cout << "uxyxz = " << uxyxz << std::endl; // print the evaluated derivative d⁴u/dx²dydz
}
/*-------------------------------------------------------------------------------------------------
=== Note ===
---------------------------------------------------------------------------------------------------
In most cases, dual can be replaced by real, as commented in other examples.
However, computing higher-order cross derivatives has definitely to be done
using higher-order dual types (e.g., dual3rd, dual4th)! This is because real
types (e.g., real2nd, real3rd, real4th) are optimally designed for computing
higher-order directional derivatives.
-------------------------------------------------------------------------------------------------*/
Higher-order directional derivatives of a scalar function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
using namespace autodiff;
// The multi-variable function for which higher-order derivatives are needed (up to 4th order)
real4th f(real4th x, real4th y, real4th z)
{
return sin(x) * cos(y) * exp(z);
}
int main()
{
real4th x = 1.0;
real4th y = 2.0;
real4th z = 3.0;
auto dfdv = derivatives(f, along(1.0, 1.0, 2.0), at(x, y, z)); // the directional derivatives of f along direction v = (1, 1, 2) at (x, y, z) = (1, 2, 3)
std::cout << "dfdv[0] = " << dfdv[0] << std::endl; // print the evaluated 0th order directional derivative of f along v (equivalent to f(x, y, z))
std::cout << "dfdv[1] = " << dfdv[1] << std::endl; // print the evaluated 1st order directional derivative of f along v
std::cout << "dfdv[2] = " << dfdv[2] << std::endl; // print the evaluated 2nd order directional derivative of f along v
std::cout << "dfdv[3] = " << dfdv[3] << std::endl; // print the evaluated 3rd order directional derivative of f along v
std::cout << "dfdv[4] = " << dfdv[4] << std::endl; // print the evaluated 4th order directional derivative of f along v
}
/*-------------------------------------------------------------------------------------------------
=== Note ===
---------------------------------------------------------------------------------------------------
This example would also work if dual was used instead of real. However, real
types are your best option for directional derivatives, as they were optimally
designed for this kind of derivatives.
-------------------------------------------------------------------------------------------------*/
Higher-order directional derivatives of a vector function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The vector function for which higher-order directional derivatives are needed (up to 4th order).
ArrayXreal4th f(const ArrayXreal4th& x, real4th p)
{
return p * x.log();
}
int main()
{
using Eigen::ArrayXd;
ArrayXreal4th x(5); // the input vector x
x << 1.0, 2.0, 3.0, 4.0, 5.0;
real4th p = 7.0; // the input parameter p = 1
ArrayXd vx(5); // the direction vx in the direction vector v = (vx, vp)
vx << 1.0, 1.0, 1.0, 1.0, 1.0;
double vp = 1.0; // the direction vp in the direction vector v = (vx, vp)
auto dfdv = derivatives(f, along(vx, vp), at(x, p)); // the directional derivatives of f along direction v = (vx, vp) at (x, p)
std::cout << std::scientific << std::showpos;
std::cout << "Directional derivatives of f along v = (vx, vp) up to 4th order:" << std::endl;
std::cout << "dfdv[0] = " << dfdv[0].transpose() << std::endl; // print the evaluated 0th order directional derivative of f along v (equivalent to f(x, p))
std::cout << "dfdv[1] = " << dfdv[1].transpose() << std::endl; // print the evaluated 1st order directional derivative of f along v
std::cout << "dfdv[2] = " << dfdv[2].transpose() << std::endl; // print the evaluated 2nd order directional derivative of f along v
std::cout << "dfdv[3] = " << dfdv[3].transpose() << std::endl; // print the evaluated 3rd order directional derivative of f along v
std::cout << "dfdv[4] = " << dfdv[4].transpose() << std::endl; // print the evaluated 4th order directional derivative of f along v
std::cout << std::endl;
double t = 0.1; // the step length along direction v = (vx, vp) used to compute 4th order Taylor estimate of f
ArrayXreal4th u = f(x + t * vx, p + t * vp); // start from (x, p), walk a step length t = 0.1 along direction v = (vx, vp) and evaluate f at this current point
ArrayXd utaylor = dfdv[0] + t*dfdv[1] + (t*t/2.0)*dfdv[2] + (t*t*t/6.0)*dfdv[3] + (t*t*t*t/24.0)*dfdv[4]; // evaluate the 4th order Taylor estimate of f along direction v = (vx, vp) at a step length of t = 0.1
std::cout << "Comparison between exact evaluation and 4th order Taylor estimate:" << std::endl;
std::cout << "u(exact) = " << u.transpose() << std::endl;
std::cout << "u(taylor) = " << utaylor.transpose() << std::endl;
}
/*-------------------------------------------------------------------------------------------------
=== Output ===
---------------------------------------------------------------------------------------------------
Directional derivatives of f along v = (vx, vp) up to 4th order:
dfdv[0] = +0.000000e+00 +4.852030e+00 +7.690286e+00 +9.704061e+00 +1.126607e+01
dfdv[1] = +7.000000e+00 +4.193147e+00 +3.431946e+00 +3.136294e+00 +3.009438e+00
dfdv[2] = -5.000000e+00 -7.500000e-01 -1.111111e-01 +6.250000e-02 +1.200000e-01
dfdv[3] = +1.100000e+01 +1.000000e+00 +1.851852e-01 +3.125000e-02 -8.000000e-03
dfdv[4] = -3.400000e+01 -1.625000e+00 -2.222222e-01 -3.906250e-02 -3.200000e-03
Comparison between exact evaluation and 4th order Taylor estimate:
u(exact) = +6.767023e-01 +5.267755e+00 +8.032955e+00 +1.001801e+01 +1.156761e+01
u(taylor) = +6.766917e-01 +5.267755e+00 +8.032955e+00 +1.001801e+01 +1.156761e+01
-------------------------------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------------------------------
=== Note ===
---------------------------------------------------------------------------------------------------
This example would also work if dual was used instead of real. However, real
types are your best option for directional derivatives, as they were optimally
designed for this kind of derivatives.
-------------------------------------------------------------------------------------------------*/
Taylor series of a scalar function along a direction¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The scalar function for which a 4th order directional Taylor series will be computed.
real4th f(const real4th& x, const real4th& y, const real4th& z)
{
return sin(x * y) * cos(x * z) * exp(z);
}
int main()
{
real4th x = 1.0; // the input vector x
real4th y = 2.0; // the input vector y
real4th z = 3.0; // the input vector z
auto g = taylorseries(f, along(1, 1, 2), at(x, y, z)); // the function g(t) as a 4th order Taylor approximation of f(x + t, y + t, z + 2t)
double t = 0.1; // the step length used to evaluate g(t), the Taylor approximation of f(x + t, y + t, z + 2t)
real4th u = f(x + t, y + t, z + 2*t); // the exact value of f(x + t, y + t, z + 2t)
double utaylor = g(t); // the 4th order Taylor estimate of f(x + t, y + t, z + 2t)
std::cout << std::fixed;
std::cout << "Comparison between exact evaluation and 4th order Taylor estimate of f(x + t, y + t, z + 2t):" << std::endl;
std::cout << "u(exact) = " << u << std::endl;
std::cout << "u(taylor) = " << utaylor << std::endl;
}
/*-------------------------------------------------------------------------------------------------
=== Output ===
---------------------------------------------------------------------------------------------------
Comparison between exact evaluation and 4th order Taylor estimate of f(x + t, y + t, z + 2t):
u(exact) = -16.847071
u(taylor) = -16.793986
-------------------------------------------------------------------------------------------------*/
Taylor series of a vector function along a direction¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
using namespace autodiff;
// The vector function for which a 4th order directional Taylor series will be computed.
ArrayXreal4th f(const ArrayXreal4th& x)
{
return x.sin() / x;
}
int main()
{
using Eigen::ArrayXd;
ArrayXreal4th x(5); // the input vector x
x << 1.0, 2.0, 3.0, 4.0, 5.0;
ArrayXd v(5); // the direction vector v
v << 1.0, 1.0, 1.0, 1.0, 1.0;
auto g = taylorseries(f, along(v), at(x)); // the function g(t) as a 4th order Taylor approximation of f(x + t·v)
double t = 0.1; // the step length used to evaluate g(t), the Taylor approximation of f(x + t·v)
ArrayXreal4th u = f(x + t * v); // the exact value of f(x + t·v)
ArrayXd utaylor = g(t); // the 4th order Taylor estimate of f(x + t·v)
std::cout << std::fixed;
std::cout << "Comparison between exact evaluation and 4th order Taylor estimate of f(x + t·v):" << std::endl;
std::cout << "u(exact) = " << u.transpose() << std::endl;
std::cout << "u(taylor) = " << utaylor.transpose() << std::endl;
}
/*-------------------------------------------------------------------------------------------------
=== Output ===
---------------------------------------------------------------------------------------------------
Comparison between exact evaluation and 4th order Taylor estimate of f(x + t·v):
u(exact) = 0.810189 0.411052 0.013413 -0.199580 -0.181532
u(taylor) = 0.810189 0.411052 0.013413 -0.199580 -0.181532
-------------------------------------------------------------------------------------------------*/
Reverse mode¶
Single-variable function¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
using namespace autodiff;
// The single-variable function for which derivatives are needed
var f(var x)
{
return 1 + x + x*x + 1/x + log(x);
}
int main()
{
var x = 2.0; // the input variable x
var u = f(x); // the output variable u
auto [ux] = derivatives(u, wrt(x)); // evaluate the derivative of u with respect to x
cout << "u = " << u << endl; // print the evaluated output variable u
cout << "ux = " << ux << endl; // print the evaluated derivative ux
}
Multi-variable function¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
using namespace autodiff;
// The multi-variable function for which derivatives are needed
var f(var x, var y, var z)
{
return 1 + x + y + z + x*y + y*z + x*z + x*y*z + exp(x/y + y/z);
}
int main()
{
var x = 1.0; // the input variable x
var y = 2.0; // the input variable y
var z = 3.0; // the input variable z
var u = f(x, y, z); // the output variable u
auto [ux, uy, uz] = derivatives(u, wrt(x, y, z)); // evaluate the derivatives of u with respect to x, y, z
cout << "u = " << u << endl; // print the evaluated output u
cout << "ux = " << ux << endl; // print the evaluated derivative ux
cout << "uy = " << uy << endl; // print the evaluated derivative uy
cout << "uz = " << uz << endl; // print the evaluated derivative uz
}
Multi-variable function with conditional¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
using namespace autodiff;
// A two-variable piecewise function for which derivatives are needed
var f(var x, var y) { return condition(x < y, x * y, x * x); }
int main()
{
var x = 1.0; // the input variable x
var y = 2.0; // the input variable y
var u = f(x, y); // the output variable u
auto [ux, uy] = derivatives(u, wrt(x, y)); // evaluate the derivatives of u
cout << "x = " << x << ", y = " << y << endl;
cout << "u = " << u << endl; // x = 1, y = 2, so x < y, so x * y = 2
cout << "ux = " << ux << endl; // d/dx x * y = y = 2
cout << "uy = " << uy << endl; // d/dy x * y = x = 1
x.update(3.0); // Change the value of x in the expression tree
u.update(); // Update the expression tree in a sweep
auto [ux2, uy2] = derivatives(u, wrt(x, y)); // re-evaluate the derivatives
cout << "x = " << x << ", y = " << y << endl;
cout << "u = " << u << endl; // Now x > y, so x * x = 9
cout << "ux = " << ux2 << endl; // d/dx x * x = 2x = 6
cout << "uy = " << uy2 << endl; // d/dy x * x = 0
// condition-associated functions
cout << "min(x, y) = " << min(x, y) << endl;
cout << "max(x, y) = " << max(x, y) << endl;
cout << "sgn(x) = " << sgn(x) << endl;
}
Multi-variable function with parameters¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
using namespace autodiff;
// A type defining parameters for a function of interest
struct Params
{
var a;
var b;
var c;
};
// The function that depends on parameters for which derivatives are needed
var f(var x, const Params& params)
{
return params.a * sin(x) + params.b * cos(x) + params.c * sin(x)*cos(x);
}
int main()
{
Params params; // initialize the parameter variables
params.a = 1.0; // the parameter a of type var, not double!
params.b = 2.0; // the parameter b of type var, not double!
params.c = 3.0; // the parameter c of type var, not double!
var x = 0.5; // the input variable x
var u = f(x, params); // the output variable u
auto [ux, ua, ub, uc] = derivatives(u, wrt(x, params.a, params.b, params.c)); // evaluate the derivatives of u with respect to x and parameters a, b, c
cout << "u = " << u << endl; // print the evaluated output u
cout << "ux = " << ux << endl; // print the evaluated derivative du/dx
cout << "ua = " << ua << endl; // print the evaluated derivative du/da
cout << "ub = " << ub << endl; // print the evaluated derivative du/db
cout << "uc = " << uc << endl; // print the evaluated derivative du/dc
}
Gradient of a scalar function¶
// C++ includes
#include <iostream>
// autodiff include
#include <autodiff/reverse/var.hpp>
#include <autodiff/reverse/var/eigen.hpp>
using namespace autodiff;
// The scalar function for which the gradient is needed
var f(const ArrayXvar& x)
{
return sqrt((x * x).sum()); // sqrt(sum([xi * xi for i = 1:5]))
}
int main()
{
using Eigen::VectorXd;
VectorXvar x(5); // the input vector x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
var y = f(x); // the output variable y
VectorXd dydx = gradient(y, x); // evaluate the gradient vector dy/dx
std::cout << "y = " << y << std::endl; // print the evaluated output y
std::cout << "dy/dx = \n" << dydx << std::endl; // print the evaluated gradient vector dy/dx
}
Hessian of a scalar function¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
#include <autodiff/reverse/var/eigen.hpp>
using namespace autodiff;
// The scalar function for which the gradient is needed
var f(const ArrayXvar& x)
{
return sqrt((x * x).sum()); // sqrt(sum([xi * xi for i = 1:5]))
}
int main()
{
VectorXvar x(5); // the input vector x with 5 variables
x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
var u = f(x); // the output variable u
Eigen::VectorXd g; // the gradient vector to be computed in method `hessian`
Eigen::MatrixXd H = hessian(u, x, g); // evaluate the Hessian matrix H and the gradient vector g of u
cout << "u = " << u << endl; // print the evaluated output variable u
cout << "g = \n" << g << endl; // print the evaluated gradient vector of u
cout << "H = \n" << H << endl; // print the evaluated Hessian matrix of u
}
Higher-order derivatives of a single-variable function¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
using namespace autodiff;
int main()
{
var x = 0.5; // the input variable x
var u = sin(x) * cos(x); // the output variable u
auto [ux] = derivativesx(u, wrt(x)); // evaluate the first order derivatives of u
auto [uxx] = derivativesx(ux, wrt(x)); // evaluate the second order derivatives of ux
cout << "u = " << u << endl; // print the evaluated output variable u
cout << "ux(autodiff) = " << ux << endl; // print the evaluated first order derivative ux
cout << "ux(exact) = " << 1 - 2*sin(x)*sin(x) << endl; // print the exact first order derivative ux
cout << "uxx(autodiff) = " << uxx << endl; // print the evaluated second order derivative uxx
cout << "uxx(exact) = " << -4*cos(x)*sin(x) << endl; // print the exact second order derivative uxx
}
/*===============================================================================
Output:
=================================================================================
u = 0.420735
ux(autodiff) = 0.540302
ux(exact) = 0.540302
uxx(autodiff) = -1.68294
uxx(exact) = -1.68294
===============================================================================*/
Higher-order derivatives of a multi-variable function¶
// C++ includes
#include <iostream>
using namespace std;
// autodiff include
#include <autodiff/reverse/var.hpp>
using namespace autodiff;
int main()
{
var x = 1.0; // the input variable x
var y = 0.5; // the input variable y
var z = 2.0; // the input variable z
var u = x * log(y) * exp(z); // the output variable u
auto [ux, uy, uz] = derivativesx(u, wrt(x, y, z)); // evaluate the derivatives of u with respect to x, y, z.
auto [uxx, uxy, uxz] = derivativesx(ux, wrt(x, y, z)); // evaluate the derivatives of ux with respect to x, y, z.
auto [uyx, uyy, uyz] = derivativesx(uy, wrt(x, y, z)); // evaluate the derivatives of uy with respect to x, y, z.
auto [uzx, uzy, uzz] = derivativesx(uz, wrt(x, y, z)); // evaluate the derivatives of uz with respect to x, y, z.
cout << "u = " << u << endl; // print the evaluated output variable u
cout << "ux = " << ux << endl; // print the evaluated first order derivative ux
cout << "uy = " << uy << endl; // print the evaluated first order derivative uy
cout << "uz = " << uz << endl; // print the evaluated first order derivative uz
cout << "uxx = " << uxx << endl; // print the evaluated second order derivative uxx
cout << "uxy = " << uxy << endl; // print the evaluated second order derivative uxy
cout << "uxz = " << uxz << endl; // print the evaluated second order derivative uxz
cout << "uyx = " << uyx << endl; // print the evaluated second order derivative uyx
cout << "uyy = " << uyy << endl; // print the evaluated second order derivative uyy
cout << "uyz = " << uyz << endl; // print the evaluated second order derivative uyz
cout << "uzx = " << uzx << endl; // print the evaluated second order derivative uzx
cout << "uzy = " << uzy << endl; // print the evaluated second order derivative uzy
cout << "uzz = " << uzz << endl; // print the evaluated second order derivative uzz
}
Integration with CMake-based projects¶
Integrating autodiff in a CMake-based project is very simple as shown next.
Let's assume our CMake-based project consists of two files: main.cpp and
CMakeLists.txt, whose contents are shown below:
main.cpp
#include <iostream>
#include <autodiff/forward/dual.hpp>
using namespace autodiff;
dual f(dual x)
{
return 1 + x + x*x + 1/x + log(x);
}
int main()
{
dual x = 1.0;
dual u = f(x);
double dudx = derivative(f, wrt(x), at(x));
std::cout << "u = " << u << std::endl;
std::cout << "du/dx = " << dudx << std::endl;
}
CMakeLists.txt
cmake_minimum_required(VERSION 3.16)
project(app)
find_package(autodiff)
add_executable(app main.cpp)
target_link_libraries(app autodiff::autodiff)
In the CMakeLists.txt file, note the use of the command:
find_package(autodiff)
to find the header files of the autodiff library, and the command:
target_link_libraries(app autodiff::autodiff)
app against the autodiff library
(autodiff::autodiff) using CMake's modern target-based design.
To build the application, do:
mkdir build && cd build
cmake .. -DCMAKE_PREFIX_PATH=/path/to/autodiff/install/dir
make
Attention
If autodiff has been installed system-wide, then the CMake argument
CMAKE_PREFIX_PATH should not be needed. Otherwise, you will need to specify
where autodiff is installed in your machine. For example:
cmake .. -DCMAKE_PREFIX_PATH=$HOME/local
assuming directory $HOME/local is where autodiff was installed to, which
should then contain the following directory:
$HOME/local/include/autodiff/
where the autodiff header files are located.
To execute the application, do:
./app