Differenziazione automatica

La differenziazione automatica (in lingua inglese automatic differentiation, AD), nota anche come differenziazione algoritmica o differenziazione computazionale,^[1] è un insieme di tecniche per il calcolo automatico delle derivate di una funzione matematica implementata da un programma informatico. La differenziazione automatica sfrutta il fatto che l'implementazione di una funzione, indipendentemente da quanto sia complessa, si riduce all'esecuzione di una serie di operazioni aritmetiche (somma, sottrazione, moltiplicazione, divisione, etc.) e funzioni elementari (esponenziale, funzioni trigonometriche, etc.). Applicando la regola della catena ripetutamente a tali operazioni, la derivata di una funzione arbitrariamente complessa può essere calcolata automaticamente, alla precisione di calcolo in uso, e velocemente, usando un numero di operazioni equivalente al più ad un fattore piccolo e costante rispetto al numero di operazioni usate nella funzione originale.

La differenziazione automatica non deve essere confusa con la differenziazione simbolica, che viene eseguita manipolando le espressioni rappresentate in forma simbolica, né con la differenziazione numerica, che approssima la derivata con una differenza finita. Il primo metodo è solitamente molto più lento e limitato dalla difficoltà di convertire funzioni arbitrarie in forma simbolica, mentre il secondo metodo introduce errori numerici di discretizzazione che tendono a propagarsi e limitano la stabilità numerica. Inoltre, entrambi i metodi sono lenti nel calcolare derivate parziali di funzioni con un numero elevato di variabili, problema comunemente incontrato in metodi di ottimizzazione basati su discesa del gradiente. Tali metodi sono alla base dell'implementazione di molte tecniche di apprendimento automatico, in particolare apprendimento profondo, e la maggior parte dei framework di apprendimento profondo implementano la differenziazione automatica per il calcolo automatico del gradiente.^[2][3][4]

La regola della catena permette di scomporre una derivata in una serie di derivate più semplici. Per una composizione ${\displaystyle y=f(g(h(x)))=f(g(h(w_{0})))=f(g(w_{1}))=f(w_{2})=w_{3}}$ la regola della catena fornisce come risultato

${\displaystyle {\frac {dy}{dx}}={\frac {dy}{dw_{2}}}{\frac {dw_{2}}{dw_{1}}}{\frac {dw_{1}}{dx}}}$

La differenziazione automatica può essere operata in due modi distinti, accumulazione in avanti (forward accumulation o forward mode) e accumulazione inversa (reverse accumulation o reverse mode). L'accumulazione in avanti attraversa la catena di operazioni dall'interno all'esterno (nell'esempio, calcolando prima ${\displaystyle dw_{1}/dx}$, poi ${\displaystyle dw_{2}/dx}$ e infine ${\displaystyle dy/dx}$), mentre l'accumulazione inversa prevede di attraversare la catena dall'esterno all'interno (calcolando prima ${\displaystyle dy/dw_{2}}$, poi ${\displaystyle dy/dw_{1}}$ e infine ${\displaystyle dy/dx}$). In sintesi, si ha che l'accumulazione in avanti calcola la relazione ricorsiva

${\displaystyle {\frac {dw_{i}}{dx}}={\frac {dw_{i}}{dw_{i-1}}}{\frac {dw_{i-1}}{dx}}}$

con ${\displaystyle w_{3}=y}$, mentre l'accumulazione inversa calcola

${\displaystyle {\frac {dy}{dw_{i}}}={\frac {dy}{dw_{i+1}}}{\frac {dw_{i+1}}{dw_{i}}}}$

con ${\displaystyle w_{0}=x}$. In generale, entrambi i metodi sono casi particolari dell'applicazione dell'operatore di composizione dei programmi, fissando opportunamente uno dei due argomenti${\displaystyle (w,y)}$.

Nell'accumulazione in avanti, vengono fissate le variabili indipendenti nelle quali la differenziazione è calcolata, quindi il calcolo delle derivate delle sotto-espressioni è eseguito ricorsivamente. Questo può essere eseguito sostituendo ripetutamente le derivate delle funzioni interne nella regola della catena:

${\displaystyle {\frac {\partial y}{\partial x}}={\frac {\partial y}{\partial w_{n-1}}}{\frac {\partial w_{n-1}}{\partial x}}={\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)={\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)=\cdots }$

Il procedimento viene generalizzato al caso in più variabili, come prodotto di matrici jacobiane.

A differenza dell'accumulazione inversa, l'accumulazione in avanti è più naturale e semplice da implementare, in quanto il flusso di informazioni nel calcolo della derivata coincide con l'ordine di valutazione, aumentando ogni variabile ${\displaystyle w}$ con il valore numerico ${\displaystyle {\dot {w}}}$ della sua derivata,

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

Le derivate sono quindi calcolate in sincronia con i passaggi di valutazione delle espressioni, e combinate con altre derivate tramite la regola della catena.

Per esempio, sia data la funzione

${\displaystyle {\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}}}$

dove, per chiarezza, ogni singola sotto-espressione è stata etichettata con una variabile ${\displaystyle w_{i}}$. La scelta delle variabili indipendenti determina i valori iniziali ${\displaystyle {\dot {w}}_{1}}$ e ${\displaystyle {\dot {w}}_{2}}$ (chiamati semi). Ad esempio, se si vuole calcolare la derivata della funzione rispetto a ${\displaystyle x_{1}}$, i valori iniziali vengono impostati a

${\displaystyle {\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}}}$

Con questi semi, le derivate vengono propagate tramite la regola della catena, come mostrato nella tabella seguente e nel grafo computazionale nella figura a lato.

${\displaystyle {\begin{array}{l|l}{\text{Operazioni per calcolare l'espressione}}&{\text{Operazioni per calcolare la derivata}}\\\hline w_{1}=x_{1}&{\dot {w}}_{1}=1{\text{ (seme)}}\\w_{2}=x_{2}&{\dot {w}}_{2}=0{\text{ (seme)}}\\w_{3}=w_{1}\cdot w_{2}&{\dot {w}}_{3}=w_{2}\cdot {\dot {w}}_{1}+w_{1}\cdot {\dot {w}}_{2}\\w_{4}=\sin w_{1}&{\dot {w}}_{4}=\cos w_{1}\cdot {\dot {w}}_{1}\\w_{5}=w_{3}+w_{4}&{\dot {w}}_{5}={\dot {w}}_{3}+{\dot {w}}_{4}\end{array}}}$

Per calcolare il gradiente della funzione, sono necessari attraversamenti addizionali del grafo computazionale, usando semi adeguati, ad esempio ${\displaystyle {\dot {w}}_{1}=0;{\dot {w}}_{2}=1}$ per calcolare la derivata rispetto a ${\displaystyle x_{2}}$.

La complessità computazionale di ogni attraversamento del grafo è proporzionale alla complessità della funzione iniziale. L'accumulazione in avanti è più efficiente rispetto a quella inversa per funzioni nella forma ${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} ^{m}}$ con ${\displaystyle m\gg n}$, in quanto sono necessari solamente ${\displaystyle n}$ attraversamenti, rispetto agli ${\displaystyle m}$ richiesti dall'accumulazione inversa.

Nell'accumulazione inversa, vengono fissate le variabili dipendenti rispetto alle quali la differenziazione è calcolata, quindi la derivata viene calcolata ricorsivamente rispetto ad ogni sotto-espressione. Equivalentemente, la derivata della funzione più esterna nella regola della catena viene sostituita

${\displaystyle {\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 }$

La quantità usata nel calcolo tramite accumulazione inversa è chiamata aggiunta (adjoint), denotata con una barra (${\displaystyle {\bar {w}}}$), ed è la derivata di una variabile dipendente rispetto alla sotto-espressione ${\displaystyle w}$:

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

L'accumulazione inversa attraversa la catena dall'esterno verso l'interno, ovvero attraversa il grafo computazionale (figura a lato) dall'alto verso il basso. È necessario un attraversamento del grafo per ogni componente della funzione, quindi solo uno nel caso di funzioni scalari. Questo genera un compromesso tra costo in spazio e in tempo rispetto all'accumulazione in avanti: mentre l'accumulazione inversa può richiedere un numero inferiore di attraversamenti del grafo, è però necessario memorizzare le variabili intermedie ${\displaystyle w_{i}}$ (tipicamente usando una lista di Wengert),^[5][6] che può rappresentare un costo insostenibile in termini di memoria nel caso il grafo computazionale sia molto grande. Tale problema è parzialmente mitigato memorizzando solo una parte dei valori intermedi, e ricostruendo i restanti tramite calcoli addizionali (metodo noto come checkpointing).

Le operazioni richieste per calcolare la derivata tramite accumulazione inversa in questo esempio sono:

${\displaystyle {\begin{array}{l}{\text{Operazioni per il calcolo della derivata}}\\\hline {\bar {w}}_{5}=1{\text{ (seme)}}\\{\bar {w}}_{4}={\bar {w}}_{5}\\{\bar {w}}_{3}={\bar {w}}_{5}\\{\bar {w}}_{2}={\bar {w}}_{3}\cdot w_{1}\\{\bar {w}}_{1}={\bar {w}}_{3}\cdot w_{2}+{\bar {w}}_{4}\cdot \cos w_{1}\end{array}}}$

Il grafo computazionale può essere manipolato per calcolare il gradiente dell'operazione originaria, aumentando ogni nodo con un nodo aggiunto, e con i vertici tra nodi aggiunti orientati in direzione opposta rispetto ai vertici tra nodi originali. Tali nodi aggiunti rappresentano la moltiplicazione per la derivata dell'espressione calcolata nel nodo originale. Per esempio, un nodo contenente la funzione ${\displaystyle y=f(x)}$ avrà l'espressione ${\displaystyle {\bar {x}}={\bar {y}}f(x)}$ nel corrispondente nodo aggiunto.

L'accumulazione inversa è più efficiente rispetto a quella in avanti per funzioni nella forma ${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} ^{m}}$ con ${\displaystyle m\ll n}$, in quanto richiede ${\displaystyle m}$ attraversamenti del grafo computazionale, rispetto ai ${\displaystyle n}$ attraversamenti richiesti dall'accumulazione in avanti.

L'accumulazione inversa è stata pubblicata per la prima volta nel 1970 da Seppo Linnainmaa nella sua tesi di M.Sc.,^[7][8][9] ed è comunemente usata per implementare la retropropagazione dell'errore nelle reti neurali artificiali.

Le accumulazioni in avanti e inversa sono solo due metodi estremi per attraversare la catena di derivate. Il problema del calcolo della matrice jacobiana con un numero minimo di operazioni è noto come problema dell'accumulazione ottimale (optimal Jacobian accumulation, OJA), ed è un problema NP-completo.^[10]

L'accumulazione in avanti può essere implementata estendendo l'algebra dei numeri reali, dove ad ogni numero è associata una componente addizionale nilpotente che rappresenta la derivata della funzione valutata in quel numero. Le operazioni aritmetiche possono essere estese a questa nuova algebra, nota come algebra dei numeri duali. Tale formulazione può essere ulteriormente generalizzata, in quella che è nota come teoria del calcolo operazionale.

L'aritmetica estesa può essere costruita sostituendo ogni numero reale ${\displaystyle x}$ con la quantità ${\displaystyle x+x'\varepsilon }$, dove ${\displaystyle x'}$ è un numero reale e ${\displaystyle \varepsilon }$ è un infinitesimo, un numero astratto nilpotente tale che ${\displaystyle \varepsilon ^{2}=0}$. Da tale definizione, si ottiene che

${\displaystyle {\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}}}$

e le altre operazioni possono essere definite analogamente.

È possibile definire polinomi in questa aritmetica aumentata. Se ${\displaystyle P(x)=p_{0}+p_{1}x+p_{2}x^{2}+\cdots +p_{n}x^{n}}$, allora ${\displaystyle {\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}}}$

dove ${\displaystyle P^{(1)}}$ è la derivata di ${\displaystyle P}$ rispetto al primo argomento, e ${\displaystyle x'}$, chiamato seme (seed), può essere assegnato arbitrariamente.

Tale nuova aritmetica consiste di coppie ordinate ${\displaystyle \langle x,x'\rangle }$, con aritmetica ordinaria nella prima componente e aritmetica della derivata prima nella seconda componente. Estendendo i risultati precedenti dai polinomi alle funzioni analitiche, si ottiene una nuova aritmetica elementare:

${\displaystyle {\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}}}$

e in generale, per una funzione ${\displaystyle g}$ si ha

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

dove ${\displaystyle g_{u}}$ e ${\displaystyle g_{v}}$ sono le derivate di ${\displaystyle g}$ rispetto al primo e al secondo argomento rispettivamente.

Un'operazione binaria elementare può essere estesa a costanti, e l'applicazione a un numero duale ${\displaystyle \langle u,u'\rangle }$ e una costante (numero reale) ${\displaystyle c}$ è eseguita promuovendo la costante al numero duale ${\displaystyle \langle c,0\rangle }$. È possibile calcolare la derivata di una funzione ${\displaystyle f:\mathbb {R} \rightarrow \mathbb {R} }$ in ${\displaystyle x_{0}}$, calcolando ${\displaystyle f(\langle x_{0},1\rangle )}$ in questa aritmetica, il cui risultato è ${\displaystyle \langle f(x_{0}),f'(x_{0})\rangle }$.

Funzioni vettoriali a più variabili possono essere trattate con analoga efficienza usando lo stesso meccanismo per le funzioni a una variabile, definendo un operatore di derivata direzionale che calcola ${\displaystyle y'=\nabla f(x)\cdot x'}$, la derivata direzionale ${\displaystyle y'\in \mathbb {R} ^{m}}$ di ${\displaystyle f:\mathbb {R} ^{n}\rightarrow \mathbb {R} ^{m}}$ at ${\displaystyle x\in \mathbb {R} ^{n}}$ in direzione ${\displaystyle x'\in \mathbb {R} ^{n}}$, che può essere calcolata come ${\displaystyle (\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 )}$ usando la stessa aritmetica definita in precedenza. Per calcolare il gradiente ${\displaystyle \nabla f}$ sono necessarie ${\displaystyle n}$ valutazioni.

L'aritmetica definita per la derivata prima può essere estesa per calcolare derivate di ordine superiore. Tuttavia, all'aumentare dell'ordine le espressioni diventano molto più complesse, e il costo è quadratico nell'ordine della derivata più alta. È possibile definire una simile algebra su polinomi di Taylor troncati, consentendo di eseguire i calcoli più velocemente. Una formulazione generale e rigorosa è definita nel contesto del calcolo operazionale.

Lo stesso argomento in dettaglio: Calcolo operazionale.

Il calcolo operazionale nello spazio dei programmi^[11] generalizza il concetto di differenziazione automatica, generalizzando l'algebra dei numeri duali con un'algebra tensoriale, e fornisce una fondazione rigorosa alla matematica usata nell'apprendimento profondo.

L'accumulazione in avanti può essere implementata tramite un'interpretazione non-standard del programma, nella quale ogni numero reale è sostituito con un numero duale, le costanti sono sostituite da numeri duali con parte infinitesima nulla, e le operazioni sui reali sono sostituite con operazioni sui duali. L'interpretazione non-standard può essere generalmente implementata con due approcci: trasformazione del codice sorgente, e overloading degli operatori.

Il codice sorgente di una funzione è sostituito con il codice, generato automaticamente, di una funzione che interlaccia le istruzioni originarie con le istruzioni aggiuntive per il calcolo delle derivate. Tale approccio può essere implementato in ogni linguaggio, e rende possibile un maggior numero di ottimizzazioni automatiche da parte del compilatore, tuttavia l'implementazione dello strumento di generazione automatica è abbastanza complesso.

Tale metodo può essere usato in linguaggi che supportano l'overloading degli operatori. Non richiede generalmente di alterare la sequenza di operazioni nel programma, ma richiede la definizione di tipi di dati numerici aggiuntivi. È più semplice da implementare, e può essere usato anche per l'accumulazione inversa, tuttavia rende più complessa l'ottimizzazione automatica del codice. Esempi di librerie di differenziazione automatica tramite overloading degli operatori sono Adept e Stan.

• Louis B. Rall, Automatic Differentiation: Techniques and Applications, Lecture Notes in Computer Science, vol. 120, Springer, 1981, ISBN 3-540-10861-0.
• Andreas Griewank e Andrea Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Other Titles in Applied Mathematics, vol. 105, 2nd, SIAM, 2008, ISBN 978-0-89871-659-7 (archiviato dall'url originale il 23 marzo 2010).
• Richard Neidinger, Introduction to Automatic Differentiation and MATLAB Object-Oriented Programming (PDF), in SIAM Review, vol. 52, n. 3, 2010, pp. 545–563, DOI:10.1137/080743627. URL consultato il 15 marzo 2013.

• (EN) www.autodiff.org
• (EN) Automatic Differentiation of Parallel OpenMP Programs
• (EN) Automatic Differentiation, C++ Templates and Photogrammetry
• (EN) Compute analytic derivatives of any Fortran77, Fortran95, or C program through a web-based interface
• (EN) JuliaDiff, strumenti di differenziazione automatica in Julia

Portale Informatica

Portale Matematica