Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
Automatic and symbolic differentiation

A group dedicated to the implementation of functions and classes that relate to automatic and symbolic differentiation. More...

Namespaces

namespace  Differentiation
 
namespace  Differentiation::AD
 
namespace  Differentiation::SD
 

Detailed Description

A group dedicated to the implementation of functions and classes that relate to automatic and symbolic differentiation.

Below we provide a very brief introduction as to what automatic and symbolic differentiation are, what variations of these computational/numerical schemes exist, and how they are integrated within deal.II's framework. The purpose of all of these schemes is to automatically compute the derivative of functions, or approximations of it, in cases where one does not want to compute them by hand. Common examples are situations in the finite element context is where one wants to solve a nonlinear problem that is given by requiring that some residual \(F(u,\nabla u)=0\) where \(F\) is a complicated function that needs to be differentiated to apply Newton's method; and situations where one is given a parameter dependent problem \({\cal A}(q,u,\nabla u) = f\) and wants to form derivatives with regards to the parameters \(q\), for example to optimize an output functional with regards to \(q\), or for a sensitivity analysis with regards to \(q\). One should think of \(q\) as design parameters: say, the width or shape of a wing, the stiffness coefficients of a material chosen to build an object, the power sent to a device, the chemical composition of the gases sent to a burner. In all of these cases, one should think of \(F\) and \(\cal A\) as complicated and cumbersome to differentiate – at least when doing it by hand. A relatively simple case of a nonlinear problem that already highlights the tedium of computing derivatives by hand is shown in step-15. However, in reality, one might, for example, think about problems such as chemically reactive flows where the fluid equations have coefficients such as the density and viscosity that depend strongly and nonlinearly on the chemical composition, temperature, and pressure of the fluid at each point; and where the chemical species react with each other based on reaction coefficients that also depend nonlinearly and in complicated ways on the chemical composition, temperature, and pressure. In many cases, the exact formulas for all of these coefficients can take several lines to write out, may include exponentials and (harmonic or geometric) averages of several nonlinear terms, and/or may contain table lookup of and interpolation between data points. Just getting these terms right is difficult enough; computing derivatives of these terms is impractical in most applications and, in reality, impossible to get right. Higher derivatives are even more impossible to do without computer aid. Automatic or symbolic differentiation is a way out of this: One only has to implement the function that computes these coefficients in terms of their inputs only once, and gets the (correct!) derivatives without further coding effort (though at a non-negligible computational cost either at run time, compile time, or both).

Automatic differentiation

Automatic differentiation (commonly also referred to as algorithmic differentiation), is a numerical method that can be used to "automatically" compute the first, and perhaps higher-order, derivatives of function(s) with respect to one or more input variables. Although this comes at a certain computational cost, the benefits to using such a tool may be significant. When used correctly the derivatives of often complicated functions can be computed to a very high accuracy. Although the exact accuracy achievable by these frameworks largely depends on their underlying mathematical formulation, some implementations compute with a precision on the order of machine accuracy. Note that this is different to classical numerical differentiation (using, for example, a finite difference approximation of a function by evaluating it at different points), which has an accuracy that depends on both the perturbation size as well as the chosen finite-difference scheme; the error of these methods is measurably larger than well-formulated automatic differentiation approaches.

Three practical examples of auto-differentiation use within a finite-element context would then be

There are quite a number of implementations for auto-differentiable numbers. They primarily fall into two broad categories, namely source code transformation and operator overloading. The first method generates new, compilable code based on some input function that, when executed, returns the derivatives of the input function. The second exploits the capability of C++ operator definitions to be overloaded for custom class types. Therefore a class that represents such an auto-differentiable number can, following each mathematical operation performed on or with it, in principle evaluate and keep track of its value as well as that of its directional derivative(s). As the libraries exclusively implementing the source code transformation approach collectively describe highly specialized tools that are to be used as function preprocessors, they have no direct support within deal.II itself. The latter, however, represent specialized number types that can be supported through the use of template metaprogramming in the appropriate context. Given the examples presented above, this means that the FEValues class (and friends), as well as the Tensor and SymmetricTensor classes should support calculations performed with these specialized numbers. (In theory an entire program could be made differentiable. This could be useful in, for example, the sensitivity analysis of solutions with respect to input parameters. However, to date this has not been tested.)

Implementations of specialized frameworks based on operator overloading typically fall into one of three categories. In each, some customized data classes representing the floating point value of an evaluated function and its derivative(s) by

  1. exploiting dual/complex-step/hyper-dual formulations (occasionally called tapeless methods),
  2. those utilizing taping strategies, and
  3. those using compile-time optimization through expression templates.

To provide some tentative insight into how these various implementations might look like in practice, we offer the following generic summary of these approaches:

  1. The first two tapeless approaches listed above (dual numbers and complex-step method) use some variation of a truncated Taylor series, along with a particular choice of definition for the perturbation parameter, to compute function derivatives using a finite-difference based approach. The "dual" number constitutes the accumulated directional derivatives computed simultaneously as the function values are evaluated; in the complex-step approach, the imaginary value effectively serves this purpose. The choice of the perturbation parameter determines the numerical qualities of the scheme, such as the influence of the truncation of the Taylor scheme; dual numbers do not contain any higher-order terms in their first derivative, while for the complex-step method these existent higher-order terms are neglected. It can be shown that both of these methods are not subject to subtractive cancellation errors and that, within their finite-difference scheme, they are not numerically sensitive to the internal step-size chosen for the numerical perturbation. The dual number approach thus produces exact first derivatives, while the complex-step approximation does not. The standard implementation of the dual numbers, however, cannot yield exact values for second derivatives. Hyper-dual numbers take a different view of this idea, with numbers being represented in a form similar to quaternions (i.e. carrying additional non-real components) and the derivatives being computed from a high-order truncation of the Taylor series all four components. The outcome is that, with the appropriate implementation, both first and second derivatives can be computed exactly.
  2. With taped approaches, a specified subregion of code is selected as one for which all operations executed with active (marked) input variables are tracked and recorded in a data structure referred to as a tape. At the end of the taped region, the recorded function(s) may be reevaluated by "replaying" the tape with a different set of input variables instead of recomputing the function directly. Assuming that the taped region represents a smooth function, arbitrarily high-order derivatives of the function then can be computed by referring to the code path tracked and stored on the tape. (This could perhaps be achieved, for example, through evaluation of the function around the point of interest.) There exist strategies to deal with situations where the taped function is not smooth at the evaluated point, or if it is not analytic. Furthermore, one might need to consider the case of branched functions, where the tape is no longer sequential, but rather forks off on a different evaluation path to that due to the original recorded inputs.
  3. Methods based on expression templates leverage the computational graph (in this case, a directed acyclic graph (DAG)), constructed from the abstract syntax tree (AST), that resolves the function output from its input values. The outermost leaves on the tree represent the independent variables or constants, and are transformed by unary operators and connected by binary operators (in the most simple case). Therefore, the operations performed on the function inputs is known at compile time, and with that the associated derivative operation can also be defined at the same time using the well-known rules of computing the derivative of an operation (such as the associativity of derivatives under addition and subtraction, the product rule, and the chain rule). The compiled output type returned by this operator need not be generic, but can rather be specialized based on the specific inputs (possibly carrying a differential history) given to that specific operator on the vertex of the DAG. In this way, a compile-time optimized set of instructions can be generated for the very specialized individual operations used to evaluate each intermediate result of the dependent function.

Each of these methods, of course, has its advantages and disadvantages, and one may be more appropriate than another for a given problem that is to be solved. As the aforementioned implementational details (and others not discussed) may be hidden from the user, it may still be important to understand the implications, run-time cost, and potential limitations, of using any one of these "black-box" auto-differentiable numbers.

In addition to the supplied linked articles, resources used to furnish the details supplied here include:

@InProceedings{Fike2011a,
author = {Fike, Jeffrey A and Alonso, Juan J},
title = {The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations},
booktitle = {49th {AIAA} Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition},
year = {2011},
volume = {886},
pages = {124},
month = {jan},
publisher = {American Institute of Aeronautics and Astronautics},
doi = {10.2514/6.2011-886},
}
@Manual{Walther2009a,
title = {Getting Started with ADOL-C},
author = {Walther, Andrea and Griewank, Andreas},
year = {2009},
booktitle = {Combinatorial scientific computing},
doi = {10.1.1.210.4834},
pages = {181--202}
}

of the chain-rule

In the most practical sense, any of the above categories exploit the chain-rule to compute the total derivative of a composite function. To perform this action, they typically use one of two mechanisms to compute derivatives, specifically

As a point of interest, the optimal Jacobian accumulation, which performs a minimal set of computations, lies somewhere between these two limiting cases. Its computation for a general composite function remains an open problem in graph theory.

With the aid of the diagram below (it and some of the listed details courtesy of this Wikipedia article), let us think about the representation of the calculation of the function \(f (\mathbf{x}) = \sin (x_{1}) + x_{1} x_{2}\) and its derivatives:

Forward mode automatic differentiation
Forward mode automatic differentiation
Reverse mode automatic differentiation
Reverse mode automatic differentiation

Specifically, we will briefly describe what forward and reverse auto-differentiation are. Note that in the diagram, along the edges of the graph in text are the directional derivative of function \(w\) with respect to the \(i\)-th variable, represented by the notation \(\dot{w} = \dfrac{d w}{d x_{i}}\). The specific computations used to render the function value and its directional derivatives for this example are tabulated in the source article. For a second illustrative example, we refer the interested reader to this article.

Consider first that any composite function \(f(x)\), here represented as having two independent variables, can be dissected into a composition of its elementary functions

\[ f (\mathbf{x}) = f_{0} \circ f_{1} \circ f_{2} \circ \ldots \circ f_{n} (\mathbf{x}) \quad . \]

As was previously mentioned, if each of the primitive operations \(f_{n}\) is smooth and differentiable, then the chain-rule can be universally employed to compute the total derivative of \(f\), namely \(\dfrac{d f(x)}{d \mathbf{x}}\). What distinguishes the "forward" from the "reverse" mode is how the chain-rule is evaluated, but ultimately both compute the total derivative

\[ \dfrac{d f (\mathbf{x})}{d \mathbf{x}} = \dfrac{d f_{0}}{d f_{1}} \dfrac{d f_{1}}{d f_{2}} \dfrac{d f_{2}}{d f_{3}} \ldots \dfrac{d f_{n} (\mathbf{x})}{d \mathbf{x}} \quad . \]

In forward-mode, the chain-rule is computed naturally from the "inside out". The independent variables are therefore fixed, and each sub-function \(f'_{i} \vert_{f'_{i+1}}\) is computed recursively and its result returned as inputs to the parent function. Encapsulating and fixing the order of operations using parentheses, this means that we compute

\[ \dfrac{d f (\mathbf{x})}{d \mathbf{x}} = \dfrac{d f_{0}}{d f_{1}} \left( \dfrac{d f_{1}}{d f_{2}} \left(\dfrac{d f_{2}}{d f_{3}} \left(\ldots \left( \dfrac{d f_{n} (\mathbf{x})}{d \mathbf{x}} \right)\right)\right)\right) \quad . \]

The computational complexity of a forward-sweep is proportional to that of the input function. However, for each directional derivative that is to be computed one sweep of the computational graph is required.

In reverse-mode, the chain-rule is computed somewhat unnaturally from the "outside in". The values of the dependent variables first get computed and fixed, and then the preceding differential operations are evaluated and multiplied in succession with the previous results from left to right. Again, if we encapsulate and fix the order of operations using parentheses, this implies that the reverse calculation is performed by

\[ \dfrac{d f (\mathbf{x})}{d \mathbf{x}} = \left( \left( \left( \left( \left( \dfrac{d f_{0}}{d f_{1}} \right) \dfrac{d f_{1}}{d f_{2}} \right) \dfrac{d f_{2}}{d f_{3}} \right) \ldots \right) \dfrac{d f_{n} (\mathbf{x})}{d \mathbf{x}} \right) \quad . \]

The intermediate values \(\dfrac{d f_{i-1}}{d f_{i}}\) are known as adjoints, which must be computed and stored as the computational graph is traversed. However, for each dependent scalar function one sweep of the computational graph renders all directional derivatives at once.

Overall, the efficiency of each mode is determined by the number of independent (input) variables and dependent (output) variables. If the outputs greatly exceed the inputs in number, then forward-mode can be shown to be more efficient than reverse-mode. The converse is true when the number of input variables greatly exceeds that of the output variables. This point may be used to help inform which number type is most suitable for which set of operations are to be performed using automatic differentiation. For example, in many applications for which second derivatives are to be computed it is appropriate to combine both reverse- and forward-modes. The former would then typically be used to calculate the first derivatives, and the latter the second derivatives.

Supported automatic differentiation libraries

We currently have validated implementations for the following number types and combinations:

Note that in the above, "dynamic memory allocation" refers to the fact that the number of independent variables need not be specified at compile time.

The ADOL-C user manual

@Manual{Walther2009a,
title = {Getting Started with ADOL-C},
author = {Walther, Andrea and Griewank, Andreas},
year = {2009},
booktitle = {Combinatorial scientific computing},
doi = {10.1.1.210.4834},
pages = {181--202},
url = {https://projects.coin-or.org/ADOL-C/browser/trunk/ADOL-C/doc/adolc-manual.pdf}
}

provides the principle insights into their taped and tapeless implementations, and how ADOL-C can be incorporated into a user code. Some further useful resources for understanding the implementation of ADOL-C, and possibilities for how it may be used within a numerical code, include:

@Article{Griewank1996a,
author = {Griewank, Andreas and Juedes, David and Utke, Jean},
title = {Algorithm 755: {ADOL-C}: a package for the automatic differentiation of algorithms written in {C/C++}},
journal = {ACM Transactions on Mathematical Software (TOMS)},
year = {1996},
volume = {22},
number = {2},
pages = {131--167},
doi = {10.1145/229473.229474},
publisher = {ACM}
}
@InCollection{Bischof2008a,
author = {Bischof, Christian and Guertler, Niels and Kowarz, Andreas and Walther, Andrea},
title = {Parallel reverse mode automatic differentiation for OpenMP programs with ADOL-C},
booktitle = {Advances in Automatic Differentiation},
publisher = {Springer},
year = {2008},
pages = {163--173}
}
@InBook{Kulshreshtha2012a,
chapter = {Computing Derivatives in a Meshless Simulation Using Permutations in {ADOL}-C},
pages = {321--331},
title = {Recent Advances in Algorithmic Differentiation},
publisher = {Springer Berlin Heidelberg},
year = {2012},
author = {Kshitij Kulshreshtha and Jan Marburger},
editor = {Forth S. and Hovland P. and Phipps E. and Utke J. and Walther A.},
series = {Lecture Notes in Computational Science and Engineering},
doi = {10.1007/978-3-642-30023-3_29},
}
@InProceedings{Kulshreshtha2013a,
author = {Kulshreshtha, Kshitij and Koniaeva, Alina},
title = {Vectorizing the forward mode of ADOL-C on a GPU using CUDA},
booktitle = {13th European AD Workshop},
year = {2013},
month = jun
}

Similarly, a selection of useful resources for understanding the implementation of Sacado number types (in particular, how expression templating is employed and exploited) include:

@InCollection{Bartlett2006a,
author = {Bartlett, R. A. and Gay, D. M. and Phipps, E. T.},
title = {Automatic Differentiation of C++ Codes for Large-Scale Scientific Computing},
booktitle = {International Conference on Computational Science {\textendash} {ICCS} 2006},
publisher = {Springer Berlin Heidelberg},
year = {2006},
editor = {Alexandrov, V.N. and van Albada, G.D. and Sloot, P.M.A. amd Dongarra, J.},
pages = {525--532},
doi = {10.1007/11758549_73},
organization = {Springer}
}
@InBook{Gay2012a,
chapter = {Using expression graphs in optimization algorithms},
pages = {247--262},
title = {Mixed Integer Nonlinear Programming},
publisher = {Springer New York},
year = {2012},
author = {Gay, D. M.},
editor = {Lee, J. and Leyffer, S.},
isbn = {978-1-4614-1927-3},
doi = {10.1007/978-1-4614-1927-3_8}
}
@InBook{Phipps2012a,
chapter = {Efficient Expression Templates for Operator Overloading-based Automatic Differentiation},
pages = {309--319},
title = {Recent Advances in Algorithmic Differentiation},
publisher = {Springer},
year = {2012},
author = {Eric Phipps and Roger Pawlowski},
editor = {Forth S. and Hovland P. and Phipps E. and Utke J. and Walther A.},
series = {Lecture Notes in Computational Science and Engineering},
volume = {73},
date = {2012-05-15},
doi = {10.1007/978-3-642-30023-3_28},
eprint = {1205.3506v1},
eprintclass = {cs.MS},
eprinttype = {arXiv}
}

The implementation of both forward- and reverse-mode Sacado numbers is quite intricate. As of Trilinos 12.12, the implementation of math operations involves a lot of preprocessor directives and macro programming. Accordingly, the code may be hard to follow and there exists no meaningful companion documentation for these classes. So, a useful resource for understanding the principle implementation of these numbers can be found at this link for the Sacado::Fad::SimpleFad class that outlines a reference (although reportedly inefficient) implementation of a forward-mode auto-differentiable number that does not use expression templates. (Although not explicitly stated, it would appear that the Sacado::Fad::SimpleFad class is implemented in the spirit of dual numbers.)

How automatic differentiation is integrated into deal.II

Since the interface to each automatic differentiation library is so vastly different, a uniform internal interface to each number will be established in the near future. The goal will be to allow some driver classes (that provide the core functionality, and will later be introduced in the next section) to have a consistent mechanism to interact with different auto-differentiation libraries. Specifically, they need to be able to correctly initialize and finalize data that is to be interpreted as the dependent and independent variables of a formula.

A summary of the files that implement the interface to the supported auto-differentiable numbers is as follows:

By using type codes for each supported number type, we artificially limit the type of auto-differentiable numbers that can be used within the library. This design choice is due to the fact that its not trivial to ensure that each number type is correctly initialized and that all combinations of nested (templated) types remain valid for all operations performed by the library. Furthermore, there are some lengthy functions within the library that are instantiated for the supported number types and have internal checks that are only satisfied when a auto-differentiable number, of which the library has knowledge, is used. This again ensures that the integrity of all computations is maintained. Finally, using a simple enumeration as a class template parameter ultimately makes it really easy to switch between the type used in production code with little to no further amendments required to user code.

User interface to the automatic differentiation libraries

The deal.II library offers a unified interface to the automatic differentiation libraries that we support. To date, the helper classes have been developed for the following contexts:

Naturally, it is also possible for users to manage the initialization and derivative computations themselves.

The most up-to-date examples of how this is done using ADOL-C can be found in

while for Sacado, illustrative examples can be found in

Symbolic expressions and differentiation

Symbolic differentiation is, in terms of its design and usage, quite different to automatic differentiation. Underlying any symbolic library is a computer algebra system (CAS) that implements a language and collection of algorithms to manipulate symbolic (or "string-like") expressions. This is most similar, from a philosophical point of view, to how algebraic operations would be performed by hand.

To help better distinguish between symbolic differentiation and numerical methods like automatic differentiation, let's consider a very simple example. Suppose that the function \(f(x,y) = [2x+1]^{y}\), where \(x\) and \(y\) are variables that are independent of one another. By applying the chain-rule, the derivatives of this function are simply \(\dfrac{d f(x,y)}{d x} = 2y[2x+1]^{y-1}\) and \(\dfrac{d f(x,y)}{d y} = [2x+1]^{y} \ln(2x+1)\). These are exactly the results that you get from a CAS after defining the symbolic variables x and y, defining the symbolic expression f = pow(2x+1, y) and computing the derivatives diff(f, x) and diff(f, y). At this point there is no assumption of what x and y represent; they may later be interpreted as plain (scalar) numbers, complex numbers, or something else for which the power and natural logarithm functions are well defined. Obviously this means that there is also no assumption about which point to evaluate either the expression or its derivatives. One could readily take the expression for \(\dfrac{d f(x, y)}{d x}\) and evaluate it at \(x=1, y=2.5\) and then later, with no recomputation of the derivative expression itself, evaluate it at \(x=3.25, y=-6\). In fact, the interpretation of any symbolic variable or expression, and the inter-dependencies between variables, may be defined or redefined at any point during their manipulation; this leads to a degree of flexibility in computations that cannot be matched by auto-differentiation. For example, one could perform the permanent substitution \(g(x) = \dfrac{d f(x, y)}{d x} \vert_{y=1}\) and then recompute \(g(x)\) for several different values of \(x\). One could also post-factum express an interdependency between x and y, such as \(y \rightarrow y(x) := 2x\). For such a case, this means that the initially computed derivatives \(\dfrac{d f(x, y)}{d x} \rightarrow \dfrac{\partial f(x, y(x))}{\partial x} = 2y(x) [2x+1]^{y(x)-1} = 4x[2x+1]^{2x-1}\) and \(\dfrac{d f(x, y)}{d y} \rightarrow \dfrac{\partial f(x, y(x))}{\partial y} = [2x+1]^{y(x)} \ln(2x+1) = [2x+1]^{2x} \ln(2x+1)\) truly represent partial derivatives rather than total derivatives. Of course, if such an inter-dependency was explicitly defined before the derivatives \(\dfrac{d f(x, y(x))}{d x}\) and \(\dfrac{d f(x, y(x))}{d y}\) are computed, then this could correspond to the total derivative (which is the only result that auto-differentiation is able to achieve for this example).

Due to the sophisticated CAS that forms the foundation of symbolic operations, the types of manipulations are not necessarily restricted to differentiation alone, but rather may span a spectrum of manipulations relevant to discrete differential calculus, topics in pure mathematics, and more. The documentation for the SymPy library gives plenty of examples that highlight what a fully-fledged CAS is capable of. Through the Differentiation::SD::Expression class, and the associated functions in the Differentiation::SD namespace, we provide a wrapper to the high-performance SymEngine symbolic manipulation library that has enriched operator overloading and a consistent interface that makes it easy and "natural" to use. In fact, this class can be used as a "drop-in" replacement for arithmetic types in many situations, transforming the operations from being numeric to symbolic in nature; this is made especially easy when classes are templated on the underlying number type. Being focused on numerical simulation of PDEs, the functionality of the CAS that is exposed within deal.II focuses on symbolic expression creation, manipulation, and differentiation.

The convenience wrappers to SymEngine functionality are primarily focused on manipulations that solely involve dictionary-based (i.e., something reminiscent of "string-based") operations. Although SymEngine performs these operations in an efficient manner, they are still known to be computationally expensive, especially when the operations are performed on large expressions. It should therefore be expected that the performance of the parts of code that perform differentiation, symbolic substitution, etc., may be a limiting factor when using this in production code. deal.II therefore provides an interface to accelerate the evaluation of lengthy symbolic expression through the BatchOptimizer class (itself often leveraging functionality provided by SymEngine). In particular, the BatchOptimizer simultaneously optimizes a collection of symbolic expressions using methods such as common subexpression elimination (CSE), as well as by generating high performance code-paths to evaluate these expressions through the use of a custom-generated std::function or by compiling the expression using the LLVM JIT compiler. The usage of the Differentiation::SD::BatchOptimizer class is exemplified in step-71.

As a final note, it is important to recognize the remaining major deficiencies in deal.II's current implementation of the interface to the supported symbolic library. The level of functionality currently implemented effectively limits the use of symbolic algebra to the traditional use case (i.e. scalar and tensor algebra, as might be useful to define constitutive relations or complex functions for application as boundary conditions or source terms). In fact, step-71 demonstrates how it can be used to implement challenging constitutive models. In the future we will also implement classes to assist in performing assembly operations in the same spirit as that which has been done in the Differentiation::AD namespace.

A summary of the files that implement the interface to the supported symbolic differentiable numbers is as follows: