Programme of the 14th Euro AD Workshop

The operator-overloading approach to performing reverse-mode automatic differentiation (AD) is more convenient and flexible for the user than source-code transformation, but current implementations are significantly slower, and typically 10 to 35 times slower than the original algorithm. In this talk, two methods to speed up reverse-mode operator-overloading AD will be discussed:

1) The expression template programming technique in C++ provides a compile-time representation of each mathematical expression as a computational graph that can be efficiently traversed in either direction, significantly speeding up the forward pass. This is the first time this technique has been used for full reverse-mode AD.

2) Use of a minimalist tape structure containing only the differential statements leads to a much faster reverse pass than libraries that store the entire algorithm, as well as being less demanding in memory usage. Moreover, computation of Jacobian matrices from this tape structure (either in forward or reverse mode) is easily parallelized.

These ideas have been implemented in the freely available C++ "Adept" (Automatic Differentiation using Expression Templates) software library. Benchmarking with four different numerical algorithms shows reverse-mode AD with Adept to be 2.6 to 9 times faster and 1.3 to 7.7 times more efficient in memory usage than the established operator-overloading libraries ADOL-C, CppAD and Sacado. For typical algorithms containing mathematical functions (exp, log etc.), Adept is found to be only 2.7 to 3.9 times the cost of the original algorithm, which is very competitive compared to source-code transformation. However, for simple loops containing no mathematical functions, it is found that the relative cost of all operator-overloading libraries (including Adept) is much greater, due to the AD impeding the aggressive optimization that compilers perform in the undifferentiated versions of such loops.

In modern simulation software the deployment of template functions and classes is a fairly common concept as it reduces the amount of repetitive code and increases the readability. For instance, it is not necessary to rewrite complete classes when the values of the simulation are changed from complex to real. This simplifies the top level programming when implementing a new simulation case. However, when incorporating a new data type originally unintended for use with the software, such as an operator-overloading based AD class, the combination of templated function signatures and operator overloading proves to be challenging. In this talk, the challenges of incorporating ADOL-C and its data type 'adouble' into a template based simulation software are presented.

An examination of the mathematical principles underlying high performance implementations of AD reveals general-purpose mechanisms for program transformation, coupled with specifics related to the particulars of taking derivatives. These program transformation mechanisms would appear to have application beyond AD, to a variety of programming contexts. Here we present preliminary work on a general calculus of program transformation based on the notion of ``lifting'' a computation from its primal category to another related category. After a generic description of the novel mechanism, we show how forward AD in particular can be expressed by specification of an appropriate algebra combined with a generic forward lifting operator, and give examples of other useful tasks which can be accomplished by using other algebras.

We present some recent developments in Chebfun for differenctial equations which give rise to systems with block structure. The blocks consists of operators, functions, functionals and scalars; AD is used to compute the blocks from the original nonlinear operator.

This talk describes an algorithm for robust optimization consisting of two main independently developed components. The algorithm combines uncertainty quantification and global optimization (minimization) to a robust optimization. First uncertainty quantification is used to produce a robust version of the functional $F: R^{n} \rightarrow R$. Subsequently a McCormick based branch and bound optimization approach is applied to this robust function.

The uncertainty quantification is based on moments method. Moments method assumes the inputs of the underlying vector function to be probabilistic because of possible measurement errors. Therefore the mean and variance of the outputs are approximated by a Taylor series. Algorithmic Differentiation (AD) is employed to compute the required derivatives either in tangent-linear or adjoint mode. The second component is a deterministic branch and bound algorithm for global optimization based on McCormick relaxations. McCormick relaxations are convex underestimators of factorable functions which are possibly nonsmooth. Their subgradients are used instead of interval analysis in the branch and bound method. The fact that the computation of subgradients of McCormick relaxations is equivalent to AD was shown in previous work.

The combination of these two procedures ensures a more stable minimum accounting the possible uncertainties in the inputs.

Over the last few decades there has been a lot of research into nonsmooth analysis. Mainly theoretical, but also with some hand coded applications. So far no AD like automization has taken place to my knowledge. I strongly feel that the AD community can no longer ignore these developments and see great opportunities ahead.

Numerical analysis has mostly been concerned with smooth problems and in optimization also with convex problems. A joint property of smooth and convex problems is 'semi-smoothness'. This concept that has risen to fame as the key assumption in the local convergence analysis for generalized Newton methods by Qi and Sun. All Lipschitz continuous compositions of smooth elementals and max/min/abs have this property. For these, as shown by Barton & Khan and myself generalized Jacobians of such functions can be evaluated by a minor extension of standard AD technology.

The key concept is piecewise linearization as a first order approximation of piecewise smooth functions. Local piecewise linear models can be used to attack the standard numerical tasks of equation solving, (un)constrained optimization and more. A package of common numerical utilities called PLAN-C for Piecewie Linear functions in Abs-Normal form in C is currently under development. Preliminary solution concepts and results are detailed in the subsequent talks.

It follows from the well known min/max representation given by Scholtes

In the smooth case the steepest descent method proved to be a reliable choice for the minimization of an unconstrained objective function. However, for non-smooth functions it was shown that this step direction in combination with an exact line-seacrh might exhibit a zigzagging convergence to non-stationary points. Therefore, we propose a bundle method based on the piecewise linearization for composite Lipschitzian optimization problems. Here the step directions are computed using a piecewise linearization of the original function plus a proximal term. The necessary information for the approximation of the objective are obtained by techniques from Automatic Differentiation.

In praxis one often encounters ODEs that have a nonsmooth right hand side. Classical methods tend to lose some accuracy especially on problems with many indifferentiabi- lities in the considered interval or may even fail to converge. Thus we want to generalize two of these methods, the midpoint and the trapezoidal rule, using a piecewise linear approximation of the right hand side of the ODE. These generalized methods can be shown to converge with second order and proved to be more stable in our numerical experiments.

An implementation of the discrete adjoint of an industrial coupled flow CFD solver was begun in 2011. The project was constrained by the infeasibility of using standard operator overloading or source transformation approaches to generate the derivative algorithms. Since a manual implementation appeared to be the only feasible approach, ways to minimise this implementation effort, maintain acceptable execution times and memory overhead were sought. One aspect of the approach used is to perform loop-local operator overloading of templated tensor types; i.e. to devise an interpreter that can evaluate the adjoint for zero, first and second order tensors as native types. The details of this, along with other pertinent aspects of the implementation, are discussed.

Secondly, given the significant developments in compile-time evaluation, is it conceivable to have the compiler generate adjoint code from a domain-specific embedded language? Some very brief considerations are presented to indicate its possibility in the D programming language.

In the context of working on adjoints for an ice-sheet model some new functionality had to be added. Most of the work went into implementing the support for the AdjoinableMPI library and making its use in the model code mostly transparent. Additional work addressed the handling of (sparse) direct solvers and the management of the Adol-C internal pseudo addresses. Without going into the implementation details I will briefly survey the last two areas while spending most of the time on the AdjoinableMPI work and lessons learned from the application to the ice sheet model.

This talk will comprise of two parts. The first part will describe the new implementation of linearisation in Abs-Normal Form for piecewise smooth functions, as proposed by Andreas Griewank. The new aspect is that the same ADOL-C trace can be used to compute Abs-Normal Form extended Jacobian as well as any traditional derivative information, using the appropriate interface functions. The second part will give a how-to for efficiently differentiating through solutions of linear systems, perhaps using LAPACK or any other linear solver, using the newly extended and debugged external function interface of ADOL-C. Such systems arise in many applications where the function evaluation itself comprises of or is computed from solutions to some intermediate linear systems. This is a frequently asked question by many users of ADOL-C.

OpenFOAM (Open Field Operation and Manipulation) is a C++ toolbox for the de- velopment of customized numerical solvers, and pre-/post-processing utilities for the solution of continuum mechanics problems, including computational fluid dynamics (CFD). The code is released as free and open source software under the GNU General Public License. We introduced a discrete adjoint version of OpenFOAM using the operator overload- ing tool dco/c++ in order to generate derivatives and to apply optimization (mainly focussing on topology optimization)[3]. A discrete adjoint implementation in general yields a significant overhead to the passive evaluation in computation time and more importantly required memory. (Intermediate values from the whole computational history have to be stored in order to evaluate the derivatives). This talk will mainly focus on how we managed to significantly reduce the memory requirements of said discrete adjoint OpenFOAM version by applying analytical knowledge about the linear solvers used to solve the underlying partial differential equations[1], thus eliminating the need to store intermediate values generated during the solution process of the linearized equations. This linear solver treatment yields a significant improvement in both runtime and memory usage and also eliminates the dependency of the memory usage on the number of linear solver iterations. The used amount of memory still remains high but is much more manageable with checkpointing techniques (e.g. [2]). Further I want to discuss some problems we encountered with our discrete implementation, in particular the known phenomenon that the adjoints of a discrete implementation of the CG-Algorithm differ from the analytically obtained ones.

References

[1] M. Giles. Collected matrix derivative results for forward and reverse mode algorithmic differentiation. Advances in Automatic Differentiation, pages 35–44, 2008.

[2] A. Griewank and A. Walther. Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation. ACM Transactions on Mathematical Software, 26(1):19–45, Mar. 2000.

[3] M. Towara and U. Naumann. A Discrete Adjoint Model for OpenFOAM. Procedia Computer Science, 18(0):429 – 438, 2013.

Adjoint methods are probably the most efficient way to obtain the gradient of a numerical simulation, more precisely the gradient of a cost function with respect to design parameters. Theses simulations often contain Fixed-Point algorithms, whose derivatives have been studied theoretically by several authors. We will focus on the methods described in articles by Griewank-Walther and by Christianson. We will show how these methods can be provided in a tool such as Tapenade. We plan to ask the audience about their effective use of Fixed-Point simulations and the benefit that could be expected from a specific adjoint.

Obtaining sensitivities with the discrete adjoint approach is a very popular approach as it is automatable through AD and as the exact transposition of the system Jacobian (to machine precision) guarantees stability. Typically in CFD design this approach is applied to steady-state flow solvers using Christianson's fixed-point accumulation: all we need to do is to snapshot the converged flowfield and compute its Jacobian.

In practice though we find that things don't work as advertised. For complex cases on fine meshes the primal does not converge to machine zero, but stalls in limit cycles: unstable modes appear, grow a bit until the non-linear discretisation of the primal adapts and ``swats'' it, only to have other unstable modes appear. This is normally not an issue for the primal, resulting fluctuations in the cost function are for the most negligible. However, blindly applying the steady-state approach means taking a snapshot at an arbitrary iteration within the limit cycle, which invariably will have non-contractive eigenvalues, the iterative adjoint solver will blow up.

The presentation does not directly deal with the application of AD, but highlights how a important a well-conditioned and stable Jacobian is for converging a discrete adjoint. We will briefly point out remedies from the literature such as RPM to tackle these unstable modes, to then show how an improved numerical stability of the primal can improve the eigenvalue spectrum to become contractive and consequently achieve adjoint convergence, making it work as advertised. Examples from turbo-machinery will be presented.

We introduce an efficient extension of a recently introduced auto-validating rejection sampler that is capable of producing indepen- dent and identically distributed (IID) samples from a large class of target densities with locally Lipschitz arithmetical expressions. Our extension is restricted to target densities that are differentiable. We use the centered form, as opposed to the natural interval extension, to get tighter range enclosures of the differentiable multivariate target density using interval extended gradient differentiation arithmetic. By using the centered form we are able to sample one hundred times faster from the posterior density over the space of phylogenetic trees with four leaves (quartets).

The crucial idea here is to obtain tighter range enclosures of the target function. A regular paving is a finite succession of bisections that partition a root box x in R^d into sub-boxes using a binary tree-based data structure. The sequence of splits that generate such a partition is given by the sub-boxes associated with the nodes of the tree. The leaf boxes, i.e., the sub-boxes associated with the leaf nodes, form a partition of x. We employ randomized algorithms to tightly enclose the range of a function g using its interval extension. Our idea is to (i) refine the regular paving partition of the domain x by splitting the leaf boxes, (ii) obtain range enclosures of g over them, (iii) propagate the range enclosures of the leaf boxes up the internal nodes of the tree and finally (iv) prune back the leaves to get a coarser partition with fewer leaf boxes but with tighter range enclosures. This approach allows one to obtain tighter range enclosures for interval inclusion functions and thereby increase the acceptance probability of rejection samplers.

Data flow reversal plays a grave role in the adjoint mode of algorithmic differentiation. To tackle the arising memory issue, checkpointing strategies are employed, which give the possibility to trade memory consumption in additional runtime. Computer programs can be interpreted as directed acyclic graphs (DAG) which transfers the problem to graphs. Finding an optimal checkpointing on DAGs is known as the DAG reversal problem and proven to be NP-complete. The problem size for the DAG reversal problem gets infeasible for large computer programs very quickly. For the call tree reversal (CTR), we therefore restrict the general checkpointing approach on the DAG to subroutine calls in a given implementation. This yields a much smaller problem size but is still an effective strategy for checkpointing. For this NP-complete combinatorial optimization problem, we develop an integer programming formulation and we succesfully apply standard algorithms for an efficient solution thereof. We also work on the toolchain from the untouched C++ computer program to the solution of CTR and the implementation of the resulting checkpointing scheme.

Large-scale combinatorial scientific computing problems arising in sparse or otherwise structured matrix computation are often expressed using an appropriate graph model, and sometime the same problem can be given in more than one graph model with similar asymptotic computational complexity. In a recent work we have proposed a new graph model, the pattern graph, to formulate combinatorial problems arising in sparse Jacobian matrix determination. In this presentation, we extend the pattern graph model to incorporate matrix symmetry.

Given our community's habit of alternating between US and European venues it is Europe's time to host the next International Conference on Algorithmic Differentiation in Summer 2016. We will present a short list of possible venues and solicit any final hosting offers with a view to a community vote in early 2014 to finalise the venue.