Programme of the Twelfth Euro AD Workshop

Three alternative methods for high-order multivariate AD were programmed in the same language (MATLAB) on the same computer (laptop) in order to compare accuracy and efficiency. In limited testing, the Direct Forward method was seen to be more accurate than two interpolation methods, the original from [GUW] and a method using nested directions [N]. Run-time measures sometimes confirmed the theoretical savings of the interpolation methods, though program structures can overshadow operation counts. In addition, the Direct Forward and [GUW] methods require one-time computation to set up the method for a specific number of variables and order of derivatives. Reasonable time for generating these reference arrays requires efficiency tricks that are not in most publications. Such a trick for the direct method will be shared.

McCormick introduced factorable programming in a technical report from 1967 which was later published in the book "Nonlinear Programming: Sequential Unconstrained Minimization Techniques" by Fiacco and McCormick from 1968. This is one of the first attempt to define structure in optimization as a tool to make more efficient implementations and new methods. In this talk we look closer at factorable programming and argue its equivalence to AD. We also look at some other classical 'structures in optimization' and the use of structure.

Allowing besides smooth elementals the functions abs, min, and max one obtains a class of piecewise differentiable functions, about which a fair amount of mathematical theory is known. They can be approximated to second order truncation errors by piecewise linearizations. We show how these may be evaluated constructively in an AD fashion and how they can be used to iteratively minimize nonsmooth functions and solve nonsmooth systems of equations.

For computing adjoints of codes including an iterative matrix-free algorithm for solving a nonlinear system (e.g. Newton/gmres), we want to give for discussion three different approaches. It is important to mantain the algorithm to be matrix-free as well as to tackle the memory consumption if trying to adjoin the algorithm black-box. The alternatives are using semi-direct (semi-analytic) approaches.

Many applications require permuting an input vector or intermediate results in a given way in order to evaluate a function, however derivatives of the result are required w.r.t. the original independents. Sometimes the permutation itself may be computed within the function evaluation and has no relevance outside of this computation. Similarly selection of certain independent or dependent variables based on some conditions may be required. Such operations may be categorised as branching in program code, but are infact simpler to handle in an AD-Tool like ADOL-C than general if-then-else branches that require retracing the function. In this talk I shall show how ADOL-C can handle such situations without the need for retracing the function.

For a fixed given integrand function f, good weights are required for construction of lattice rules for integration in high--dimensions. The weights define the embedding of f on a precise weighted unanchored Sobolev space. 'Optimal weights' are chosen to minimize certain upper bounds for the worst case integration error over these spaces. We use the new algorithmic differentiation technique ("Cross-Differentiation") for a low cost computation of 'optimal weights', based on the low--effective dimension of the integrands in some examples related to finance.

In this presentation I will introduce the Taylor series method for ordinary differential equations. It works by calculating the Taylor coefficients of the solution up to arbitrary order and by using these coefficients to compute the next step of the solution. Due to its high order, the method is especially suited for problems where very high precision is needed, like in astrophysical applications or simulations of chaotic dynamical systems. Along with the method a very efficient computation scheme is introduced which is based on auto-differentiation. This scheme exploits the power of the C++ template system and is uses C++ meta-programming and expression templates. Contrary to existing implementation of the Taylor series method, an additional preprocessing step is not necessary because the ODE can be completely expressed in terms of C++ templates. Furthermore, the performance of this method is as good as specialized code for an individual ODE.

Lie derivatives are derivatives of tensor fields along the solution of a given dynamic system. They occur in many fields such as general relativity, non-relativistic mechanics and control theory. We are working toward an efficient implementation of the calculation of Lie derivatives of scalar, vector and covector fields based on the low-level C drivers of ADOL-C.

TBA

Given a function, or more generally an operator $F: u \mapsto F(u)$, the question ``Is $F$ linear?" seems simple. In many applications of scientific computing it might be worth determining the answer to that question in an automated way; some functionality is only defined for linear operators and in other problems, time saving is available if it is known that the problem being solved is linear. However, implementing such an automated detection is not as straightforward as one might expect. This talk describes how linearity detection is implemented within the Chebfun system. The key ingredient is that information obtained via automatic differentiation is propagated as derivatives are evaluated.

In adjoint mode AD, the Call Tree Reversal problem is known to be NP-complete (U. Naumann. Call Tree Reversal is NP-complete. In C. Bischof et. al., Advances in Automatic Differentiation, pages 13–22. Springer, 2008). Based on our overloading AD tool dco/c++ (derivative code by overloading for c++), we discuss a method for implementing arbitrary (optimized) Call Tree Reversal schemes in a semi automatic fashion.

This project extends the ADOL-C package such that the differentiation of MPI-parallel simulation codes becomes possible. Using MPI the user explicitly decides which calculation is done by which process. Currently the simulation to be differentiated with ADOL-C can include the MPI routines Send, Receive, Barrier, Broadcast, Scatter, Gather and Reduce. For the implementation of these routines each MPI-function has its own wrapper containing appropriate algorithms for sending/receiving data using the active variables. For the subsequent derivative calculation, additional send- and receive-type operations were added to take different modes like first-order-scalar, higher-order-vector, etc. into account. Furthermore, to perform data transactions in forward or reverse mode, strategies were included to take the corresponding order into account.

The automatic differentiation of a huge C++ package, using the AD tool DCO, with many dependencies and underlying libraries has its own challenges. We will report about the problems encountered during the automatic differentiation and how we could solve them. Finally, we present the verification of the calculated gradients and give an overview of the amount of work it took to differentiate the code.

For efficient detailed aerodynamic design as well as optimal unsteady flow control, the use of adjoint approaches is a first essential ingredient. We compare continuous and discrete adjoint approaches. For the generation of discrete adjoint solvers, we discuss the use of Automatic Differentiation (AD) and for unsteady optimal flow control its combination with checkpointing techniques. In the case of detailed aerodynamic shape optimization, we discuss so-called one-shot methods. Here, one achieves simultaneously convergence of the primal state equation, the adjoint state equation as well as the design equation. The direction and size of the one-shot optimization steps are determined by a carefully selected design space preconditioner. It turns out, that one-shot methods enable aerodynamic shape designs for the computational effort of a small, constant multiple of the effort of an aerodynamic simulation. Integral part of these approaches are, next to suitable preconditioners for the coupled one-shot loop, gradient smoothing and shape derivatives.

Ned Nedialkov and I are the authors of the C++ code DAETS for solving DAE initial value problems, using Pryce's structural analysis (SA) approach, and a Taylor series expansion of the solution at each time-step.

Ned is talking about DAESA, which is essentially the SA front-end of DAETS, rewritten as a stand-alone Matlab tool with extra features to help exploratory analysis.

My talk is about the synergy, during the SA, between block triangularization and the detection of quasilinearity. Our algorithms determine a block-lower-triangular (BLT) sequence of solving for variables and their derivatives in the correct order. This is done as part of consistent initialization, and also at each integration step to project the solution on the consistent manifold.

The equations concerned are nonlinear in general, and a mix of square (as many equations as unknowns) and under-constrained (fewer equations than unknowns). Our semi-symbolic processing lets us determine whether the unknowns to solve for, in some block subsystem, actually occur linearly: we call this quasilinearity.

When a block is both linear and square, it needs no initial values or trial guesses for solution. In many of the test examples we have tried, this reduces the amount of initial data the user must submit, and has the potential to seriously speed up solution. For the Chemical Akzo problem for instance, an apparently nonlinear 8x8 system reduces to eight linear 1x1 blocks that are trivial to solve.

The quasilinearity analysis currently in DAETS is fairly primitive. We are putting an improved version into DAESA before including it in DAETS.

I show what improvements to quasilinearity analysis seem desirable and practical, with a number of examples.

A regularity region describes the local characteristics of a Differential-Algebraic Equation (DAE). We determine regularity regions of DAEs by means of sequences of continuous matrix functions. The matrix sequence is built step-by-step by certains admissible projector functions starting with the Jacobian matrices of the DAE data. For time-dependent and nonlinear DAEs the sequence contains a differentiation of a projector function. Beside common linear algebra tools such as matrix factorizations and generalized inverses, widely orthogonal projector functions are applied and algorithmic differentiation technices are used to compute the needed Jacobian matrices and to realize the differentiation. A comparison with other approaches, the discussion of numerical experiments and open problems completes the paper.

Differential-algebraic equations (DAEs) present today the state-of-the-art in dynamical systems arising from automated modularized modeling in almost all areas of science and engineering. While the modeling becomes more and more convenient, the resulting models are typically not easy to treat with current numerical simulation, control and optimization methods. In many cases a reformulation of the models or even a regularization is necessary to avoid failure of the computational methods. In this contribution we will discuss general DAE control problems and how (based on derivative arrays and smooth factorizations) they can be systematically reformulated and regularized so that the resulting system can be used in simulation, control and optimization procedures without much further difficulty.

Ned Nedialkov and John Pryce are the authors of DAETS, a C++ code for solving high-index, any order differential-algebraic equations (DAEs). It is based on Pryce’s structural analysis (SA) theory and a Taylor series expansion of the solution. Although DAETS can be used for SA, our goal is to provide a "lighter", stand-alone tool for SA of DAEs. In the last few years, we have been using Matlab to investigate the structure of numerous DAE problems, which resulted in the DAESA package.

Using DAESA, a user can specify a DAE in a general form: it can be of high-index, fully implicit, and contain derivatives of order higher than one. DAESA pre-processes the DAE through operator overloading to determine what variables and derivatives of them occur in each equation, and to determine if the problem is linear in highest order derivatives. Then DAESA finds the structural index of a DAE, its degrees of freedom, constraints, and what variables and derivatives need to be initialized. It also constructs a block-triangular structure of the problem (where such exists), and prescribes a solution scheme exploiting this structure.

Joint work in John Pryce and Guangning Tan.

The calculation of derivatives appears at different levels in numerical solution methods for DAE optimal control problems ranging from the computation of iteration matrices in BDF methods for coupled DAE-PDE systems to checkpointing strategies in semi-smooth Newton methods. Some applications of AD are discussed.

The rising demand in the simulation and optimization of complex processes in different fields (e.g. automotive industry, energy supplying systems, electronic industry, medicine) yields more and more often the task to solve differential-algebraic equations (DAEs) in an efficient and robust way. Solving differential-algebraic equations (DAEs) numerically, one is confronted with integration tasks but also with differentiation tasks. These differentiations may cause substantial stability problems for numerical methods applied to DAEs. We will explain these stability problems and discuss some strategies to avoid them. We will see that the model formulation may play a significant role. Furthermore, we will address the problem that the inherent differentiation tasks usually are not given explicitly.

Our AD tool for Matlab and Octave, ADiMat, has been provided with a new high level user interface. Derivatives can now be computed in forward and reverse mode by a single command, specifying a function handle and arguments and the seed matrix. Matlab and Octave can both solve index-1 DAEs using so called direct methods. Matlab provides the integrator ode15s to that end and Octave has interfaces to the integrators DASPK and DASSL. We present the new ADiMat user interface and give examples of how ADiMat supplied derivatives can be used with these integrators.

Differential-algebraic equations (DAEs) are mainly used for the modeling of dynamical processes. For instance, the dynamical behavior of mechanical systems, electrical circuits, chemical reactions and many other are often described by DAEs, in particular, of quasi-linear structure. It is well known that the numerical treatment of DAEs is nontrivial in general and more complicated than the one of ordinary differential equations. Arising effects in the numerical treatment of DAEs are for example drift, instabilities, convergence problems, or inconsistencies. In this talk we will discuss quasi-linear DAEs of higher index, with respect to their regularization in view of an efficient and robust numerical simulation. We will present two iterative procedures which provides a general tool for the regularization of quasi-linear DAEs of an arbitrary index. One of this procedures is of rank inflation type while the other one is of rank deflation type. Both procedures regularize a quasi-linear DAE by lowering the index while all constraints are maintained, in particular, the hidden constraints. The procedures end with the projected-strangeness-free form of the quasi-linear DAE which can be used as basis for numerical simulations or further numerical investigations of the dynamical system.

Of all the DAE examples in the CWI test suite of stiff ODE and DAE, only the transistor amplifier and the ring modulator prove to be non-trivial under the structural analysis of Pryce/Nedialkov. The talk will discuss methods to minimally modify those DAE systems so that structural analysis leads to a valid result. The transistor amplifier is transformed to a regular form by identification of common sub-expressions. For the ring modulator, further study of the computational graph is necessary.