3.1 Algorithmic Differentiation

Participants: Laurent Hascoët, Valérie Pascual.

Algorithmic Differentiation (AD) differentiates programs. The input of AD is a source program $P$ that, given some $X\in {ℝ}^n$, returns some $Y=F\left(X\right)\in {ℝ}^m$, for a differentiable $F$. AD generates a new source program $P^{\text{'}}$ that, given $X$, computes some derivatives of $F$ 4.

Any execution of $P$ amounts to a sequence of instructions, which is identified with a composition of vector functions. Thus, if

where each $f_k$ is the elementary function implemented by instruction $I_k$. AD applies the chain rule to obtain derivatives of $F$. Calling $X_k$ the values of all variables after instruction $I_k$, i.e. $X_0=X$ and $X_k=f_k\left(X_{k-1}\right)$, the Jacobian of $F$ is

which can be mechanically written as a sequence of instructions $I_k^{\text{'}}$. This can be generalized to higher level derivatives, Taylor series, etc. Combining the $I_k^{\text{'}}$ with the control of $P$ yields $P^{\text{'}}$, and therefore this differentiation is piecewise.

The above computation of $F^{\text{'}}\left(X\right)$, albeit simple and mechanical, can be prohibitively expensive on large codes. In practice, many applications only need cheaper projections of $F^{\text{'}}\left(X\right)$ such as:

• Sensitivities, defined for a given direction $\dot{X}$ in the input space as:

  $$F^{\text{'}}\left(X\right).\dot{X}=f_p^{\text{'}}\left(X_{p-1}\right).f_{p-1}^{\text{'}}\left(X_{p-2}\right).\cdots .f_1^{\text{'}}\left(X_0\right).\dot{X}.$$ 3

  This expression is easily computed from right to left, interleaved with the original program instructions. This is the tangent mode of AD.

• Adjoints, defined after transposition ($F^{\text{'}*}$), for a given weighting $\overline{Y}$ of the outputs as:

  $$F^{\text{'}*}\left(X\right).\overline{Y}=f_1^{\text{'}*}\left(X_0\right).f_2^{\text{'}*}\left(X_1\right).\cdots .f_{p-1}^{\text{'}*}\left(X_{p-2}\right).f_p^{\text{'}*}\left(X_{p-1}\right).\overline{Y}.$$ 4

  This expression is most efficiently computed from right to left, because matrix$\times$vector products are cheaper than matrix$\times$matrix products. This is the adjoint mode of AD, most effective for optimization, data assimilation  25, adjoint problems  19, or inverse problems.

Adjoint AD builds a very efficient program  21, which computes the gradient in a time independent from the number of parameters $n$. In contrast, computing the same gradient with the tangent mode would require running the tangent differentiated program $n$ times.

However, the $X_k$ are required in the inverse of their computation order. If the original program overwrites a part of $X_k$, the differentiated program must restore $X_k$ before it is used by $f_{k+1}^{\text{'}*}\left(X_k\right)$. Therefore, the central research problem of adjoint AD is to make the $X_k$ available in reverse order at the cheapest cost, using strategies that combine storage, repeated forward computation from available previous values, or even inverted computation from available later values.

Another research issue is to make the AD model cope with the constant evolution of modern language constructs. From the old days of Fortran77, novelties include pointers and dynamic allocation, modularity, structured data types, objects, vectorial notation and parallel programming. We keep developing our models and tools to handle these new constructs.

3.2 Static Analysis and Transformation of programs

Participants: Laurent Hascoët, Valérie Pascual.

The most obvious example of a program transformation tool is certainly a compiler. Other examples are program translators, that go from one language or formalism to another, or optimizers, that transform a program to make it run better. AD is just one such transformation. These tools share the technological basis that lets them implement the sophisticated analyses  12 required. In particular there are common mathematical models to specify these analyses and analyze their properties.

An important principle is abstraction: the core of a compiler should not bother about syntactic details of the compiled program. The optimization and code generation phases must be independent from the particular input programming language. This is generally achieved using language-specific front-ends, language-independent middle-ends, and target-specific back-ends. In the middle-end, analysis can concentrate on the semantics of a reduced set of constructs. This analysis operates on an abstract representation of programs made of one call graph, whose nodes are themselves flow graphs whose nodes (basic blocks) contain abstract syntax trees for the individual atomic instructions. To each level are attached symbol tables, nested to capture scoping.

Static program analysis can be defined on this internal representation, which is largely language independent. The simplest analyses on trees can be specified with inference rules  15, 22, 13. But many data-flow analyses are more complex, and better defined on graphs than on trees. Since both call graphs and flow graphs may be cyclic, these global analyses will be solved iteratively. Abstract Interpretation  16 is a theoretical framework to study complexity and termination of these analyses.

Data flow analyses must be carefully designed to avoid or control combinatorial explosion. At the call graph level, they can run bottom-up or top-down, and they yield more accurate results when they take into account the different call sites of each procedure, which is called context sensitivity. At the flow graph level, they can run forwards or backwards, and yield more accurate results when they take into account only the possible execution flows resulting from possible control, which is called flow sensitivity.

Even then, data flow analyses are limited, because they are static and thus have very little knowledge of actual run-time values. Far before reaching the very theoretical limit of undecidability, one reaches practical limitations to how much information one can infer from programs that use arrays  28, 17 or pointers. Therefore, conservative over-approximations must be made, leading to derivative code less efficient than ideal.

3.3 Algorithmic Differentiation and Scientific Computing

Participants: Alain Dervieux, Laurent Hascoët, Bruno Koobus, Matthieu Gschwend, Stephen Wornom.

Scientific Computing provides reliable simulations of complex systems. For example it is possible to simulate the steady or unsteady 3D air flow around a plane that captures the physical phenomena of shocks and turbulence. Next comes optimization, one degree higher in complexity because it repeatedly simulates and applies gradient-based optimization steps until an optimum is reached. The next sophistication is robustness, that detects undesirable solutions which, although maybe optimal, are very sensitive to uncertainty on design parameters or on manufacturing tolerances. This makes second derivatives come into play. Similarly Uncertainty Quantification can use second derivatives to evaluate how uncertainty on the simulation inputs imply uncertainty on its outputs.

To obtain this gradient and possibly higher derivatives, we advocate adjoint AD (cf3.1) of the program that discretizes and solves the direct system. This gives the exact gradient of the discrete function computed by the program, which is quicker and more sound than differentiating the original mathematical equations  19. Theoretical results  18 guarantee convergence of these derivatives when the direct program converges. This approach is highly mechanizable. However, it requires careful study and special developments of the AD model  23, 26 to master possibly heavy memory usage. Among these additional developments, we promote in particular specialized AD models for Fixed-Point iterations  20, 14, efficient adjoints for linear algebra operators such as solvers, or exploitation of parallel properties of the adjoint code.

4.1 Algorithmic Differentiation

Algorithmic Differentiation of programs gives sensitivities or gradients, useful for instance for :

• optimum shape design under constraints, multidisciplinary optimization, and more generally any algorithm based on local linearization,
• inverse problems, such as parameter estimation and in particular 4Dvar data assimilation in climate sciences (meteorology, oceanography),
• first-order linearization of complex systems, or higher-order simulations, yielding reduced models for simulation of complex systems around a given state,
• adaptation of parameters for classification tools such as Machine Learning systems, in which Adjoint Differentiation is also known as backpropagation.
• mesh adaptation and mesh optimization with gradients or adjoints,
• equation solving with the Newton method,
• sensitivity analysis, propagation of truncation errors.

4.2 Multidisciplinary optimization

A CFD program computes the flow around a shape, starting from a number of inputs that define the shape and other parameters. On this flow one can define optimization criteria e.g. the lift of an aircraft. To optimize a criterion by a gradient descent, one needs the gradient of the criterion with respect to all inputs, and possibly additional gradients when there are constraints. Adjoint AD is the most efficient way to compute these gradients.

4.3 Inverse problems and Data Assimilation

Inverse problems aim at estimating the value of hidden parameters from other measurable values, that depend on the hidden parameters through a system of equations. For example, the hidden parameter might be the shape of the ocean floor, and the measurable values of the altitude and velocities of the surface. Figure 1 shows an example of an inverse problem using the glaciology code ALIF (a pure C version of ISSM  24) and its AD-adjoint produced by Tapenade.

One particular case of inverse problems is data assimilation  25 in weather forecasting or in oceanography. The quality of the initial state of the simulation conditions the quality of the prediction. But this initial state is not well known. Only some measurements at arbitrary places and times are available. A good initial state is found by solving a least squares problem between the measurements and a guessed initial state which itself must verify the equations of meteorology. This boils down to solving an adjoint problem, which can be done though AD  27. The special case of 4Dvar data assimilation is particularly challenging. The 4th dimension in “4D” is time, as available measurements are distributed over a given assimilation period. Therefore the least squares mechanism must be applied to a simulation over time that follows the time evolution model. This process gives a much better estimation of the initial state, because both position and time of measurements are taken into account. On the other hand, the adjoint problem involved is more complex, because it must run (backwards) over many time steps. This demanding application of AD justifies our efforts in reducing the runtime and memory costs of AD adjoint codes.

4.4 Linearization

Simulating a complex system often requires solving a system of Partial Differential Equations. This can be too expensive, in particular for real-time simulations. When one wants to simulate the reaction of this complex system to small perturbations around a fixed set of parameters, there is an efficient approximation: just suppose that the system is linear in a small neighborhood of the current set of parameters. The reaction of the system is thus approximated by a simple product of the variation of the parameters with the Jacobian matrix of the system. This Jacobian matrix can be obtained by AD. This is especially cheap when the Jacobian matrix is sparse. The simulation can be improved further by introducing higher-order derivatives, such as Taylor expansions, which can also be computed through AD. The result is often called a reduced model.

4.5 Mesh adaptation

Some approximation errors can be expressed by an adjoint state. Mesh adaptation can benefit from this. The classical optimization step can give an optimization direction not only for the control parameters, but also for the approximation parameters, and in particular the mesh geometry. The ultimate goal is to obtain optimal control parameters up to a precision prescribed in advance.

5.1 Impact of research results

Our research has an impact on environmental questions through two of its application domains

• in CFD, we provide gradients that are used for design and shape optimization, and in many cases the optimization criterion is to reduce fuel or energy consumption. This is the case in particular for automotive or aircraft design.
• in Earth sciences, our gradients are used in inverse problems, to determine key properties in oceanography, glaciology, or climate models. For instance they determine basal friction coefficients of glaciers that are necessary to simulate their future evolution. Another example is to locate sources and sinks of CO2 by coupling atmospheric models and remote measurements.

6.1 New software

7.1 Algorithmic Differentiation of OpenMP

Participants: Laurent Hascoët, Jan Hueckelheim.

An article has been written and is currently under review with ACM TOMS.

7.2 Application to large industrial codes

Participants: Valérie Pascual, Laurent Hascoët, Hervé Guillard, Bruno Maugars, Sébastien Bourasseau, Cédric Content.

On the academic side, we supported usage of Tapenade for adjoint differentiation of several large codes:

• The plasma code SOLPS-ITER, differentiated by Stefano Carli at KU Leuven.
• The plasma code CTFEM, differentiated by Ali Elarif of the INRIA team Castor (Supervisor Hervé Guillard)
• The glaciology code Sicopolis, differentiated by Shreyas Suni Gaikwad at U. Texas Austin (Supervisor Patrick Heimbach)

7.3 Aeroacoustics

Participants: Alain Dervieux, Matthieu Gschwend, Bruno Koobus, Florian Miralles, Stephen Wornom, Tanya Kozubskaya.

The four-year program Norma associating IMAG of Montpellier university (B. Koobus, coordinator), CAALAB of Keldysh Institute of Moscow (T. Kozubskaya) and our team Ecuador, is supported by ANR and Russian Science Foundation.

See https://imag.umontpellier.fr/~koobus/norma.html

Norma is a cooperation on the subject of the extension of these methods to simulate the noise emission of rotating machines (helicopters, future aerial vehicles, unmanned aerial vehicles, wind turbines...).

7.4 Turbulence models

Participants: Alain Dervieux, Bruno Koobus, Florian Miralles, Stephen Wornom, Tanya Kozubskaya.

The development of hybrid models, in particular DES in the litterature, has raised the question of the domain of validity of these models. According to theory, these models should not be applied to flow involving laminar boundary layers (BL). But industrial flows are complex flows and often present regions of laminar BL, regions of fully developed turbulent BL and regions of non-equilibrium vortical BL. It is then mandatory for industrial use that the new hybrid models give a reasonable prediction for all these types of flow. We concentrated on evaluating the behavior of hybrid models for laminar BL and for vortical wakes. While less predictive than pure LES on laminar BL, some hybrid models still give reasonable predictions for rather low Reynolds numbers.

During the first phase of Norma, Montpellier and Moscow are computing a series of initial test cases in order to control the consistancy of the results produced by the two platforms of CFD, namely Noisette for Moscow, and Aironum for Montpellier.

7.5 Rotating machines

Participants: Alain Dervieux, Didier Chargy, Matthieu Gschwend, Bruno Koobus, Florian Miralles, Tanya Kozubskaya.

7.6 High order approximations

Participants: Alain Dervieux, Matthieu Gschwend, Bruno Koobus, Stephen Wornom, Tanya Kozubskaya.

7.7 Control of approximation errors

Participants: Alain Dervieux, Matthieu Gschwend, Bruno Koobus, Adrien Loseille, Frédéric Alauzet.

This year, a novel a priori estimate has been developed for incompressible flows, in cooperation with Lemma. We are also starting a new analysis for h-p anisotropic mesh adaptation. The monography on mesh adaptation currently being written by Alauzet, Loseille, Koobus and Dervieux now has been extended to 16 chapters.

8.1 Promoting scientific activities

9.1 Major publications

• 1 articleF. Courty, A. Dervieux, B. Koobus and L. Hascoët. Reverse automatic differentiation for optimum design: from adjoint state assembly to gradient computationOptimization Methods and Software1852003, 615-627
• 2 inproceedings B. Dauvergne and L. Hascoët. The Data-Flow Equations of Checkpointing in reverse Automatic Differentiation International Conference on Computational Science, ICCS 2006, Reading, UK 2006
• 3 articleDanielD. Goldberg, Sri Hari KrishnaS. Narayanan, LaurentL. Hascoët and JeanJ. Utke. An optimized treatment for algorithmic differentiation of an important glaciological fixed-point problemGeoscientific Model Development952016, 27
  HAL
• 4 incollectionL. Hascoët. Adjoints by Automatic DifferentiationAdvanced data assimilation for geosciencesOxford University Press2014, URL: https://hal.inria.fr/hal-01109881
  back to text
• 5 article L. Hascoët, U. Naumann and V. Pascual. ``To Be Recorded'' Analysis in Reverse-Mode Automatic Differentiation Future Generation Computer Systems 21 8 2004
• 6 articleL. Hascoët and V. Pascual. The Tapenade Automatic Differentiation tool: Principles, Model, and SpecificationACM Transactions On Mathematical Software3932013, URL: http://dx.doi.org/10.1145/2450153.2450158
• 7 articleL. Hascoët, J. Utke and U. Naumann. Cheaper Adjoints by Reversing Address ComputationsScientific Programming1612008, 81--92
• 8 articleLaurentL. Hascoët and JeanJ. Utke. Programming language features, usage patterns, and the efficiency of generated adjoint codeOptimization Methods and Software312016, 885--903
  HALDOI
• 9 articleL. Hascoët, M. Vázquez, B. Koobus and A. Dervieux. A Framework for Adjoint-based Shape Design and Error ControlComputational Fluid Dynamics Journal1642008, 454-464
• 10 article Jan ChristianJ. Hueckelheim, LaurentL. Hascoët and Jens-DominikJ.-D. Müller. Algorithmic differentiation of code with multiple context-specific activities ACM Transactions on Mathematical Software 2016
  HAL

9.2 Publications of the year

9.3 Cited publications

• 12 book A. Aho, R. Sethi and J. Ullman. Compilers: Principles, Techniques and Tools Addison-Wesley 1986
  back to text
• 13 techreportI. Attali, V. Pascual and C. Roudet. A language and an integrated environment for program transformations3313INRIA1997, URL: http://hal.inria.fr/inria-00073376
  back to text
• 14 articleB. Christianson. Reverse accumulation and implicit functionsOptimization Methods and Software941998, 307--322
  back to text
• 15 articleD. Clément, J. Despeyroux, L. Hascoët and G. Kahn. Natural semantics on the computerProceedings, France-Japan AI and CS Symposium, ICOTAlso, Information Processing Society of Japan, Technical Memorandum PL-86-6. Also INRIA research report # 4161986, 49-89URL: http://hal.inria.fr/inria-00076140
  back to text
• 16 articleP. Cousot. Abstract InterpretationACM Computing Surveys2811996, 324-328
  back to text
• 17 articleB. Creusillet and F. Irigoin. Interprocedural Array Region AnalysesInternational Journal of Parallel Programming2461996, 513--546
  back to text
• 18 articleJ.C.J. Gilbert. Automatic differentiation and iterative processesOptimization Methods and Software11992, 13--21
  back to text
• 19 inproceedings M.-B. Giles. Adjoint methods for aeronautical design Proceedings of the ECCOMAS CFD Conference Swansea, U.K. 2001
  back to text back to text
• 20 articleA. Griewank and C. Faure. Reduced Gradients and Hessians from Fixed Point Iteration for State EquationsNumerical Algorithms30(2)2002, 113--139
  back to text
• 21 book A. Griewank and A. Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation SIAM, Other Titles in Applied Mathematics 2008
  back to text
• 22 phdthesis L. Hascoët. Transformations automatiques de spécifications sémantiques: application: Un vérificateur de types incremental Université de Nice Sophia-Antipolis 1987
  back to text
• 23 techreport P. Hovland, B. Mohammadi and C. Bischof. Automatic Differentiation of Navier-Stokes computations MCS-P687-0997 Argonne National Laboratory 1997
  back to text
• 24 articleE. Larour, J. Utke, B. Csatho, A. Schenk, H. Seroussi, M. Morlighem, E. Rignot, N. Schlegel and A. Khazendar. Inferred basal friction and surface mass balance of the Northeast Greenland Ice Stream using data assimilation of ICESat (Ice Cloud and land Elevation Satellite) surface altimetry and ISSM (Ice Sheet System Model)Cryosphere862014, 2335-2351URL: http://www.the-cryosphere.net/8/2335/2014/
  DOIback to text
• 25 articleF.-X. Le Dimet and O. Talagrand. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspectsTellus38A1986, 97-110
  back to textback to text
• 26 article B. Mohammadi. Practical application to fluid flows of automatic differentiation for design problems Von Karman Lecture Series 1997
  back to text
• 27 phdthesis N. Rostaing. Différentiation Automatique: application à un problème d'optimisation en météorologie université de Nice Sophia-Antipolis 1993
  back to text
• 28 inproceedings R. Rugina and M. Rinard. Symbolic Bounds Analysis of Pointers, Array Indices, and Accessed Memory Regions Proceedings of the ACM SIGPLAN'00 Conference on Programming Language Design and Implementation ACM 2000
  back to text