123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313 |
- [/ Copyright Matthew Pulver 2018 - 2019.
- // Distributed under the Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt or copy at
- // https://www.boost.org/LICENSE_1_0.txt)]
- [section:autodiff Automatic Differentiation]
- [template autodiff_equation[name] '''<inlinemediaobject><imageobject><imagedata fileref="../equations/autodiff/'''[name]'''"></imagedata></imageobject></inlinemediaobject>''']
- [h1:synopsis Synopsis]
- #include <boost/math/differentiation/autodiff.hpp>
-
- namespace boost {
- namespace math {
- namespace differentiation {
-
- // Function returning a single variable of differentiation. Recommended: Use auto for type.
- template <typename RealType, size_t Order, size_t... Orders>
- autodiff_fvar<RealType, Order, Orders...> make_fvar(RealType const& ca);
-
- // Function returning multiple independent variables of differentiation in a std::tuple.
- template<typename RealType, size_t... Orders, typename... RealTypes>
- auto make_ftuple(RealTypes const&... ca);
-
- // Type of combined autodiff types. Recommended: Use auto for return type (C++14).
- template <typename RealType, typename... RealTypes>
- using promote = typename detail::promote_args_n<RealType, RealTypes...>::type;
-
- namespace detail {
-
- // Single autodiff variable. Use make_fvar() or make_ftuple() to instantiate.
- template <typename RealType, size_t Order>
- class fvar {
- public:
- // Query return value of function to get the derivatives.
- template <typename... Orders>
- get_type_at<RealType, sizeof...(Orders) - 1> derivative(Orders... orders) const;
-
- // All of the arithmetic and comparison operators are overloaded.
- template <typename RealType2, size_t Order2>
- fvar& operator+=(fvar<RealType2, Order2> const&);
-
- fvar& operator+=(root_type const&);
-
- // ...
- };
-
- // Standard math functions are overloaded and called via argument-dependent lookup (ADL).
- template <typename RealType, size_t Order>
- fvar<RealType, Order> floor(fvar<RealType, Order> const&);
-
- template <typename RealType, size_t Order>
- fvar<RealType, Order> exp(fvar<RealType, Order> const&);
-
- // ...
-
- } // namespace detail
-
- } // namespace differentiation
- } // namespace math
- } // namespace boost
- [h1:description Description]
- Autodiff is a header-only C++ library that facilitates the [@https://en.wikipedia.org/wiki/Automatic_differentiation
- automatic differentiation] (forward mode) of mathematical functions of single and multiple variables.
- This implementation is based upon the [@https://en.wikipedia.org/wiki/Taylor_series Taylor series] expansion of
- an analytic function /f/ at the point ['x[sub 0]]:
- [/ Thanks to http://www.tlhiv.org/ltxpreview/ for LaTeX-to-SVG conversions. ]
- [/ \Large\begin{align*}
- f(x_0+\varepsilon) &= f(x_0) + f'(x_0)\varepsilon + \frac{f''(x_0)}{2!}\varepsilon^2 + \frac{f'''(x_0)}{3!}\varepsilon^3 + \cdots \\
- &= \sum_{n=0}^N\frac{f^{(n)}(x_0)}{n!}\varepsilon^n + O\left(\varepsilon^{N+1}\right).
- \end{align*} ]
- [:[:[autodiff_equation taylor_series.svg]]]
- The essential idea of autodiff is the substitution of numbers with polynomials in the evaluation of ['f(x[sub 0])].
- By substituting the number ['x[sub 0]] with the first-order polynomial ['x[sub 0]+\u03b5], and using the same
- algorithm to compute ['f(x[sub 0]+\u03b5)], the resulting polynomial in ['\u03b5] contains the function's derivatives
- ['f'(x[sub 0])], ['f''(x[sub 0])], ['f\'\'\'(x[sub 0])], ... within the coefficients. Each coefficient is equal to the
- derivative of its respective order, divided by the factorial of the order.
- In greater detail, assume one is interested in calculating the first /N/ derivatives of /f/ at ['x[sub 0]]. Without
- loss of precision to the calculation of the derivatives, all terms ['O(\u03b5[super N+1])] that include powers
- of ['\u03b5] greater than /N/ can be discarded. (This is due to the fact that each term in a polynomial depends
- only upon equal and lower-order terms under arithmetic operations.) Under these truncation rules, /f/ provides a
- polynomial-to-polynomial transformation:
- [/ \Large$$f \qquad : \qquad x_0+\varepsilon \qquad \mapsto \qquad
- \sum_{n=0}^Ny_n\varepsilon^n=\sum_{n=0}^N\frac{f^{(n)}(x_0)}{n!}\varepsilon^n.$$ ]
- [:[:[autodiff_equation polynomial_transform.svg]]]
- C++'s ability to overload operators and functions allows for the creation of a class `fvar` ([_f]orward-mode autodiff
- [_var]iable) that represents polynomials in ['\u03b5]. Thus the same algorithm /f/ that calculates the numeric value
- of ['y[sub 0]=f(x[sub 0])], when written to accept and return variables of a generic (template) type, is also used
- to calculate the polynomial ['\u03a3[sub n]y[sub n]\u03b5[super n]=f(x[sub 0]+\u03b5)]. The derivatives
- ['f[super (n)](x[sub 0])] are then found from the product of the respective factorial ['n!] and coefficient
- ['y[sub n]]:
- [/ \Large$$\frac{d^nf}{dx^n}(x_0)=n!y_n.$$ ]
- [:[:[autodiff_equation derivative_formula.svg]]]
- [h1:examples Examples]
- [h2:example-single-variable Example 1: Single-variable derivatives]
- [h3 Calculate derivatives of ['f(x)=x[super 4]] at /x/=2.]
- In this example, `make_fvar<double, Order>(2.0)` instantiates the polynomial 2+['\u03b5]. The `Order=5`
- means that enough space is allocated (on the stack) to hold a polynomial of up to degree 5 during the proceeding
- computation.
- Internally, this is modeled by a `std::array<double,6>` whose elements `{2, 1, 0, 0, 0, 0}` correspond
- to the 6 coefficients of the polynomial upon initialization. Its fourth power, at the end of the computation, is
- a polynomial with coefficients `y = {16, 32, 24, 8, 1, 0}`. The derivatives are obtained using the formula
- ['f[super (n)](2)=n!*y[n]].
- #include <boost/math/differentiation/autodiff.hpp>
- #include <iostream>
-
- template <typename T>
- T fourth_power(T const& x) {
- T x4 = x * x; // retval in operator*() uses x4's memory via NRVO.
- x4 *= x4; // No copies of x4 are made within operator*=() even when squaring.
- return x4; // x4 uses y's memory in main() via NRVO.
- }
-
- int main() {
- using namespace boost::math::differentiation;
-
- constexpr unsigned Order = 5; // Highest order derivative to be calculated.
- auto const x = make_fvar<double, Order>(2.0); // Find derivatives at x=2.
- auto const y = fourth_power(x);
- for (unsigned i = 0; i <= Order; ++i)
- std::cout << "y.derivative(" << i << ") = " << y.derivative(i) << std::endl;
- return 0;
- }
- /*
- Output:
- y.derivative(0) = 16
- y.derivative(1) = 32
- y.derivative(2) = 48
- y.derivative(3) = 48
- y.derivative(4) = 24
- y.derivative(5) = 0
- */
- The above calculates
- [/ \Large\begin{alignat*}{3}
- {\tt y.derivative(0)} &=& f(2) =&& \left.x^4\right|_{x=2} &= 16\\
- {\tt y.derivative(1)} &=& f'(2) =&& \left.4\cdot x^3\right|_{x=2} &= 32\\
- {\tt y.derivative(2)} &=& f''(2) =&& \left.4\cdot 3\cdot x^2\right|_{x=2} &= 48\\
- {\tt y.derivative(3)} &=& f'''(2) =&& \left.4\cdot 3\cdot2\cdot x\right|_{x=2} &= 48\\
- {\tt y.derivative(4)} &=& f^{(4)}(2) =&& 4\cdot 3\cdot2\cdot1 &= 24\\
- {\tt y.derivative(5)} &=& f^{(5)}(2) =&& 0 &
- \end{alignat*} ]
- [:[:[autodiff_equation example1.svg]]]
- [h2:example-multiprecision
- Example 2: Multi-variable mixed partial derivatives with multi-precision data type]
- [/ \Large$\frac{\partial^{12}f}{\partial w^{3}\partial x^{2}\partial y^{4}\partial z^{3}}(11,12,13,14)$]
- [/ \Large$f(w,x,y,z)=\exp\left(w\sin\left(\frac{x\log(y)}{z}\right)+\sqrt{\frac{wz}{xy}}\right)+\frac{w^2}{\tan(z)}$]
- [h3 Calculate [autodiff_equation mixed12.svg] with a precision of about 50 decimal digits,
- where [autodiff_equation example2f.svg].]
- In this example, `make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14)` returns a `std::tuple` of 4
- independent `fvar` variables, with values of 11, 12, 13, and 14, for which the maximum order derivative to
- be calculated for each are 3, 2, 4, 3, respectively. The order of the variables is important, as it is the same
- order used when calling `v.derivative(Nw, Nx, Ny, Nz)` in the example below.
- #include <boost/math/differentiation/autodiff.hpp>
- #include <boost/multiprecision/cpp_bin_float.hpp>
- #include <iostream>
-
- using namespace boost::math::differentiation;
-
- template <typename W, typename X, typename Y, typename Z>
- promote<W, X, Y, Z> f(const W& w, const X& x, const Y& y, const Z& z) {
- using namespace std;
- return exp(w * sin(x * log(y) / z) + sqrt(w * z / (x * y))) + w * w / tan(z);
- }
-
- int main() {
- using float50 = boost::multiprecision::cpp_bin_float_50;
-
- constexpr unsigned Nw = 3; // Max order of derivative to calculate for w
- constexpr unsigned Nx = 2; // Max order of derivative to calculate for x
- constexpr unsigned Ny = 4; // Max order of derivative to calculate for y
- constexpr unsigned Nz = 3; // Max order of derivative to calculate for z
- // Declare 4 independent variables together into a std::tuple.
- auto const variables = make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14);
- auto const& w = std::get<0>(variables); // Up to Nw derivatives at w=11
- auto const& x = std::get<1>(variables); // Up to Nx derivatives at x=12
- auto const& y = std::get<2>(variables); // Up to Ny derivatives at y=13
- auto const& z = std::get<3>(variables); // Up to Nz derivatives at z=14
- auto const v = f(w, x, y, z);
- // Calculated from Mathematica symbolic differentiation.
- float50 const answer("1976.319600747797717779881875290418720908121189218755");
- std::cout << std::setprecision(std::numeric_limits<float50>::digits10)
- << "mathematica : " << answer << '\n'
- << "autodiff : " << v.derivative(Nw, Nx, Ny, Nz) << '\n'
- << std::setprecision(3)
- << "relative error: " << (v.derivative(Nw, Nx, Ny, Nz) / answer - 1) << '\n';
- return 0;
- }
- /*
- Output:
- mathematica : 1976.3196007477977177798818752904187209081211892188
- autodiff : 1976.3196007477977177798818752904187209081211892188
- relative error: 2.67e-50
- */
- [h2:example-black_scholes
- Example 3: Black-Scholes Option Pricing with Greeks Automatically Calculated]
- [h3 Calculate greeks directly from the Black-Scholes pricing function.]
- Below is the standard Black-Scholes pricing function written as a function template, where the price, volatility
- (sigma), time to expiration (tau) and interest rate are template parameters. This means that any greek based on
- these 4 variables can be calculated using autodiff. The below example calculates delta and gamma where the variable
- of differentiation is only the price. For examples of more exotic greeks, see `example/black_scholes.cpp`.
- ```
- #include <boost/math/differentiation/autodiff.hpp>
- #include <iostream>
- using namespace boost::math::constants;
- using namespace boost::math::differentiation;
- // Equations and function/variable names are from
- // https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
- // Standard normal cumulative distribution function
- template <typename X>
- X Phi(X const& x) {
- return 0.5 * erfc(-one_div_root_two<X>() * x);
- }
- enum class CP { call, put };
- // Assume zero annual dividend yield (q=0).
- template <typename Price, typename Sigma, typename Tau, typename Rate>
- promote<Price, Sigma, Tau, Rate> black_scholes_option_price(CP cp,
- double K,
- Price const& S,
- Sigma const& sigma,
- Tau const& tau,
- Rate const& r) {
- using namespace std;
- auto const d1 = (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
- auto const d2 = (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
- switch (cp) {
- case CP::call:
- return S * Phi(d1) - exp(-r * tau) * K * Phi(d2);
- case CP::put:
- return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1);
- }
- }
- int main() {
- double const K = 100.0; // Strike price.
- auto const S = make_fvar<double, 2>(105); // Stock price.
- double const sigma = 5; // Volatility.
- double const tau = 30.0 / 365; // Time to expiration in years. (30 days).
- double const r = 1.25 / 100; // Interest rate.
- auto const call_price = black_scholes_option_price(CP::call, K, S, sigma, tau, r);
- auto const put_price = black_scholes_option_price(CP::put, K, S, sigma, tau, r);
- std::cout << "black-scholes call price = " << call_price.derivative(0) << '\n'
- << "black-scholes put price = " << put_price.derivative(0) << '\n'
- << "call delta = " << call_price.derivative(1) << '\n'
- << "put delta = " << put_price.derivative(1) << '\n'
- << "call gamma = " << call_price.derivative(2) << '\n'
- << "put gamma = " << put_price.derivative(2) << '\n';
- return 0;
- }
- /*
- Output:
- black-scholes call price = 56.5136
- black-scholes put price = 51.4109
- call delta = 0.773818
- put delta = -0.226182
- call gamma = 0.00199852
- put gamma = 0.00199852
- */
- ```
- [h1 Advantages of Automatic Differentiation]
- The above examples illustrate some of the advantages of using autodiff:
- * Elimination of code redundancy. The existence of /N/ separate functions to calculate derivatives is a form
- of code redundancy, with all the liabilities that come with it:
- * Changes to one function require /N/ additional changes to other functions. In the 3rd example above,
- consider how much larger and inter-dependent the above code base would be if a separate function were
- written for [@https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
- each Greek] value.
- * Dependencies upon a derivative function for a different purpose will break when changes are made to
- the original function. What doesn't need to exist cannot break.
- * Code bloat, reducing conceptual integrity. Control over the evolution of code is easier/safer when
- the code base is smaller and able to be intuitively grasped.
- * Accuracy of derivatives over finite difference methods. Single-iteration finite difference methods always include
- a ['\u0394x] free variable that must be carefully chosen for each application. If ['\u0394x] is too small, then
- numerical errors become large. If ['\u0394x] is too large, then mathematical errors become large. With autodiff,
- there are no free variables to set and the accuracy of the answer is generally superior to finite difference
- methods even with the best choice of ['\u0394x].
- [h1 Manual]
- Additional details are in the [@../differentiation/autodiff.pdf autodiff manual].
- [endsect]
|