We consider implementations of multivariate vector functions $$ f : D^{n} \times D^{n'} \rightarrow D^m \times D^{m'} : (\Y, \Y') = f(\X,\X') $$ as computer programs over some base data type $D$ (for example, single or higher precision floating-point data, intervals, convex/concave relaxations, vectors/ensembles).

We assume that the arithmetic inside $f$ is completely
defined (through overloading of the *elemental functions*; see below)
.

The $n$ active (also: independent)
and $n'$ passive inputs are mapped onto $m$ active
(also: dependent) and $m'$ passive outputs. The given
implementation is assumed to be $k$ times continuously
differentiable at all points of interest implying the
existence and finiteness of the Jacobian
$$
\nabla f = \nabla_\X f(\X,\X') \equiv
\frac{\partial \Y}{\partial \X}(\X,\X') \in D^{m \times n},
$$
the Hessian
$$
\nabla^2 f = \nabla^2 f(\X,\X') \equiv
\frac{\partial^2 \Y}{\partial \X^2}(\X,\X') \in D^{m \times n \times n},
$$
if $k\geq2,$ and of potentially higher derivative tensors
$$
\nabla^k f = \nabla^k f(\X,\X') \equiv
\frac{\partial^k \Y}{\partial \X^k}(\X,\X') \in D^{m \times n \times \underset{k \text{times}}{\ldots} \times n}.
$$
We denote $$\nabla^k f=\left(\left [\nabla^k f \right ]_i^{j_1,\ldots,j_k} \right)_{i=0,\ldots,m-1}^{j_1,\ldots,j_k=0,\ldots,n-1},$$
and we use $*$ to denote the entire range of an index.
For example,
$\left [\nabla f \right ]_i^*$ denotes the $i$th row
and
$\left [\nabla f \right ]_*^j$ the $j$th column of the
Jacobian, respectively.
Algorithmic differentiation is implemented by the overloading
of elemental functions including the built-in functions
and operators of C++ as well as user-defined higher-level
elemental functions.
Chunks are written to and read from disk.

The global tapes are now thread-local by default allowing an easier switch to multi-threading.

Data types for values, partial derivatives and adjoints can be set individually when defining the differentiation mode.

Multiple “blob” or “chunk” tapes can be allocated, for example, to implement thread-safe adjoints.

- checkpointing
- symbolic adjoints of numerical methods
- combinations of tangent and adjoint modes
- preaccumulation of local derivative tensors
- creation of collections of domain-specific higher-level elemental functions

An easy-to-use interface for robust preaccumulation of local Jacobians is provided. Aggressive reduction of the tape size is facilitated.

Several (concurrent) interpretations of the same tape are enabled through separate (thread-local) adjoint vectors.

The vector of adjoints is compressed by analysing the maximum number of required distinct adjoint memory locations. During interpretation, adjoint memory, which is no longer required, is overwritten and thus reused by indexing through modulo operations. This feature is especially useful for iterative algorithms (e.g. time iterations). The required memory for the vector of adjoints usually stays constant, independent of the number of iterations. Combined with the use of disk tape, almost arbitrarily sized tapes can be generated, which might be especially of interest for prototyping or validation purposes.

The adjoint interpretation of the tape can omit propagation along edges when the corresponding adjoint to be propagated is zero. This might be of use when NaNs or Infs occur as local partial derivatives (e.g. when computing square root of zero), but this local result is only used for subsequent computations which are not relevant for the overall output. This feature might have a performance impact on tape interpretation and should therefore be considered for debugging configuration only.

Instantiation of the generic first-order vector tangent data type with the generic first-order scalar tangent data type over a non-derivative base data type yields a second-order tangent data type. Similarly, instantiation of the generic first-order scalar tangent data type with the generic first-order vector tangent data type over a non-derivative base data type yields a second-order tangent data type.

Instantiation of the generic first-order scalar adjoint data type with the generic first-order scalar tangent data type over a non-derivative base data type yields a second-order scalar adjoint data type. It enables the computation of scalar projections of the Hessian $\nabla^2 f \in D^{m \times n \times n}$ in its image and one of its domain dimensions as $$ \X_{(1)}^{(2)}=\langle\Y_{(1)},\nabla^2 f,\X^{(2)}\rangle \equiv \left(\Y_{(1)}^T \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot \X^{(2)} \right)_{j=0,\ldots,n-1} \in D^n $$ for $\Y_{(1)} \in D^m$ and $\X^{(2)} \in D^n$ over arbitrary base data types $D$. The computational cost of accumulating the whole Hessian over $D$ becomes $O(m\cdot n) \cdot \text{Cost}(f).$ with $\Y_{(1)}$ and $\X^{(2)}$ ranging over the Cartesian basis vectors in $D^m$ and $D^n$, respectively.

Instantiation of the generic first-order vector adjoint data type with the generic first-order vector tangent data type over a non-derivative base data type yields a second-order vector adjoint data type. It enables the computation of vector projections of the Hessian $\nabla^2 f \in D^{m \times n \times n}$ in its image and one of its domain dimensions as $$ X_{(1)}^{(2)}=\langle Y_{(1)},\nabla^2 f,X^{(2)}\rangle \equiv \left(Y_{(1)}^T \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot X^{(2)} \right)_{j=0,\ldots,n-1} \in D^{l_1 \times n \times l_2} $$ for $Y_{(1)} \in D^{m \times l_1}$ and $X^{(2)} \in D^{n \times l_2}$ over arbitrary base data types $D$. The computational cost of accumulating the whole Hessian over $D$ remains equal to $O(m\cdot n) \cdot \text{Cost}(f).$ with $Y_{(1)}$ and $X^{(2)}$ set equal to the identities in $D^m$ and $D^n$, respectively.

Instantiation of the generic first-order vector adjoint data type with the generic first-order scalar tangent data type over a non-derivative base data type yields a second-order adjoint data type. Similarly, instantiation of the generic first-order scalar adjoint data type with the generic first-order vector tangent data type over a non-derivative base data type yields a second-order adjoint data type.

Symmetry of the Hessian in its two domain dimensions yields the following additional second-order adjoint data types: Instantiation of a generic first-order tangent data type with a generic first-order adjoint data type over a non-derivative base data type yields a second-order adjoint data type for computing $$ \langle\Y_{(2)}^{(1)},\nabla^2 f,\X^{(1)}\rangle \equiv \left(\Y_{(2)}^{(1)^T} \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot \X^{(1)} \right)_{j=0,\ldots,n-1} \in D^n. $$ Similarly, instantiation of a generic first-order adjoint data type with a generic first-order adjoint data type over a non-derivative base data type yields a second-order adjoint data type for computing $$ \langle\X_{(1,2)},\Y_{(1)},\nabla^2 f\rangle \equiv \left(\Y_{(1)}^T \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot \X_{(1,2)} \right)_{j=0,\ldots,n-1} \in D^n. $$

For making use of vectorization, a respective vector data type can be used as base type for tangent and adjoint modes.

For enabling internal debugging checks, a compile time flag was added.