 |
Reference documentation for deal.II version 9.3.3
|
|
Reference documentation for deal.II version 9.3.3
|
A module 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 |
A module 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 (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
To provide some tentative insight into how these various implementations might look like in practice, we offer the following generic summary of these approaches:
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:
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 represention of the calculation of the function \(f (\mathbf{x}) = \sin (x_{1}) + x_{1} x_{2}\) and its derivatives:
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.
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.
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:
Similarly, a selection of useful resources for understanding the implementation of Sacado number types (in particular, how expression templating is employed and exploited) include:
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.)
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.
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 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:
1.9.4