Dérivation automatique

En mathématique et en calcul formel, la dérivation automatique (DA), également appelé dérivation algorithmique, dérivation formelle, ou auto-dérivation est un ensemble de techniques d'évaluation de la dérivée d'une fonction par un programme informatique.

Généralités

La dérivation automatique exploite le fait que chaque programme informatique, aussi compliqué soit-il, exécute une séquence d'opérations arithmétiques élémentaires (addition, soustraction, multiplication, division, etc.) et des fonctions élémentaires (exp, log,sin, cos, etc.). En appliquant de façon répétée la règle de dérivation des fonctions composées (règle de la chaîne) à ces opérations, les dérivées d'ordre arbitraire peuvent être calculées automatiquement et avec autant de précision que souhaitée. De plus calculer la dérivation automatique de la fonction ne modifie son temps de calcul qu'au plus d'un petit facteur constant, rendant ces techniques très efficaces en pratique.

Figure 1 : Lien entre la dérivation automatique et la dérivation symbolique.

La dérivation automatique est distincte de la dérivation symbolique et de la dérivation numérique (méthode des différences finies). La dérivation symbolique peut conduire à un code inefficace et se heurte à la difficulté de convertir un programme informatique en une seule expression, tandis que la dérivation numérique peut introduire des erreurs d'arrondi lors de la discrétisation et de l'annulation. Ces deux méthodes classiques ont par ailleurs des problèmes lors du calcul de dérivées d'ordre plus élevées, car la complexité et/ou les erreurs augmentent. Enfin, ces deux méthodes sont lentes pour calculer les dérivées partielles d'une fonction lorsque cette dernière dépend de nombreuses variables, comme cela est nécessaire pour les algorithmes d'optimisation par descente de gradient. La dérivation automatique résout tous ces problèmes.

La règle de la chaîne : accumulation avant et arrière

La dérivation automatique repose fondamentalement sur la décomposition des dérivées fournie par la règle de la chaîne. Pour la composition simple

$${\begin{aligned}y&=f(g(h(x)))=f(g(h(w_{0})))=f(g(w_{1}))=f(w_{2})=w_{3}\\w_{0}&=x\\w_{1}&=h(w_{0})\\w_{2}&=g(w_{1})\\w_{3}&=f(w_{2})=y\end{aligned}}$$

la règle de la chaîne donne alors:

$${\frac {dy}{dx}}={\frac {dy}{dw_{2}}}{\frac {dw_{2}}{dw_{1}}}{\frac {dw_{1}}{dx}}={\frac {df(w_{2})}{dw_{2}}}{\frac {dg(w_{1})}{dw_{1}}}{\frac {dh(w_{0})}{dx}}$$

Habituellement, deux modes de DA distincts sont présentés, l'accumulation directe (ou mode avant) et l'accumulation inverse (ou mode inverse). L'accumulation directe spécifie que l'on parcourt la règle de la chaîne de l'intérieur vers l'extérieur (c'est-à-dire le premier calcul $dw_{1}/dx$ et puis $dw_{2}/dw_{1}$ et enfin $dy/dw_{2}$), tandis que l'accumulation inverse a la traversée de l'extérieur vers l'intérieur (premier calcul $dy/dw_{2}$ et puis $dw_{2}/dw_{1}$ et enfin $dw_{1}/dx$). Plus succinctement,

1. l'accumulation directe calcule la relation récursive : ${\frac {dw_{i}}{dx}}={\frac {dw_{i}}{dw_{i-1}}}{\frac {dw_{i-1}}{dx}}$ avec $w_{3}=y$, et,
2. l'accumulation inverse calcule la relation récursive : ${\frac {dy}{dw_{i}}}={\frac {dy}{dw_{i+1}}}{\frac {dw_{i+1}}{dw_{i}}}$ avec $w_{0}=x$.

Accumulation directe

Figure 2 : Exemple d'accumulation directe avec un graphe de calcul

Lors d'une DA par accumulation directe, on fixe d'abord le variable indépendante par rapport à laquelle la dérivation est effectuée et on calcule la dérivée de chaque sous-expression récursivement. Dans un calcul papier-crayon, cela implique de substituer à plusieurs reprises la dérivée de la fonction intérieure dans la règle de chaîne :

$${\begin{aligned}{\frac {\partial y}{\partial x}}&={\frac {\partial y}{\partial w_{n-1}}}{\frac {\partial w_{n-1}}{\partial x}}\\[6pt]&={\frac {\partial y}{\partial w_{n-1}}}\left({\frac {\partial w_{n-1}}{\partial w_{n-2}}}{\frac {\partial w_{n-2}}{\partial x}}\right)\\[6pt]&={\frac {\partial y}{\partial w_{n-1}}}\left({\frac {\partial w_{n-1}}{\partial w_{n-2}}}\left({\frac {\partial w_{n-2}}{\partial w_{n-3}}}{\frac {\partial w_{n-3}}{\partial x}}\right)\right)\\[6pt]&=\cdots \end{aligned}}$$

Cela peut être généralisé à plusieurs variables en considérant le produit de matrices jacobiennes.

Par rapport à l'accumulation inverse, l'accumulation directe est naturelle et facile à mettre en œuvre car le flux d'information des dérivées coïncide avec l'ordre d'évaluation. Chaque variable w est augmentée de sa dérivée ẇ (stockée sous forme de valeur numérique et non d'expression symbolique),

$${\dot {w}}={\frac {\partial w}{\partial x}}$$

comme indiqué par le point. Les dérivées sont ensuite calculées de façon synchronisée lors des étapes d'évaluation du programme et combinées avec les autres dérivées via la règle de chaîne.

À titre d'exemple, considérons la fonction :

$${\begin{aligned}z&=f(x_{1},x_{2})\\&=x_{1}x_{2}+\sin x_{1}\\&=w_{1}w_{2}+\sin w_{1}\\&=w_{3}+w_{4}\\&=w_{5}\end{aligned}}$$

Pour plus de clarté, les sous-expressions individuelles ont été étiquetées avec les variables w_i.

Le choix de la variable indépendante à laquelle la différenciation est effectuée affecte les valeurs initiales de ẇ₁ et ẇ₂. Si l'on s'intéresse à la dérivée de cette fonction par rapport à x₁, les valeurs de départ doivent être définies par :

$${\begin{aligned}{\dot {w}}_{1}={\frac {\partial x_{1}}{\partial x_{1}}}=1\\{\dot {w}}_{2}={\frac {\partial x_{2}}{\partial x_{1}}}=0\end{aligned}}$$

Une fois les valeurs de départ ainsi définies, les valeurs se propagent à l'aide de la règle de chaîne comme indiqué. La figure 2 montre une représentation picturale de ce processus sous forme de graphique de calcul.

Pour calculer le gradient de la fonction de cet exemple, il est nécessaire de calculer les dérivées de f par rapport à x₁ et x₂, un parcours supplémentaire est effectué sur le graphique de calcul avec les es valeurs de départ ${\dot {w}}_{1}=0;{\dot {w}}_{2}=1$.

La complexité informatique du calcul de DA lors d'un parcours de l'accumulation directe est proportionnelle à la complexité du code d'origine.

L'accumulation directe est plus efficace que l'accumulation inverse pour les fonctions f : Rⁿ → R^m avec m ≫ n puisque seul n parcours sont nécessaires, par rapport à m parcours pour l'accumulation inverse.

Accumulation inverse

Figure 3 : Exemple d'accumulation inverse avec son graphe de calcul.

Lors d'une DA par accumulation inverse, la variable dépendante à dériver est fixe et la dérivée est calculée récursivement pour chaque sous-expression . Dans un calcul papier-crayon, la dérivée des fonctions extérieures sont substituées à plusieurs reprises dans la règle de chaîne :

$${\frac {\partial y}{\partial x}}={\frac {\partial y}{\partial w_{1}}}{\frac {\partial w_{1}}{\partial x}}=\left({\frac {\partial y}{\partial w_{2}}}{\frac {\partial w_{2}}{\partial w_{1}}}\right){\frac {\partial w_{1}}{\partial x}}=\left(\left({\frac {\partial y}{\partial w_{3}}}{\frac {\partial w_{3}}{\partial w_{2}}}\right){\frac {\partial w_{2}}{\partial w_{1}}}\right){\frac {\partial w_{1}}{\partial x}}=\cdots$$

Dans l'accumulation inverse, la quantité d'intérêt est l'adjoint, noté avec une barre (w̄); c'est une dérivée d'une variable dépendante choisie par rapport à une sous-expression w :

$${\bar {w}}={\frac {\partial y}{\partial w}}$$

L'accumulation inverse traverse la règle de la chaîne de l'extérieur vers l'intérieur, ou dans le cas du graphique de calcul de la figure 3, de haut en bas. L'exemple de fonction est à valeur scalaire, et il n'y a donc qu'une seule graine pour le calcul de la dérivée, et un seul parcours du graphe de calcul est nécessaire pour calculer le gradient (qui a pourtant deux composantes). C'est seulement la moitié du travail par rapport à l'accumulation directe, mais l'accumulation inverse nécessite le stockage des variables intermédiaires w_i ainsi que les instructions qui les ont produites dans une structure de données connue sous le nom de liste de Wengert (ou "bande"), ce qui peut consommer une mémoire importante si le graphe de calcul est volumineux. Cela peut être atténué dans une certaine mesure en ne stockant qu'un sous-ensemble des variables intermédiaires, puis en reconstruisant les variables de travail nécessaires en répétant les évaluations, une technique connue sous le nom de rematérialisation. Des Points de Contrôle peuvent également être utilisés pour enregistrer les états intermédiaires.

Les opérations pour calculer la dérivée en utilisant l'accumulation inverse sont présentées dans le tableau ci-dessous (notez l'ordre inversé).

Le graphe de flux des données d'un calcul peut être manipulé pour calculer le gradient de son calcul d'origine. En ajoutant un nœud adjoint pour chaque nœud primaire, et en reliant par des liens adjoints parallèles aux bords primaires mais circulant dans le sens opposé. Les nœuds du graphe adjoint représentent la multiplication par les dérivées des fonctions calculées par les nœuds du primaire. Par exemple, l'addition dans le primaire provoque la sortie dans l'adjoint ; la sortie dans le primaire provoque l'ajout dans l'adjoint ; a unaire fonction y = f (x) dans les causes primaires x̄ = ȳ f ′(x) dans l'adjoint, etc.

L'accumulation inverse est plus efficace que l'accumulation directe pour les fonctions f : Rⁿ → R^m avec m ≪ n puisque seuls m parcours sont alors nécessaires (par rapport à n parcours pour l'accumulation directe).

La DA par accumulation inverse a été publiée pour la première fois en 1976 par Seppo Linnainmaa.

La rétropropagation des erreurs dans les perceptrons multicouches, une technique centrale dans l'apprentissage automatique, est un cas particulier de mode inverse de DA.

Au-delà de l'accumulation avant et arrière

Les accumulations avant et arrière ne sont que deux façons (extrêmes) de traverser la règle de la chaîne. Le problème du calcul d'un Jacobien complet de f : Rⁿ → R^m avec un nombre minimum d'opérations arithmétiques est connu comme le problème d'accumulation optimale du Jacobien, qui est NP-complet. Au cœur de cette preuve est l'idée que des dépendances algébriques peuvent exister entre les dérivées partielles locales qui étiquettent les arêtes du graphe. En particulier, deux étiquettes (ou plus) peuvent être reconnues comme égales. La complexité du problème reste ouverte si l'on suppose que toutes les étiquettes de bord sont uniques et algébriquement indépendantes.

Dérivation automatique à l'aide de nombres duaux

La dérivation automatique directe est réalisée en augmentant l'algèbre des nombres réels pour obtenir une nouvelle arithmétique. Une composante supplémentaire est ajouté à chaque nombre pour représenter la dérivée d'une fonction de ce nombre, et tous les opérateurs arithmétiques sont étendus pour l'algèbre augmentée. L'algèbre augmentée est l'algèbre des nombres duaux.

On obtient cette algèbre en remplaçant chaque nombre $\,x$ par le nombre $x+x'\varepsilon$, où $x'$ est un nombre réel, mais $\varepsilon$ est un nombre abstrait avec la propriété $\varepsilon ^{2}=0$ (un infinitésimal ; voir Analyse infinitésimale fluide). En ajoutant uniquement cela, l'arithmétique régulière donne alors :

$${\begin{aligned}(x+x'\varepsilon )+(y+y'\varepsilon )&=x+y+(x'+y')\varepsilon \\(x+x'\varepsilon )\cdot (y+y'\varepsilon )&=xy+xy'\varepsilon +yx'\varepsilon +x'y'\varepsilon ^{2}=xy+(xy'+yx')\varepsilon \end{aligned}}$$

et de même pour la soustraction et la division.

Maintenant, les polynôme peuvent être calculés dans cette arithmétique augmentée. Si $P(x)=p_{0}+p_{1}x+p_{2}x^{2}+\cdots +p_{n}x^{n}$, alors

$${\begin{aligned}P(x+x'\varepsilon )&=p_{0}+p_{1}(x+x'\varepsilon )+\cdots +p_{n}(x+x'\varepsilon )^{n}\\&=p_{0}+p_{1}x+\cdots +p_{n}x^{n}+p_{1}x'\varepsilon +2p_{2}xx'\varepsilon +\cdots +np_{n}x^{n-1}x'\varepsilon \\&=P(x)+P^{(1)}(x)x'\varepsilon \end{aligned}}$$

où $P^{(1)}$ désigne la dérivée de $P$ en ce qui concerne son premier argument, et $x'$, (appelé la graine, ou valeur initiale), peut être choisi arbitrairement.

Cette nouvelle arithmétique se compose de couples, éléments écrits $\langle x,x'\rangle$, avec l'arithmétique ordinaire sur la première composante, et l'arithmétique de différenciation du premier ordre sur la deuxième composante, comme décrit ci-dessus. Étendre les résultats ci-dessus sur les polynômes à fonctions analytiques donne une liste de l'arithmétique de base et de certaines fonctions standard pour la nouvelle arithmétique :

$${\begin{aligned}\left\langle u,u'\right\rangle +\left\langle v,v'\right\rangle &=\left\langle u+v,u'+v'\right\rangle \\\left\langle u,u'\right\rangle -\left\langle v,v'\right\rangle &=\left\langle u-v,u'-v'\right\rangle \\\left\langle u,u'\right\rangle *\left\langle v,v'\right\rangle &=\left\langle uv,u'v+uv'\right\rangle \\\left\langle u,u'\right\rangle /\left\langle v,v'\right\rangle &=\left\langle {\frac {u}{v}},{\frac {u'v-uv'}{v^{2}}}\right\rangle \quad (v\neq 0)\\\sin \left\langle u,u'\right\rangle &=\left\langle \sin(u),u'\cos(u)\right\rangle \\\cos \left\langle u,u'\right\rangle &=\left\langle \cos(u),-u'\sin(u)\right\rangle \\\exp \left\langle u,u'\right\rangle &=\left\langle \exp u,u'\exp u\right\rangle \\\log \left\langle u,u'\right\rangle &=\left\langle \log(u),u'/u\right\rangle \quad (u>0)\\\left\langle u,u'\right\rangle ^{k}&=\left\langle u^{k},ku^{k-1}u'\right\rangle \quad (u\neq 0)\\\left|\left\langle u,u'\right\rangle \right|&=\left\langle \left|u\right|,u'{\mbox{sign}}u\right\rangle \quad (u\neq 0)\end{aligned}}$$

et en général pour la fonction primitive $g$,

$$g(\langle u,u'\rangle ,\langle v,v'\rangle )=\langle g(u,v),g_{u}(u,v)u'+g_{v}(u,v)v'\rangle$$

où $g_{u}$ et $g_{v}$ sont les dérivés de $g$ en ce qui concerne ses premier et deuxième arguments, respectivement.

Lorsqu'une opération arithmétique binaire de base est appliquée à des arguments mixtes - la paire $\langle u,u'\rangle$ et le nombre réel $c$—le nombre réel est d'abord porté à $\langle c,0\rangle$. La dérivée d'une fonction $f:\mathbb {R} \to \mathbb {R}$ au point $x_{0}$ est maintenant trouvé en calculant $f(\langle x_{0},1\rangle )$ en utilisant l'arithmétique ci-dessus, ce qui donne pour résultat $\langle f(x_{0}),f'(x_{0})\rangle$.

Arguments et fonctions vectoriels

Les fonctions multivariées peuvent être gérées avec les mêmes mécanismes (et la même efficacité) que les fonctions univariées en adoptant un opérateur de dérivé directionnel. Autrement dit, il suffit de calculer $y'=\nabla f(x)\cdot x'$, la dérivée directionnelle $y'\in \mathbb {R} ^{m}$ de $f:\mathbb {R} ^{n}\to \mathbb {R} ^{m}$ à $x\in \mathbb {R} ^{n}$ dans la direction $x'\in \mathbb {R} ^{n}$, et cela peut être calculé comme $(\langle y_{1},y'_{1}\rangle ,\ldots ,\langle y_{m},y'_{m}\rangle )=f(\langle x_{1},x'_{1}\rangle ,\ldots ,\langle x_{n},x'_{n}\rangle )$ en utilisant la même arithmétique que ci-dessus. Si tous les éléments de $\nabla f$ sont désirés, alors $n$ évaluations de la fonction seront nécessaires. Notez que dans de nombreuses applications d'optimisation, la dérivée directionnelle est effectivement suffisante.

Ordre élevé et variables nombreuses

L'arithmétique ci-dessus peut être généralisée pour calculer les dérivées du second ordre (et supérieures) des fonctions multivariées. Cependant, les règles arithmétiques se compliquent rapidement : la complexité est quadratique au degré de dérivée le plus élevé. Au lieu de cela, l'algèbre tronquée des polynômes de Taylor peut être utilisée. L'arithmétique résultante, définie sur des nombres duaux généralisés, permet un calcul efficace en utilisant des fonctions comme s'il s'agissait d'un type de données. Une fois que le polynôme de Taylor d'une fonction est connu, les dérivées sont facilement extraites.

Implémentation

La DA en mode direct est généralement implémentée par une interprétation non standard du programme dans lequel les nombres réels sont remplacés par des nombres duaux, les constantes sont remplacées par des nombres doubles avec un coefficient epsilon nul, et les primitives numériques sont remplacées pour fonctionner sur des nombres duaux. Cette interprétation non standard est généralement mise en œuvre à l'aide de l'une des deux stratégies suivantes : transformation du code source ou surcharge de l'opérateur.

Transformation du code source (TCS)

Figure 4 : Exemple de fonctionnement de la transformation du code source

Le code source d'une fonction est remplacé par un code source généré automatiquement qui comprend des instructions de calcul des dérivées, entrelacés avec les instructions d'origine.

La transformation du code source peut être implémenté dans tous les langages de programmation, et elle rend également plus facile pour le compilateur l'optimisation du temps de compilation. Cependant, la mise en œuvre de la DA elle-même est plus difficile.

Surcharge d'opérateur (SO)

Figure 5 : Exemple de fonctionnement de la surcharge d'opérateur

La surcharge d'opérateur est une possibilité lorsque le code source est écrit dans un langage le permettant. Les objets pour les nombres réels et les opérations mathématiques élémentaires doivent être surchargés pour répondre à l'arithmétique augmentée décrite ci-dessus. Cela ne nécessite aucun changement dans la forme ou la séquence des opérations ni dans le code source d'origine pour que la fonction soit automatiquement différenciée. Mais cela nécessite souvent des changements dans les types de données de base pour les nombres et les vecteurs afin de supporter la surcharge et implique également souvent l'insertion d'opérations de marquage spéciales.

La surcharge d'opérateur pour l'accumulation directe est facile à mettre en œuvre, mais elle est également possible pour l'accumulation inverse. Pour cette dernière cependant, les compilateurs actuels sont en retard dans l'optimisation du code par rapport à l'accumulation directe.

La surcharge d'opérateur, pour l'accumulation avant et arrière, peut être bien adaptée aux applications où les objets sont des vecteurs de nombres réels plutôt que des scalaires. En effet, la bande comprend alors des opérations vectorielles ; et cela peut faciliter des implémentations efficaces sur le plan informatique où chaque opération vectorielle effectue de nombreuses opérations scalaires. Des techniques de dérivation algorithmique par adjoint vectoriel peuvent être utilisées, par exemple, pour dériver des valeurs calculées par simulation Monte-Carlo.

Des exemples d'implémentations de différenciation automatique par surcharge d'opérateurs en C++ sont disponibles dans les bibliothèques Adepte et Stan.

Notes et références

Annexes

Bibliographie

• Louis B. Rall, Automatic Differentiation: Techniques and Applications, vol. 120, Springer, coll. « Lecture Notes in Computer Science », 1981 (ISBN 978-3-540-10861-0)
• Andreas Griewank et Andrea Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, vol. 105, SIAM, coll. « Other Titles in Applied Mathematics », 2008 (ISBN 978-0-89871-659-7, lire en ligne)
• Neidinger, « Introduction to Automatic Differentiation and MATLAB Object-Oriented Programming », SIAM Review, vol. 52, n^o 3,‎ 2010, p. 545–563 (DOI 10.1137/080743627, lire en ligne, consulté le 15 mars 2013)
• Uwe Naumann, The Art of Differentiating Computer Programs, SIAM, coll. « Software-Environments-tools », 2012 (ISBN 978-1-611972-06-1)
• Marc Henrard, Algorithmic Differentiation in Finance Explained, Palgrave Macmillan, coll. « Financial Engineering Explained », 2017 (ISBN 978-3-319-53978-2)

Liens externes

• www.autodiff.org, Un « site d'entrée à tout ce que vous voulez savoir sur la différenciation automatique »
• Différenciation automatique des programmes OpenMP parallèles
• Différenciation automatique, modèles c++ et photogrammétrie
• Différenciation automatique, approche de surcharge de l'opérateur
• Calculez les dérivées analytiques de tout programme Fortran77, Fortran95 ou C via une interface Web ; différenciation automatique des programmes Fortran
• Description et exemple de code pour la différenciation automatique en avant dans Scala
• extensions de différenciation automatique finmath-lib Différenciation automatique pour les variables aléatoires (implémentation Java de la différenciation automatique stochastique)
• Différenciation algorithmique adjointe : calibration et théorème des fonctions implicites
• Article de différenciation automatique basé sur un modèle C++ et application
• Tangente Dérivés déboguables source à source
• Grecs exacts du premier et du second ordre par différenciation algorithmique
• Différenciation algorithmique adjointe d'une application accélérée par GPU
• Méthodes adjointes en finance informatique support d'outils logiciels pour la différenciation algorithmiqueop