
Programme of the 14th Euro AD Workshop
Monday, December 9, 2013
 9^{00}  9^{30} Arrival
 9^{30}  11^{00} Operator overloading and program transformation
 Robin Hogan (University of Reading)
Can operatoroverloading ever have a speed approaching sourcecode transformation for reversemode automatic differentiation?
The operatoroverloading approach to performing reversemode automatic differentiation (AD) is more convenient and flexible for the user than sourcecode 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 reversemode operatoroverloading AD will be discussed:
1) The expression template programming technique in C++ provides a compiletime 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 reversemode 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 reversemode 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 operatoroverloading libraries ADOLC, 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 sourcecode transformation. However, for simple loops containing no mathematical functions, it is found that the relative cost of all operatoroverloading libraries (including Adept) is much greater, due to the AD impeding the aggressive optimization that compilers perform in the undifferentiated versions of such loops.
 Benjamin Jurgelucks (Universität Paderborn)
Challenges of incorporating ADOLC into template based simulation software
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 operatoroverloading based AD class, the combination of templated function signatures and operator overloading proves to be challenging.
In this talk, the challenges of incorporating ADOLC and its data type 'adouble' into a template based simulation software are presented.
 David Rush (National University of Ireland Maynooth)
Refactoring AD: Plumbing vs. Numerics
An examination of the mathematical principles underlying high performance
implementations of AD reveals generalpurpose 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.
 11^{00}  11^{30} Coffee
 11^{30}  12^{30} Applications of AD
 Asgeir Birkisson (University of Oxford)
Block systems in Chebfun
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.
 Markus Beckers (RWTH Aachen)
Global robust optimization using McCormick relaxations
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 tangentlinear 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.
 12^{30}  14^{00} Lunch
 14^{00}  16^{30} Nonsmooth numerics via ADbased piecewise linearization
 Andreas Griewank (HumboldtUniversität zu Berlin)
Stable Generalized differentiation and piecewise linearization via AD (1 hour)
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 'semismoothness'. 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 PLANC for Piecewie Linear functions in
AbsNormal form in C is currently under development. Preliminary solution
concepts and results are detailed in the subsequent talks.
 Tom Streubel (HumboldtUniversität zu Berlin)
Representation and analysis of continuous piecewise linear functions in absnormal form
It follows from the well known min/max representation given by Scholtes
 Torsten Bosse (HumboldtUniversität zu Berlin)
(Un)constrained optimization of PL functions
In the smooth case the steepest descent method proved to be a reliable
choice for the minimization of an unconstrained objective function.
However, for nonsmooth functions it was shown that this step direction in
combination with an exact lineseacrh might exhibit a zigzagging
convergence to nonstationary 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.
 Richard Hasenfelder (HumboldtUniversität zu Berlin)
Numerical integration of some Lipschitzian ODEs using PL models
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.
 16^{30}  19^{30} End
 19^{30} Conference dinner at Pierre Victoire, Little Clarendon St

Tuesday, December 10, 2013
 9^{00}  11^{00} Applications and tutorials
 Dominic Jones (CDadapco)
Discrete adjoint; an industrial perspective
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 looplocal 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 compiletime evaluation, is it conceivable to have the compiler generate adjoint code from a domainspecific embedded language? Some very brief considerations are presented to indicate its possibility in the D programming language.
 Jean Utke (Argonne National Laboratory)
Things I recently did to AdolC
In the context of working on adjoints for an icesheet 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 AdolC 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.
 Kshitij Kulshreshtha (Universität Paderborn)
A new development and a howto for ADOLC
This talk will comprise of two parts. The first part will describe the new
implementation of linearisation in AbsNormal Form for piecewise smooth
functions, as proposed by Andreas Griewank. The new aspect is that the
same ADOLC trace can be used to compute AbsNormal Form extended Jacobian
as well as any traditional derivative information, using the appropriate
interface functions. The second part will give a howto 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 ADOLC. 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 ADOLC.
 Markus Towara (RWTH Aachen)
An effective discrete adjoint model for OpenFOAM
OpenFOAM (Open Field Operation and Manipulation) is a C++ toolbox for the de
velopment of customized numerical solvers, and pre/postprocessing 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 CGAlgorithm 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.
 11^{00}  11^{30} Coffee
 11^{30}  12^{30} Fixedpoint programs
 Ala Taftaf (INRIA SophiaAntipolis, France)
Adjoints of fixedpoint programs
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 FixedPoint algorithms, whose derivatives have been studied theoretically by several authors. We will focus on the methods described in articles by GriewankWalther 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 FixedPoint simulations and the benefit that could be expected from a specific adjoint.
 JensDominik Müller (Queen Mary University of London)
Converging discrete adjoint solutions: influence of primal stability
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 steadystate flow
solvers using Christianson's fixedpoint 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 nonlinear 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 steadystate
approach means taking a snapshot at an arbitrary iteration within the
limit cycle, which invariably will have noncontractive 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 wellconditioned 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 turbomachinery will be presented.
 12^{30}  14^{00} Lunch
 14^{00} Algorithmic differentiation
 Raazesh Sainudiin (University of Canterbury, New Zealand)
An autovalidating rejection sampler for differentiable arithmetical expressions using randomized priority queues of mapped regular pavings
We introduce an efficient extension of a recently introduced
autovalidating 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 subboxes using a binary treebased data structure. The sequence of splits that
generate such a partition is given by the subboxes associated with the nodes of the tree.
The leaf boxes, i.e., the subboxes 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.
 Johannes Lotz (RWTH Aachen)
An Integer Programming Approach to the Call Tree Reversal Problem
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 NPcomplete. 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 NPcomplete 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.
 Shahadat Hossain (University of Lethbridge)
Pattern graph for sparse Jacobian matrix determination
Largescale 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.
 Shaun Forth (Cranfield University)
Plans for AD2016
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.

