**Diferenciación Automática**
La *diferenciación algorítmica* es el proceso de derivación automática de gradientes para algoritmos (programas) numéricos. Comúnmente conocida como **diferenciación automática** (**AD**), oor supuesto, no es el único método existente, y podemos encontrar las más conocidas, la **diferenciación numérica** y la **diferenciación simbólica**. Las diferencias más destacadas con estos, más clásicos, son:
1. La AD se deriva de la observación de que cualquier algoritmo numérico es finalmente una combinación de operaciones numéricas elementales que se pueden obtener por medio de la composición de funciones más simples por lo que, dadas las reglas de diferenciación elementales, las derivadas pueden propagarse numéricamente a través de las distintas fases del algoritmo utilizando la regla de la cadena.
2. La AD devuelve valores de gradiente exactos, no aproximados (como pasa con la diferenciación numérica), debido a la implementación de reglas de derivación de las funciones elementales involucradas.
3. La AD evita el engrosamiento de la expresión simbólica (como ocurre con la diferenciación simbólica) haciendo uso de cómputos numéricos intermedios (exactos) que reutiliza en múltiples ramas de cálculo.
Vamos a ver algunos detalles sobre estas afirmaciones:
Aunque la AD propaga valores numéricos, es claramente distinta de la diferenciación numérica. La diferenciación numérica utiliza algoritmos para estimar numéricamente la derivada de una función basándose en algunos valores de la función u otra información sobre ella. Por ejemplo, los algoritmos numéricos suelen basarse en la definición de la derivada de $f$ en el punto $x$: $$f'(x)=\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}$$
Sin embargo, como al establecer $h = 0$ se obtiene directamente la división por cero, es de esperar un error de estimación que, para valores pequeños de $h$, generan errores de redondeo en coma flotante que pueden convertirse en un problema de difícil solución.
La diferencia entre la diferenciación algorítmica y la simbólica es menos clara. La diferenciación simbólica manipula la expresión de la función para generar una expresión exacta de la función derivada. Se aplican las reglas fundamentales de la diferenciación, como la regla del producto, del logaritmo, o la regla de la cadena, y el resultado se simplifica hasta un nivel aceptable (cuando es posible). La diferenciación simbólica es exacta, pero es lenta cuando se aplica a una función de muchas variables de entrada, como es el caso de todas las aplicaciones de aprendizaje profundo. Pero quizás, la principal limitación es que la diferenciación simbólica requiere que el programa esté expresado en forma cerrada, es decir, como una expresión pura matemática, y no admite estructuras tales como los controles de flujo y bucles habituales en la programación, por lo que para aplicar la diferenciación simbólica primero habría que resolver la tarea de expresar todo el programa como una única expresión matemática. En cambio, el objetivo de la diferenciación algorítmica es hallar la derivada de un bloque arbitrario de código que calcula una función numérica. En el núcleo de ambas metodologías se encuentra la regla de la cadena de la diferenciación pero, a diferencia de la simbólica, la AD no pretende manipular las expresiones simbólicamente con esta regla, sino evaluar numéricamente los resultados intermedios y propagar el gradiente a través del programa. Con operaciones elementales que tienen reglas de diferenciación definidas, la AD es capaz de descomponer un bloque de código en sus operaciones básicas (muchas veces, expresadas explícitamente en el propio algoritmo), encontrar su gradiente y, a continuación, propagar los valores del gradiente con la regla de la cadena.
Podemos reconocer dos modos principales para la tarea de AD:
* **Modo directo**, que calcula las derivadas en paralelo con la ejecución del programa original, y
* **Modo inverso**, que propaga las derivadas a través de variables intermedias una vez completada la ejecución del programa original.
Podemos destacar algunas razones de peso para dedicar una entrada a un tema que, a priori, no sería más que una visión distinta del tema común de diferenciación en un curso básico de Análisis Matemático:
!!!side:1
Julia ofrece muchas implementaciones (algunas de ellas ya no se mantienen), y ha demostrado que implementar una AD (simple) es relativamente sencillo.
* Es el soporte de gran parte del aprendizaje automático moderno al permitir la diferenciación rápida de funciones matemáticas complejas. Los métodos de optimización de primer orden son omnipresentes en la búsqueda de parámetros de funciones (no solo en el aprendizaje profundo).
* Es interesante estudiar la AD tanto desde el punto de vista matemático como de la implementación, ya que los distintos enfoques conllevan diferentes ventajas y desventajas [1].
* Desde un punto de vista metodológico, es positivo entender (al menos aproximadamente), cómo funcionan los métodos AD con el fin de utilizarlos eficazmente en las tareas de ML que tan comunes y extendidas están hoy en día.
* Julia ofrece una opción única en el esfuerzo de separar las definiciones de las reglas AD de los motores AD que usan esas reglas y el backend que ejecuta las reglas. Esta capacidad facilita la tarea de añadir nuevas reglas que pueden ser compatibles con muchos frameworks de optimización sin la costosa tarea de crear interfaces o modificar los códigos previamente creados.
Pero, antes de comenzar a ver cómo funcionan estos mecanismos, veremos, aunque sea brevemente, algo acerca de los Grafos de Computación, que son las estructuras de datos funcionales necesarias para implementar cómodamente los algoritmos y procesos involucrados.
# Grafos de Computación
Los **Grafos de Computación** representan un fascinante punto de convergencia entre la **Teoría de Grafos** y la **Computación Funcional**, dos campos aparentemente dispares que encuentran un terreno de intersección en la modelización y resolución de problemas computacionales complejos.
A diferencia de los grafos matemáticos convencionales, donde los nodos y las aristas pueden representar entidades abstractas y relaciones entre ellas, los grafos de computación se centran en la ejecución de operaciones funcionales en los nodos, lo que permite representar y abordar una amplia gama de problemas algorítmicos de manera eficiente y elegante.
En términos sencillos,
!!!def:Grafo de Computación
Un grafo de computación consiste en un conjunto de nodos interconectados por aristas, donde cada nodo representa una función, o una operación computacional, y las aristas definen las dependencias entre estas operaciones, el flujo de datos que se trasladan entre nodos y determinan tanto la salida como la entrada de las operaciones que contienen.
Esta estructura de grafo permite modelar y visualizar de manera clara y concisa la secuencia de operaciones necesarias para realizar un cálculo o resolver un problema específico.
 Una de las características más interesantes de los grafos de computación es su capacidad para reflejar de forma natural la paralelización y la distribución de tareas en sistemas computacionales modernos, algo que difícilmente se puede conseguir con la representación símbolica clásica. Por ejemplo, al organizar las operaciones en nodos independientes y definir las dependencias entre ellos mediante aristas direccionales, los grafos de computación permiten representar la esencia de los recursos de hardware disponibles, como múltiples núcleos de procesamiento o incluso sistemas distribuidos.
Podemos encontrar aplicaciones de estas estructuras en una amplia variedad de áreas, desde el procesamiento de datos y el aprendizaje automático hasta la optimización de algoritmos y la simulación de sistemas complejos. En el contexto del aprendizaje automático, por ejemplo, los grafos de computación se utilizan para representar y entrenar redes neuronales profundas, donde cada nodo del grafo corresponde a una operación matemática, como una multiplicación de matrices o una función de activación, y las aristas definen el flujo de datos a través de la red.
Otro aspecto destacado de los grafos de computación es su compatibilidad con enfoques de programación funcional, como el lambda-cálculo y la programación basada en funciones de orden superior. Al modelar las operaciones computacionales como funciones puras que transforman datos de entrada en resultados de salida, los grafos de computación fomentan la modularidad, la reutilización y la composición de código, lo que facilita el desarrollo y mantenimiento de sistemas de alta complejidad.
Desde el punto de vista expresivo, los grafos de computación son una forma interesante de entender las expresiones matemáticas. Por ejemplo, consideremos la expresión $z= \sin(x)+xy$. Hay tres operaciones: una suma, una multiplicación y un seno. Para facilitar su representación, introduzcamos dos variables intermedias, $h_1$ y $h_2$, de forma que cada salida de la función tenga una variable. Ahora tenemos:
$$\begin{split}
h_1(x)&=\sin(x)\\
h_2(x,y)&=xy\\
z(h_1,h_2)&=h_1+h_2
\end{split}
$$
Para crear un grafo de computación, convertimos cada una de estas operaciones, junto con las variables de entrada, en nodos. Cuando el resultado de una operación sirve de entrada de otra, se coloca una flecha del nodo que representa la primera al nodo de la segunda:
Podemos evaluar la expresión ajustando las variables de entrada a ciertos valores y calculando los nodos a través del grafo. Por ejemplo, si establecemos $x=\pi$ y $y=1$, la expresión anterior se evalúa en $\pi$.
# Diferenciación Automática
La diferenciación es un proceso rutinario, normalmente basado en descomponer funciones más complejas en pequeños trozos más simples que sabemos diferenciar y, a partir de ahí, ensamblar las piezas obtenidas para reconstruir (calcular) el gradiente de la función original.
Por ello, la pieza esencial es la diferenciación de la función compuesta, que se realiza por la llamada **regla de la cadena**.
!!!def:Regla de la Cadena
Si tenemos $f(x)=(f_1\circ f_2)(x)$, entonces la regla de la cadena propociona la siguiente igualdad:
$$
\left.\frac{\partial f}{\partial x}\right|_{x_0} = \left.\frac{\partial f_1}{\partial y_1}\right|_{f_2(x_0)} \cdot \left.\frac{\partial f_2}{\partial y_2}\right|_{x_0}
$$
donde $y_1$ e $y_2$ son las entradas respectivas de $f_1$ y $f_2$.
En el caso de tener $f(x) = (f_1 \circ f_2 \circ \ldots \circ f_n)(x)$, la expresión anterior se generaliza de la siguiente forma:
$$
\left.\frac{\partial f}{\partial x}\right|_{x_0} = \left.\frac{\partial f_1}{\partial y_1}\right|_{y_1^0} \cdot \left.\frac{\partial f_2}{\partial y_2}\right|_{y_2^0} \cdot \ldots \cdot \left.\frac{\partial f_n}{\partial y_n}\right|_{y_n^0}
$$
donde $y_i$ denota la entrada de la función $f_i$ y
$$
\begin{alignat*}{2}
y_n^0 = &\ x_0 \\
y_i^0 = &\ \left(f_{i+1} \circ \ldots \circ f_n\right) (x_0) \\
\end{alignat*}
$$
La forma concreta de $\left.\frac{\partial f_i}{\partial y_i}\right|_{y_i^0}$ depende de la dimensión de los espacios de partida y llegada:
* Si $f_i: \mathbb{R} \rightarrow \mathbb{R}$, entonces $\frac{\partial f_i}{\partial y_i} \in \mathbb{R}$ y todo lo que se necesita es multiplicar números.
* Si $f_i: \mathbb{R}^{m_i} \rightarrow \mathbb{R}^{n_i}$, entonces $\left.\frac{\partial f_i}{\partial y_i}\right|_{y_i^0} =\mathbf{J}_i \in \mathbb{R}^{n_i\times m_i}$ es una matriz con $m_i$ filas y $n_i$ columnas.
Por tanto, el cálculo del gradiente $\frac{\partial f}{\partial x}$ *teóricamente* se consigue con el siguiente procedimiento:
1. Se calculan los Jacobianos $\left\{\mathbf{J}_i\right\}_{i=1}^n$,
2. Se multiplican los Jacobianos: $\left.\frac{\partial f}{\partial x}\right|_{y_0} = \mathbf{J}_1 \times \mathbf{J}_2 \times \ldots \times \mathbf{J}_n$.
!!!side:2
Tiene una complejidad de al menos $O(n^{2,3728596})$, que representa la cota mínima conocida para el producto de matrices. Ver [aquí](https://arxiv.org/abs/2010.05846).
!!!side:3
Para ciertas funciones puede ser recomendable un modo mixto, en el que se combinan ambos casos, donde algunas partes se calculan usando el modo directo, y otras el inverso.
Por tanto, la complejidad del cálculo (al menos, una parte de él) viene determinada por la multiplicación de matrices, que suele ser cara computacionalmente [2]. El orden en que se multiplican los jacobianos tiene, por tanto, un profundo efecto sobre la complejidad del motor AD. Aunque determinar el orden óptimo de multiplicación de la secuencia de matrices es costoso, en la práctica, reconocemos dos casos importantes [3]:
1. Los Jacobianos se multiplican de derecha a izquierda: $J_1 \times (J_2 \times ( \ldots \times (J_{n-1} \times J_n) \ldots))$, que tiene ventajas cuando la dimensión de entrada de $f: \mathbb{R}^k \rightarrow \mathbb{R}^m$ es menor que la dimensión de salida, $k < m$. Se denomina **Modo Directo**.
2. Los Jacobianos se multiplican de izquierda a derecha: $( \ldots ((J_1 \times J_2) \times J_3) \times \ldots ) \times J_n$, que tiene ventajas cuando la dimensión de entrada de $f: \mathbb{R}^k \rightarrow \mathbb{R}^m$ es mayor que la dimensión de salida, $k > m$. Se denomina **Modo Inverso**.
!!!Ejemplo
Consideremos el ejemplo:
$$ f(x,y) = xy + sin(x)$$
En este caso: $(f_1 \circ f_2 \circ f_3)(x,y) = xy + sin(x)$ donde:
$$
\begin{alignat*}{6}
f_1:&\mathbb{R}^2 \rightarrow \mathbb{R} \quad&f_1(y_1)& = y_{1,1} + y_{1,2} \quad & y_0 = & (xy + \sin(x)) \\
f_2:&\mathbb{R}^3 \rightarrow \mathbb{R}^2 \quad&f_2(y_2)& = (y_{2,1}y_{2,2}, y_{2,3}) \quad & y_1 = & (xy, \sin(x))&\\
f_3:& \mathbb{R}^2 \rightarrow \mathbb{R}^3 \quad&f_3(y_3)& = (y_{3,1}, y_{3,2}, \sin(y_{3,1})) \quad & y_2 =& (x, y, \sin(x))\\
\end{alignat*}
$$
Los correspondientes Jacobianos serían:
$$
\begin{alignat*}{4}
f_1(y_1) & = y_{1,1} + y_{1,2} \quad & \mathbf{J}_1& = \begin{bmatrix} 1 \\ 1 \end{bmatrix} \\
f_2(y_2) & = (y_{2,1}y_{2,2}, y_{2,3}) \quad & \mathbf{J}_2& = \begin{bmatrix} y_{2, 2} & 0 \\ y_{2,1} & 0 \\ 0 & 1 \end{bmatrix}\\
f_3(y_3) & = (y_{3,1}, y_{3,2}, \sin(y_{3,1})) \quad & \mathbf{J}_3 & = \begin{bmatrix} 1 & 0 & \cos(y_{3,1}) \\ 0 & 1 & 0 \end{bmatrix} \\
\end{alignat*}
$$
y para el gradiente tendríamos:
$$
\begin{bmatrix} \frac{\partial f(x, y)}{\partial{x}} \\ \frac{\partial f(x,y)}{\partial{y}} \end{bmatrix} = \mathbf{J}_3 \times \mathbf{J}_2 \times \mathbf{J}_1 = \\ = \begin{bmatrix} 1 & 0 & \cos(x) \\ 0 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} y & 0 \\ x & 0 \\ 0 & 1 \end{bmatrix} \times \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} y & \cos(x) \\ x & 0 \end{bmatrix} \times \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} y + \cos(x) \\ x \end{bmatrix}
$$
Nótese que, desde el punto de vista teórico, la expresión de una función como composición no es única, sino que normalmente viene dada por el grafo computacional en un lenguaje concreto.
En las siguientes secciones analizaremos con cierta profundidad (y con una implementación inicial en Julia) los dos modos posibles.
## Modo Directo
!!!side:4
El método que se propone es un cálculo idealizado. La aplicación real es un poco diferente, como veremos más adelante.
En teoría, podemos calcular el gradiente utilizando el modo directo de la siguiente forma [4]:
1. Inicializar el Jacobiano de $y_n$ con respecto a $x$ a una matriz identidad: como $y^0_n = x$, se tiene que $\frac{\partial y_n}{\partial x} = \mathbb{I}$.
2. Iterar $i$ desde $n$ hasta 1 hacer:
* Calcular la siguiente salida intermedia como: $y^0_{i-1} = f_i({y^0_i})$.
* Calcular el Jacobiano: ${\mathbf J}_i = \left.\frac{f_i}{\partial y_i}\right|_{y^0_i}$
* Adelantar el gradiente por medio de : $\left.\frac{\partial y_{i-1}}{\partial x}\right|_x = {\mathbf J}_i \times \left.\frac{\partial y_i}{\partial x}\right|_x$.
Obsérvese que:
* Al final, nos quedamos con $y = y^0_0$ y $\frac{\partial y_0}{\partial x}$, que es el gradiente que queremos calcular.
* Si $y$ es un escalar, entonces $\frac{\partial y_0}{\partial x}$ es una matriz formada por una fila simple.
* El Jacobiano y la salida de la función se calculan en una sola pasada.
### Números Duales
El modo directo necesita llevar la cuenta de la salida de la función y de la derivada en cada paso del cálculo de la función $f$. Esto se puede realizar elegantemente con la introducción de los [**números duales**](https://en.wikipedia.org/wiki/Dual_number), que son conceptualmente similares a los números complejos, pero en lugar del número imaginario $i$ utilizan $\varepsilon$ en su segunda componente:
$$x = v + \dot v \varepsilon$$
!!!side:4
En vez de $i^2=-1$ como en los números complejos. La colección de números duales forma un álgebra asociativa conmutativa y unitaria bidimensional particular sobre los números reales. Fueron introducidos en 1873 por William Clifford, y utilizados a principios del siglo XX por el matemático alemán Eduard Study para estudiar la posición relativa de dos rectas en el espacio.
donde $(v,\dot v) \in \mathbb R^2$ y por definición $\varepsilon^2=0$ [4].
!!!teorema:Propiedades de los Números Duales
Comencemos calculando cómo se comportan con respecto a las operaciones aritméticas básicas:
$$
\begin{array}{rll}
(v + \dot v \varepsilon) + (u + \dot u \varepsilon) &= (v + u) + (\dot v + \dot u)\varepsilon \\
(v + \dot v \varepsilon)(u + \dot u \varepsilon) &= vu + (u\dot v + \dot u v)\varepsilon + \dot v \dot u \varepsilon^2 &= vu + (u\dot v + \dot u v)\varepsilon \\
{\displaystyle \frac{v + \dot v \varepsilon}{u + \dot u \varepsilon}} &= {\displaystyle \frac{v + \dot v \varepsilon}{u + \dot u \varepsilon} \frac{u - \dot u \varepsilon}{u - \dot u \varepsilon}} &= {\displaystyle \frac{v}{u} - \frac{(\dot u v - u \dot v)}{u^2}\varepsilon}
\end{array}
$$
Es común, al igual que con los números complejos, pasar de la notación anterior a la notación de pares, por lo que tendríamos:
$$
\begin{array}{ll}
(v, \dot v) + (u, \dot u) &= (v + u,\ \dot v + \dot u)\\
(v, \dot v)(u, \dot u ) &= (vu,\ u\dot v + \dot u v)\\
{\displaystyle \frac{(v, \dot v)}{(u,\dot u)}} &= {\displaystyle \Big(\frac{v}{u},\ - \frac{\dot u v - u \dot v}{u^2}\Big)}
\end{array}
$$
Si interpretamos un número dual como si fuera la representación de una parte funcional y su derivada (en $(u+\dot u \varepsilon)$, $u$ sería la función, y $\dot u$ sería su derivada), podemos ver que la *dualidad* representa las propiedades de derivación básicas.
En particular, si evaluamos la ecuaciones anteriores en $(v, \dot v) = (v, 1)$ y $(u, \dot u) = (u, 0)$:
$$
\begin{array}{ll}
(v, 1) + (u,0) &= (v + u,\ 1) \\
(v, 1)(u,0) &= (vu,\ u)\\
{\displaystyle \frac{(v, 1)}{(u,0)}} &= {\displaystyle \Big(\frac{v}{u},\ \frac{1}{u}\Big)}
\end{array}
$$
los términos de la parte dual: $(1, u, \frac{1}{u})$, se corresponden con el gradiente de las funciones originales $(u+v, uv, \frac{v}{u})$ con respecto a $v$.
Si repetimos los cálculos con $(v, \dot v) = (v, 0)$ y $(u, \dot u) = (u, 1)$:
$$
\begin{array}{ll}
(v,0)+(u,1) &= (v + u,\ 1)\\
(v ,0)(u,1) &= (vu,\ v)\\
{\displaystyle \frac{(v,0)}{(u,1)} } &= {\displaystyle (\frac{v}{u}, - \frac{v}{u^2})}
\end{array}
$$
obtenemos en las partes duales los gradientes con respecto a $u$.
El siguiente resultado muestra que los números duales pueden usarse para calcular el gradiente de polinomios:
!!!teorema
Si $p(v) = \sum_{i=1}^n p_iv^i$, entonces
$$p(v + \dot{v} \varepsilon)= p(v) + p'(v) \dot v\varepsilon$$
!!!alg:Demostración
Dado el polinomio: $p(v) = \sum_{i=1}^n p_iv^i$:
$$
\begin{split}
p(v + \dot{v} \varepsilon) &=
\sum_{i=0}^n p_i(v + \dot{v} \varepsilon )^i =
\sum_{i=0}^n \left[p_i \sum_{j=0}^{n}\binom{i}{j}v^{i-j}(\dot v \varepsilon)^{i}\right] =
p_0 + \sum_{i=1}^n \left[p_i \sum_{j=0}^{1}\binom{i}{j}v^{i-j}(\dot v \varepsilon)^{j}\right] = \\
&= p_0 + \sum_{i=1}^n p_i(v^i + i v^{i-1} \dot v \varepsilon )
= p(v) + \left(\sum_{i=1}^n ip_i v^{i-1}\right) \dot v \varepsilon
\end{split}
$$
donde el término dual: $\sum_{i=1}^n ip_i v^{i - 1}$, se corresponde con la derivada de $p(v)$ con respecto a $v$.
Y, gracias a la expansión de Taylor, podemos generalizar este resultado a cualquier función suficientemente suave:
!!!teorema
Si tenemos una función genérica $f:\mathbb{R} \rightarrow \mathbb{R}$. Su valor en el punto $v + \dot v \varepsilon$ puede ser aproximado usando la expansión de Taylor en el punto $v$ como:
$$
f(v+\dot v \varepsilon) = \sum_{i=0}^\infty \frac{f^i(v)\dot v^i\varepsilon^n}{i!}
= f(v) + f'(v)\dot v\varepsilon
$$
Ya que todos los términos de orden superior se eliminan porque $\varepsilon^i=0$ para $i>1$.
### Implementación en Julia
Para demostrar la simplicidad de los números duales, considera la siguiente definición que hace uso de un nuevo tipo (derivado de `Number`) y sobrecarga las funciones aritméticas en Julia:
```julia
struct Dual{T<:Number} <: Number
x::T
d::T
end
Base.:+(a::Dual, b::Dual) = Dual(a.x+b.x, a.d+b.d)
Base.:-(a::Dual, b::Dual) = Dual(a.x-b.x, a.d-b.d)
Base.:/(a::Dual, b::Dual) = Dual(a.x/b.x, (a.d*b.x - a.x*b.d)/b.x^2)
Base.:*(a::Dual, b::Dual) = Dual(a.x*b.x, a.d*b.x + a.x*b.d)
# Reglas de Promoción
Dual(x::S, d::T)
where {S<:Number, T<:Number} = Dual{promote_type(S, T)}(x, d)
Dual(x::Number) = Dual(x, zero(typeof(x)))
Dual{T}(x::Number)
where {T} = Dual(T(x), zero(T))
Base.promote_rule(::Type{Dual{T}}, ::Type{S})
where {T<:Number,S<:Number} = Dual{promote_type(T,S)}
Base.promote_rule(::Type{Dual{T}}, ::Type{Dual{S}})
where {T<:Number,S<:Number} = Dual{promote_type(T,S)}
# Añadimos una API para la diferenciación directa
forward_diff(f::Function, x::Real) = _dual(f(Dual(x,1.0)))
_dual(x::Dual) = x.d
_dual(x::Vector) = _dual.(x)
```
Un ejemplo de uso de las funciones definidas podría ser:
```julia
julia> forward_diff(x -> [1 + x, 5x, 5/x], 2)
3-element Vector{Float64}:
1.0
5.0
-1.25
```
Pero lo realmente interesante de esta implementación no es que calcule las derivadas de expresiones matemáticas, sino que también puede calcular funciones que vienen dadas como algoritmos numéricos:
!!!ejemplo:Raíz cuadrada Babilónica
Vamos a comprobar cómo funcionan los números duales sobre un algoritmo que aproxima $\sqrt x$, de la que sabemos que:
$$f(x) = \sqrt{x}, \qquad f'(x) = \frac{1}{2\sqrt{x}}$$
En particular, vamos a usar el algoritmo babilónico, que se centra en el hecho de que cada lado de un cuadrado es la raíz cuadrada del área, y fue usado durante muchos años para calcular raíces cuadradas a mano debido a su gran eficacia y rapidez.
Si definimos la función auxiliar (que se define por recurrencia):
$$h(x,t,n) = \begin{cases}
t & \text{, si } n=0\\
h(x, \frac{1}{2}(t+\frac{x}{t}), n-1)) & \text{, si } n>0\\
\end{cases}$$
entonces la raíz cuadrada se puede aproximar por ($n=10$ determina una suficiente cantidad de iteraciones para obtener un resultado suficientemente aproximado):
$$\sqrt{x}\ ≈ \ h(x,\frac{1+x}{2},10)$$
En Julia, podríamos escribirlo como:
```julia
julia> babsqrt(x, t=(1+x)/2, n=10) = n==0 ? t : babysqrt(x, (t+x/t)/2, n-1)
babsqrt (generic function with 3 methods)
```
Y podemos comprobar los valores de aproximación que se consiguen, tanto de la función, como de su derivada en un punto concreto:
```julia
julia> babsqrt(2) ≈ sqrt(2)
true
julia> forward_diff(babsqrt, 2)
0.35355339059327373
julia> forward_diff(babsqrt, 2) ≈ 1/(2sqrt(2))
true
```
Vamos a comparar la solución que devuelve `forward_diff` con la que devuelve un método estándar como es el de diferencias finitas y la solución analítica exacta:
```julia
julia> using FiniteDifferences
julia> forward_dsqrt(x) = forward_diff(babsqrt,x)
forward_dsqrt (generic function with 1 method)
julia> analytc_dsqrt(x) = 1/(2babsqrt(x))
analytc_dsqrt (generic function with 1 method)
julia> forward_dsqrt(2.0)
0.35355339059327373
julia> analytc_dsqrt(2.0)
0.3535533905932738
julia> central_fdm(5, 1)(babsqrt, 2.0)
0.3535533905932651
```
Y podemos dar una representación gráfica para ver hasta qué punto la aproximación de la derivada obtenida por AD en modo directo y por medios analíticos coinciden:
```julia
plot(0.0:0.01:2, babsqrt, label="babsqrt(x)", lw=3)
plot!(0.1:0.01:2, analytc_dsqrt, label="Analytic f'", ls=:dot, lw=3)
plot!(0.1:0.01:2, forward_dsqrt, label="Dual Forward Mode f'", lw=3, ls=:dash)
```
!!!nota
Hay algunas notas interesantes que podemos extraer del ejemplo anterior:
1. La función $f'$ que devuelve el modo directo se obtiene simplemente pasándole un dual a la función original.
2. Para implementar el método AD en modo directo en Julia sólo necesitamos _sobrecargar_ unos pocos _operadores_.
3. Se puede extender adecuadamente a funciones vectoriales trabajando con una idea similar que se llama [**Números Hiperduales**](http://adl.stanford.edu/hyperdual/).
5. Aunque hay versiones más modernas, la que hemos visto aquí es la que se implementa en la librería [`ForwardDiff.jl`](https://github.com/JuliaDiff/ForwardDiff.jl), y es bastante estable y general.
6. En la actualidad, las implementaciones directas más modernas hacen uso del Jacobiano, y no de los duales.
## Modo Inverso
!!!side:5
La AD en modo inverso fue publicada por primera vez en 1976 por Seppo Linnainmaa, un informático finlandés. Se popularizó a finales de los 80 cuando se aplicó al entrenamiento de perceptrones multicapa, lo que dio lugar al famoso algoritmo de **backpropagation**, que es un caso especial de AD en modo inverso.
* Linnainmaa, S. (1976). Taylor expansion of the accumulated rounding error. *BIT Numerical Mathematics*, 16(2), 146-160.
* Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986). Learning representations by back-propagating errors. *Nature*, 323, 533--536.
En el modo inverso la computación del gradiente sigue el orden contrario: comenzamos la computación haciendo $\mathbf{J}_0 = \frac{\partial y}{\partial y_0},$ que es de nuevo la matriz identidad. Entonces, calculamos los Jacobianos y las multiplicaciones en el orden inverso. El problema es que para calcular $J_i$ necesitamos conocer el valor de $y_i^0$, que no puede ser calculado en el pase inverso, por lo que hemos de hacer primero un pase directo, donde se calculan todos los $\{y_i^0\}_{i=1}^n$ [5].
Por tanto, el algoritmo completo del modo inverso funciona como sigue:
1. Paso directo: itera $i$ desde 1 hasta $n$ y hacer:
* Calcular la siguiente salida intermedia como $y^0_{i-1} = f_i(y^0_i)$.
2. Paso inverso: itera $i$ desde $n$ hasta 1 y hacer:
* Calcular el Jacobiano: $J_i = \left.\frac{f_i}{\partial y_i}\right|_{y_i^0}$ at point $y_i^0$.
* *Empujar* hacia atrás el gradiente por medio de: $\left.\frac{\partial f(x)}{\partial y_{i}}\right|_{y^0_i} = \left.\frac{\partial y_0}{\partial y_{i-1}}\right|_{y^0_{i-1}} \times J_i$.
La necesidad de almacenar salidas intermedias tiene un gran impacto en los requisitos de memoria, que particularmente en la GPU es un gran problema.
!!!alg:Trucos para reducir la memoria
Se han propuesto (y usado) varios métodos que pueden servir para no hacer un uso tan excesivo de la memoria y así facilitar su implementación en sistemas con menos recursos (sobre todo, ahora que el tamaño de las redes es tan grande). Entre las que podemos destacar la definición de **reglas personalizadas** que *compriman* grandes bloques funcionales; cachear salidas intermedias de las funciones que sean invertibles; o no almacenar los resultados intermedios después de una secuencia larga de operaciones, volviendo a calcularlas si fuera necesario en otros pasos adelante (es lo que se denomina **Checkpointing**).
La mayoría de los motores AD de modo inverso no soportan la mutación de valores de matrices, porque habitualmente después de cada mutación se necesitaría guardar la matriz completa (con el consiguiente gasto de tiempo y memoria). Existe un paquete, [`Enzyme`](https://github.com/wsmoses/Enzyme.jl), que hace AD directamente sobre el código LLVM, ya que en LLVM cada variable se asigna una sola vez. Los métodos directos no sufren este problema, ya que el gradiente se calcula en el momento de la evaluación de la función.
El AD en modo inverso necesita registrar las operaciones sobre las variables cuando calcula el valor de una función derivada, de forma que pueda volver atrás cuando calcule el gradiente. Aunque a este registro se le denomina **tape**, en realidad es un grafo de computación como los que hemos presentado anteriormente. La construcción de este grafo puede ser *explícita* o *implícita*, y el código que calcula el gradiente puede obtenerse mediante técnicas de *sobrecarga de operadores* o de *reescritura de código*. Esta combinatoria da lugar a cuatro versiones diferentes de AD, y Julia tiene bibliotecas para las cuatro:
| Librería | Construcción Grafo | Cálculo Gradiente |
|----------|:-------------------:|:------------------:|
|[`Yota.jl`](https://github.com/dfdx/Yota.jl) |explícita|reescritura de código|
|[`ReverseDiff.jl`](https://github.com/JuliaDiff/ReverseDiff.jl)|explícita|sobrecarga de operadores|
|[`Zygote.jl`](https://github.com/FluxML/Zygote.jl)|implícita|reescritura de código|
|[`Tracker.jl`](https://github.com/FluxML/Tracker.jl), [`AutoGrad.jl`](https://github.com/denizyuret/AutoGrad.jl)|implícita|sobrecarga de operadores|
A continuación vamos a dar algunos detalles de los distintas aproximaciones en la construcción del grafo.
### AD basada en grafos
En el enfoque basado en grafos partimos de un conocimiento completo del grafo de computación (que se conoce en muchos casos, como las redes neuronales clásicas) y lo aumentamos con nodos que representan la computación del cálculo del gradiente (camino hacia atrás). Debemos tener cuidado de añadir todas las aristas que representan el flujo de información necesario para calcular el gradiente. Una vez aumentado el grafo de computación, es fácil determinar el subgrafo necesario para calcular los gradientes del nodo (o nodos) deseados.
Vamos a mostrar el proceso sobre un caso concreto. La siguiente figura muestra el grafo de computación asociado a la función $f(x, y) = \sin(x) + xy$ que presentamos al principio de esta entrada:
En el grafo, las flechas $\rightarrow$ denotan el flujo de operaciones, se denota la salida de la función $f$ por $z$, y las salidas de los nodos intermedios por $h_i$ (que significa *oculto*).
Comenzamos en la parte superior y añadimos un nodo que calcula $\frac{\partial z}{\partial h_3}$, que es la identidad, necesaria para disparar el inicio de la diferenciación.
Lo conectamos con la salida de $h_3$, aunque técnicamente en este caso no es necesario, ya que $z = h_3$.
A continuación, calculamos los padres de $h_3$, que son $h_1$ y $h_2$.
Comenzando por $h_2$, añadimos un nodo que calcula $\frac{\partial h_3}{\partial h_2}$ para el que sólo necesitamos información sobre $h_2$ y lo marcamos en el grafo (de nuevo, esta arista puede ser teóricamente eliminada debido a que es igual a 1 independientemente de las entradas). Siguiendo la regla de la cadena, tenemos que combinar $\frac{\partial h_3}{\partial h_2}$ con $\frac{\partial z}{\partial h_3}$ para calcular $\frac{\partial z}{\partial h_2}$ que también anotamos en el grafo.
Realizamos el mismo proceso con $h_1$, y calculamos $\frac{\partial h_3}{\partial h_1}$, que volvemos a combinar con $\frac{\partial z}{\partial h_3}$ para obtener $\frac{\partial z}{\partial h_1}$.
Continuando con este proceso de forma iterada, subiendo por los padres de cada nodo, obtenemos el grafo final siguiente:
que contiene los nodos buscados: $\frac{\partial z}{\partial x}$ y $\frac{\partial z}{\partial y}$. Este grafo computacional puede pasarse al compilador para que calcule los valores deseados cada vez que sea necesario.
Para aquellos que sean más cercanos a Python, diremos que este enfoque ha sido adoptado, por ejemplo, por [Theano](https://github.com/Theano/Theano) y por [TensorFlow](https://www.tensorflow.org/). En Tensorflow cuando usas funciones como `tf.mul(a,b)` o `tf.add(a,b)`, no estás realizando el cálculo en Python, sino que estás construyendo el grafo computacional anterior. A continuación, se pueden calcular los valores utilizando `tf.run` con las entradas deseadas, pero en realidad se calculan los valores en un intérprete/compilador diferente al de Python.
Esta aproximación presenta sus ventajas e inconvenientes:
!!!def:Ventajas
* Conocer el grafo de computación de antemano es genial, ya que se pueden realizar pasos de optimización para simplificar el grafo y reducir costes de computación posteriores.
* El grafo de computación tiene una semántica simple (soporte limitado para bucles, ramas, sin objetos), y el compilador de grafos es, por tanto, más simple que el compilador de lenguajes completos.
* Como el cálculo del gradiente aumenta el grafo, se puede volver a ejecutar el proceso para obtener gradientes de orden superior a costa de añadir unos pocos nodos adicionales.
* TensorFlow permite prefijar los tamaños de los tensores involucrados, lo que significa que sabe con precisión cuánta memoria se necesitará y dónde, lo que disminuye el número de asignaciones (algo de especial importancia en un dispositivo tan rígido como una GPU).
!!!def:Inconvenientes
* Está restringido a un grafo de computación fijo. Generalmente es difícil implementar `if` o `while` y, por lo tanto, cambiar la computación de acuerdo a los valores calculados durante el paso hacia adelante.
* El desarrollo y la depuración pueden ser difíciles, ya que no se está desarrollando el grafo de computación en el lenguaje anfitrión en el que está definido el proceso de optimización que hace uso de los gradientes calculados.
* Explotar el paralelismo dentro del grafo de computación puede ser difícil.
### AD basada en tracking
Una alternativa a los métodos basados en grafos estáticos, como el explicado en la sección anterior, son los métodos que construyen el grafo durante la llamada a las funciones y luego utilizan este grafo construido dinámicamente para saber cómo calcular el gradiente. Como hemos dicho anteriormente (y que repetimos por el uso extendido de este término en esta aproximación), el grafo construido dinámicamente también se suele denominar *tape*. Este enfoque es utilizado por bibliotecas muy populares, como [**_PyTorch_**](https://pytorch.org/), [**_AutoGrad_**](https://github.com/HIPS/autograd), y [**_Chainer_**](https://chainer.org/) en el ecosistema Python, o por [**_Tracker.jl_**](https://github. com/FluxML/Tracker.jl) (antiguo backend AD de `Flux.jl`), [**_ReverseDiff.jl_**](https://github.com/JuliaDiff/ReverseDiff.jl), y [**_AutoGrad.jl_**](https://github.com/denizyuret/AutoGrad.jl) (backend AD de `Knet.jl`) en Julia. Este tipo de sistemas AD también suelen hacer uso de la *sobrecarga de operadores*, ya que para registrar las operaciones realizadas sobre los argumentos necesitan reemplazar/envolver la implementación original.
Esta aproximación también presenta ventajas e inconvenientes:
!!!def:Ventajas
* La depuración y el desarrollo son más sencillos, ya que la AD está implementada en el mismo lenguaje.
* El grafo de computación es dinámico, lo que hace más sencillo calcular el gradiente en presencia del uso de instrucciones como `if` y `while`.
!!!def:Inconvenientes
* El grafo de computación se crea y modifica durante cada computación, lo que puede resultar costoso. En la mayoría de las aplicaciones de aprendizaje profundo, esta sobrecarga es insignificante en comparación con el tiempo necesario para realizar el resto de operaciones involucradas.
* Como el grafo de computación es dinámico, no se puede optimizar como el grafo estático, lo mismo ocurre con las asignaciones de memoria.
Un ejemplo más completo que permite entrenar una red neuronal feed-forward en la GPU se puede encontrar [aquí](ffnn.jl).
!!!alg
La diferencia entre los sistemas AD basados en tracking y los basados en grafos es conceptualmente similar a los lenguajes de programación interpretados y compilados. Los sistemas AD de tracking interpretan la función original y el grafo mientras calculan el gradiente, mientras que los sistemas AD basados en grafos compilan el cálculo del gradiente.
### Source-to-source AD
Hay un tercer mecanismo que es especialmente prometedor en lenguajes como Julia (y de difícil aplicación en lenguajes como Python). La idea se basa en que este tipo de lenguajes generan un código intermedio de compilación que se asemeja mucho al grafo de computación requerido (de hecho, es equivalente o muy cercano por medio de transformaciones sencillas). Por ejemplo, la función:
```julia
f(x,y) = x*y + sin(x)
```
se compila en:
```julia
julia> @code_lowered f(1.0, 1.0)
CodeInfo(
1 ─ %1 = x * y
│ %2 = Main.sin(x)
│ %3 = %1 + %2
└── return %3
)
```
Esta forma es particularmente adecuada para la diferenciación automática, ya que tenemos en el lado izquierdo de las asignaciones siempre una sola variable, lo que significa que el compilador nos ha proporcionado una forma en la que sabemos cómo aplicar las reglas AD.
## ChainRules
De lo que hemos hablado hasta ahora sobre los sistemas AD se desprende que, aunque la parte básica, el *motor*, es relativamente sencilla, lo difícil es escribir las reglas que indican cómo realizar el cálculo de los gradientes. Estas reglas son necesarias para todos los sistemas, ya estén basados en grafos, o en trazas.
!!!side:6
Iniciada por Catherine Frames White y respaldada por [Invenia](https://github.com/invenia).
!!!side:7
Hasta ahora, están soportadas por `Zygote.jl`, `Nabla.jl`, `ReverseDiff.jl` y `Diffractor.jl`, lo que sugiere que el enfoque de unificación está funcionando (pero no por `Enzyme.jl`).
Obviamente, es un esfuerzo adicional (e indeseable) que cada sistema AD tenga su propio conjunto de reglas. Por ello, la comunidad [6] ha empezado a trabajar en un sistema unificado para expresar las reglas de diferenciación, de forma que puedan compartirse entre sistemas [7]: **ChainRules**.
# Ejemplo completo
!!!note-right
Esta sección está extraída de [Rufflewind](https://rufflewind.com/2016-12-30/reverse-mode-automatic-differentiation) (con adaptaciones a Julia).
Supongamos que queremos encontrar la derivada $\frac{df}{dx}$ de la función $f:\mathbb R^2 \rightarrow \mathbb R$:
```julia
f(x,y) = x*y + sin(x)
```
Si tenemos reglas para `*`, `+`, y `sin` podríamos simplemente *alimentar* la función con `Dual(x, one(x))` y leer la derivada $\frac{df}{dx}$ del `Dual` que devuelve `f`. Si también estamos interesados en la derivada $\frac{df}{dy}$ tendremos que ejecutar `f` de nuevo, esta vez alimentando el segundo argumento con `Dual(y, one(y))`.
```julia
dfdx = f(Dual(x,one(x)), Dual(y,zero(y)))
dfdy = f(Dual(x,zero(x)), Dual(y,one(y)))
```
Por lo tanto, tenemos que evaluar `f` *dos veces* si queremos derivadas con respecto a sus dos argumentos, lo que significa que la diferenciación hacia adelante escala como $O(n)$ donde $n$ es el número de entradas a `f`.
El modo inverso de AD puede calcular gradientes de funciones con muchas entradas y una salida de una sola vez. Esto es genial porque muy a menudo queremos optimizar las funciones de pérdida que son exactamente eso: *funciones con muchas variables de entrada y una salida* (el valor numérico de la pérdida).
Con funciones $f:\mathbb R^n\rightarrow\mathbb R^m$ y $g:\mathbb R^k\rightarrow \mathbb R^n$ con un vector de entrada ${\mathbf x}$ podemos definir la composición de $f$ y $g$ como
$${\mathbf z} = (f \circ g)({\mathbf x}), \qquad \text{donde} \qquad {\mathbf y}=g({\mathbf x}), \qquad {\mathbf z} = f({\mathbf y})$$
!!!side:7
Puedes consultar [aquí](https://math.stackexchange.com/questions/3785018/intuitive-proof-of-the-multivariable-chain-rule) o [aquí](https://people.math.harvard.edu/~shlomo/docs/Advanced_Calculus.pdf).
La regla de la cadena multivariable nos diría que [7]:
$$
\left.\frac{\partial z_i}{\partial x_j}\right|_{{\mathbf x}} =
\sum_{k=1}^n \left.\frac{\partial z_i}{\partial y_k}\right|_{{\mathbf y}}
\left.\frac{\partial y_k}{\partial x_i}\right|_{{\mathbf x}}
$$
que es, esencialmente, una fila de la *matriz jacobiana* $J$.
Debe tenerse en cuenta que, para calcular la derivada, siempre tenemos que conocer la entrada a la función respectiva, porque sólo podemos calcular la derivada *en un punto específico* (denotado por la notación $|_x$ $_{}$). En nuestro ejemplo, si consideramos $g(x,y)=xy$ y $h(x)=\sin(x)$ obtenemos
$$
\left.{\frac {df}{dx}}\right|_{x,y}
= \left.{\frac {df}{dg}}\right|_{g(x,y)}\cdot \left.{\frac {dg}{dx}}\right|_{x,y}
+ \left.{\frac {df}{dh}}\right|_{h(x)}\cdot \left.{\frac {dh}{dx}}\right|_{x}
= 1 \cdot y |_{y} + 1\cdot\cos(x)|_{x}
$$
Se puede ver que, con el fin de implementar el modo inverso AD tenemos que rastrear y recordar todas las entradas de nuestras funciones intermedias durante el paso hacia adelante de tal manera que podamos calcular sus gradientes durante el paso hacia atrás. La forma más sencilla de hacerlo es construyendo dinámicamente un grafo de computación que rastree cómo afecta cada variable de entrada a sus variables de salida. El siguiente grafo representa el cálculo de la función `f` considerada en esta entrada:
```julia
z = x*y + sin(x)
# Nodos del Grafo # Derivadas parciales
h1 = x*y # dh1/dx = y; dh1/dy = x
h2 = sin(x) # dh2/dx = cos(x)
z = h1 + h2 # dz/dh1 = 1; dz/dh2 = 1
```
En el grafo se puede ver que la variable `x` puede afectar directamente a `b` y `a`. Por lo tanto, `x` tiene dos hijos `a` y `b`. Durante el paso hacia adelante construimos el grafo, haciendo un seguimiento de qué entrada afecta a qué salida. Además, incluimos las derivadas locales correspondientes (que ya podemos calcular).
Para implementar un grafo construido dinámicamente podemos introducir un nuevo tipo de número, `TrackedReal`, que tiene tres campos:
* `data`: contiene el valor de este nodo en el grafo de computación, tal y como se obtuvo en el paso hacia delante.
* `grad`: se inicializa a `nothing` y más tarde contendrá los gradientes acumulados (la suma en la regla de la cadena multivariante)
* `children`: es un `Dict` que mantiene la traza de las variables de salida que se ven afectadas por el nodo actual, y también almacena las correspondientes derivadas locales $\frac{\partial f}{\partial g_k}$.
```julia
mutable struct TrackedReal{T<:Real}
data::T
grad::Union{Nothing,T}
children::Dict
# Este campo solo se necesita para imprimir el grafo
name::String
end
track(x::Real,name="") = TrackedReal(x,nothing,Dict(),name)
function Base.show(io::IO, x::TrackedReal)
t = isempty(x.name) ? "(tracked)" : "(tracked $(x.name))"
print(io, "$(x.data) $t")
end
```
El paso hacia atrás no es más que la aplicación de la regla de la cadena para calcular la derivada. Suponiendo que sabemos cómo calcular las *derivadas locales*, $\frac{\partial f}{\partial g_k}$, para funciones simples como `+`, `*`, y `sin`, podemos escribir una función simple que implementa la acumulación del gradiente de arriba a través de la regla de la cadena:
$$
\left.\frac{\partial f}{\partial x_i}\right|_{{\mathbf x}} =
\sum_{k=1}^n \left.\frac{\partial f}{\partial g_k}\right|_{\mathbf g({\mathbf x})}
\left.\frac{\partial g_k}{\partial x_i}\right|_{{\mathbf x}}
$$
Solo tenemos que hacer un bucle sobre todos los hijos, recoger las derivadas locales, y utilizar la recursión:
```julia
function accum!(x::TrackedReal)
if isnothing(x.grad)
x.grad = sum(w*accum!(v) for (v,w) in x.children)
end
x.grad
end
```
donde `w` corresponde a $\frac{\partial f}{\partial g_k}$ y `accum!(v)` corresponde a $\frac{\partial g_k}{\partial x_i}$. Llegados a este punto, ya hemos implementado la funcionalidad principal de nuestro primer AD de modo inverso. Lo único que queda por hacer es implementar las reglas inversas para las funciones básicas. A través de la recursividad de la regla de la cadena se aplica hasta que llegamos a la salida final, `z`. Esta salida final tiene que ser alimentada (al igual que con el modo hacia adelante) con $\frac{\partial z}{\partial z}=1$.
## Implementación de Reglas Inversas
Empecemos sobrecargando las tres funciones `+`, `*`, y `sin` que necesitamos para construir nuestro grafo de computación. En primer lugar, tenemos que rastrear el cálculo hacia adelante y luego *registrar* la salida `z` como un hijo de sus entradas mediante el uso de `z` como una clave en el diccionario de los hijos. El valor correspondiente contiene las derivadas, en el caso de la multiplicación simplemente tenemos:
$$z = a \cdot b$$
cuyas derivadas son:
$$
\frac{\partial z}{\partial a}=b, \qquad
\frac{\partial z}{\partial b}=a.
$$
Conociendo las derivadas de `*` en un punto dado podemos escribir nuestra regla inversa:
```julia
function Base.:*(a::TrackedReal, b::TrackedReal)
z = track(a.data * b.data, "*")
a.children[z] = b.data # dz/da=b
b.children[z] = a.data # dz/db=a
z
end
```
Veamos cómo se comporta lo que acabamos de definir. Crear dos números de rastreo y sumarlos da como resultado:
```julia
julia> x = track(2.0)
2.0 (tracked)
julia> y = track(3.0)
3.0 (tracked)
julia> z = x*y
6.0 (tracked *)
julia> x.children
Dict{Any, Any} with 1 entry:
6.0 (tracked *) => 3.0
julia> y.children
Dict{Any, Any} with 1 entry:
6.0 (tracked *) => 2.0
```
La implementación para las operaciones `+` y `sin` quedaría:
```julia
function Base.:+(a::TrackedReal{T}, b::TrackedReal{T}) where T
z = track(a.data + b.data, "+")
a.children[z] = one(T)
b.children[z] = one(T)
z
end
function Base.sin(x::TrackedReal)
z = track(sin(x.data), "sin")
x.children[z] = cos(x.data)
z
end
```
## Paso adelante y atrás
!!!side:8
Para poder realizar esta tarea, trabajaremos con los árboles abstractos que proporciona Julia. En cualquier caso, nos interesa la observación sobre la comparativa, no el mecanismo que realizamos para poder hacerla:
```julia
using AbstractTrees
AbstractTrees.children(v::TrackedReal) = v.children |> keys |> collect
function AbstractTrees.printnode(io::IO,v::TrackedReal)
print(io,"$(v.name) data: $(round(v.data,digits=2)) grad: $(v.grad)")
end
```
Para comprobar que con el modo inverso realmente ahorramos computación podemos representar el grafo de computación en diferentes etapas. Empezamos con el paso hacia delante y observamos `x` [8]:
```julia
julia> x = track(2.0,"x");
julia> y = track(3.0,"y");
julia> a = x*y;
julia> print_tree(x)
x data: 2.0 grad: nothing
└─ * data: 6.0 grad: nothing
```
Podemos ver que `x` ahora tiene un hijo `a` que tiene el valor `2.0*3.0==6.0`. Todos los gradientes siguen siendo `nothing`. Si calculamos otro valor que dependa de `x` añadiremos otro hijo:
```julia
julia> b = sin(x)
0.9092974268256817 (tracked sin)
julia> print_tree(x)
x data: 2.0 grad: nothing
├─ sin data: 0.91 grad: nothing
└─ * data: 6.0 grad: nothing
```
En el último paso calculamos `z`, que no muta los hijos de `x` porque no depende directamente de él. El resultado `z` se añade como hijo de `a` y `b`:
```julia
julia> z = a + b
6.909297426825682 (tracked +)
julia> print_tree(x)
x data: 2.0 grad: nothing
├─ sin data: 0.91 grad: nothing
│ └─ + data: 6.91 grad: nothing
└─ * data: 6.0 grad: nothing
└─ + data: 6.91 grad: nothing
```
Para el paso hacia atrás tenemos que alimentar el valor inicial del gradiente de `z` y llamar a `accum!` en la variable que nos interesa:
```julia
julia> z.grad = 1.0
1.0
julia> dx = accum!(x)
2.5838531634528574
julia> dx ≈ y.data + cos(x.data)
true
```
Al acumular los gradientes para `x`, se evaluarán los gradientes en el subárbol conectado a `x`. Las partes del árbol que sólo están conectadas a `y` permanecen intactas:
```julia
julia> print_tree(x)
x data: 2.0 grad: 2.5838531634528574
├─ sin data: 0.91 grad: 1.0
│ └─ + data: 6.91 grad: 1.0
└─ * data: 6.0 grad: 1.0
└─ + data: 6.91 grad: 1.0
julia> print_tree(y)
y data: 3.0 grad: nothing
└─ * data: 6.0 grad: 1.0
└─ + data: 6.91 grad: 1.0
```
Si ahora acumulamos los gradientes sobre `y` reutilizamos los gradientes ya calculados. En cálculos más grandes, esto ahorrará *mucha* computación.
## Optimización de funciones 2D
Vamos a comenzar implementando una función `gradient(f, args::Real...)` que toma una función `f` y sus correspondientes argumentos (como números `Real`) y devuelve los gradientes correspondientes:
```julia
function gradient(f, args::Real...)
ts = track.(args)
y = f(ts...)
y.grad = 1.0
accum!.(ts)
end
nothing # hide
```
que ahora podemos usar de la siguiente forma:
```julia
f(x,y) = x*y + sin(x)
gradient(f, 2.0, 3.0)
```
Como ejemplo, podemos encontrar un mínimo local de la función `g` (ligeramente modificada para que tenga un mínimo):
```julia
g(x,y) = y*y + sin(x)
using Plots
color_scheme = cgrad(:RdYlBu_5, rev=true)
contour(-4:0.1:4, -2:0.1:2, g, fill=true, c=color_scheme, xlabel="x", ylabel="y")
```
Podemos encontrar un mínimo local de $g$ partiendo de un punto inicial $(x_0,y_0)$ y dando pequeños pasos en la dirección opuesta al gradiente (el algoritmo típico de **descenso de gradiente**):
$$
\begin{align}
x_{i+1} &= x_i - \lambda \frac{\partial f}{\partial x_i} \\
y_{i+1} &= y_i - \lambda \frac{\partial f}{\partial y_i},
\end{align}
$$
donde $\lambda$ es la tasa de aprendizaje que debe ser ajustada manualmente.
Podemos implementar una función `descend` que realiza un paso de Descenso de Gradiente (GD) en una función `f` con un número arbitrario de entradas y especificando la tasa de aprendizaje $\lambda$:
```julia
function descend(f::Function, λ::Real, args::Real...)
Δargs = gradient(f, args...)
args .- λ .* Δargs
end
nothing # hide
```
Si se ejecuta un paso `descend` se obtendrán dos nuevas entradas con una salida menor para `g`:
```julia
julia> g(1.0, 1.0)
1.8414709848078965
julia> (x,y) = descend(g, 0.2, 1.0, 1.0)
(0.891939538826372, 0.8)
julia> g(x,y)
1.4182910542267546
```
Se puede minimizar `g` partiendo de un valor inicial y repitiendo el algoritmo anterior un número necesario de pasos.
A continuación se muestra un fragmento de código que realiza una serie de pasos `descend` en dos puntos iniciales diferentes y crea una animación de cada paso del algoritmo:
```julia
function minimize(f::Function, args::T...; niters=20, λ=0.01) where T<:Real
paths = ntuple(_->Vector{T}(undef,niters), length(args))
for i in 1:niters
args = descend(f, λ, args...)
@info f(args...)
for j in 1:length(args)
paths[j][i] = args[j]
end
end
paths
end
xs1, ys1 = minimize(g, 1.5, -2.4, λ=0.2, niters=34)
xs2, ys2 = minimize(g, 1.8, -2.4, λ=0.2, niters=16)
p1 = contour(-4:0.1:4, -2:0.1:2, g, fill=true, c=color_scheme, xlabel="x", ylabel="y")
scatter!(p1, [xs1[1]], [ys1[1]], mc=:black, marker=:star, ms=7, label="Minimum")
scatter!(p1, [xs2[1]], [ys2[1]], mc=:black, marker=:star, ms=7, label=false)
scatter!(p1, [-π/2], [0], mc=:red, marker=:star, ms=7, label="Initial Point")
scatter!(p1, xs1[1:1], ys1[1:1], mc=:black, label="GD Path", xlims=(-4,4), ylims=(-2,2))
@gif for i in 1:max(length(xs1), length(xs2))
if i <= length(xs1)
scatter!(p1, xs1[1:i], ys1[1:i], mc=:black, lw=3, xlims=(-4,4), ylims=(-2,2), label=false)
end
if i <= length(xs2)
scatter!(p1, xs2[1:i], ys2[1:i], mc=:black, lw=3, label=false)
end
p1
end
```
## AD inversa vectorizada
Una solución ingenua para utilizar el tipo de número `TrackedReal` para diferenciar funciones que operan sobre vectores es utilizar simplemente `Array{<:TrackedReal}`. Desafortunadamente, esto significa que tenemos que reemplazar las rápidas operaciones matriciales BLAS con nuestros propios métodos de multiplicación matricial que sepan cómo tratar con `TrackedReal`. Por ello, vamos a implementar una solución más inteligente a este problema:
```julia
using LinearAlgebra
Base.zero(::TrackedReal{T}) where T = TrackedReal(zero(T))
LinearAlgebra.adjoint(x::TrackedReal) = x
track(x::Array) = track.(x)
accum!(xs::Array{<:TrackedReal}) = accum!.(xs)
const VecTracked = AbstractVector{<:TrackedReal}
const MatTracked = AbstractMatrix{<:TrackedReal}
LinearAlgebra.dot(xs::VecTracked, ys::VecTracked) = mapreduce(*, +, xs, ys)
Base.:*(X::MatTracked, y::VecTracked) = map(x->dot(x,y), eachrow(X))
Base.:*(X::MatTracked, Y::MatTracked) = mapreduce(y->X*y, hcat, eachcol(Y))
Base.sum(xs::AbstractArray{<:TrackedReal}) = reduce(+,xs)
function reset!(x::TrackedReal)
x.grad = nothing
reset!.(keys(x.children))
x.children = Dict()
end
X = rand(2,3)
Y = rand(3,2)
function run()
Xv = track(X)
Yv = track(Y)
z = sum(Xv * Yv)
z.grad = 1.0
accum!(Yv)
end
```
Una ejecución de este método daría como resultado:
```julia
julia> using BenchmarkTools
julia> @benchmark run()
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 44.838 μs … 8.404 ms ┊ GC (min … max): 0.00% … 98.78%
Time (median): 48.680 μs ┊ GC (median): 0.00%
Time (mean ± σ): 53.048 μs ± 142.403 μs ┊ GC (mean ± σ): 4.61% ± 1.71%
Memory estimate: 26.95 KiB, allocs estimate: 392.
```
## AD inversa con `TrackedArray`
Para hacer uso de los métodos BLAS, mucho más rápidos, tenemos que implementar un tipo de array personalizado que descargará las pesadas multiplicaciones matriciales a los métodos matriciales normales. Empezaremos con un **REPL nuevo** y posiblemente un **nuevo fichero** que sólo contenga la definición de nuestro `TrackedReal`:
```julia
mutable struct TrackedReal{T<:Real}
data::T
grad::Union{Nothing,T}
children::Dict
end
track(x::Real) = TrackedReal(x, nothing, Dict())
```
Definimos un nuevo tipo, `TrackedArray`, subtipo de `AbstractArray{T,N}` y que contiene los tres campos: `data`, `grad`, y `children`. Además, definimos también `track(x::Array)`, y envíamos `size`, `length`, y `eltype` a `x.data`:
```julia
mutable struct TrackedArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
data::A
grad::Union{Nothing,A}
children::Dict
end
track(x::Array) = TrackedArray(x, nothing, Dict())
track(x::Union{TrackedArray,TrackedReal}) = x
for f in [:size, :length, :eltype]
eval(:(Base.$(f)(x::TrackedArray, args...) = $(f)(x.data, args...)))
end
# only needed for hashing in the children dict...
Base.getindex(x::TrackedArray, args...) = getindex(x.data,args...)
# pretty print TrackedArray
Base.show(io::IO, x::TrackedArray) = print(io, "Tracked $(x.data)")
Base.print_array(io::IO, x::TrackedArray) = Base.print_array(io, x.data)
```
A partir de estas definiciones, la creación de un `TrackedArray` debería funcionar así:
```julia
julia> track(rand(2,2))
2×2 Main.TrackedArray{Float64, 2, Matrix{Float64}}:
0.0369491 0.682941
0.317325 0.7683
```
Y podemos adaptar las funciones anteriores:
```julia
function accum!(x::Union{TrackedReal,TrackedArray})
if isnothing(x.grad)
x.grad = sum(λ(accum!(Δ)) for (Δ,λ) in x.children)
end
x.grad
end
```
Para aplicar la primera regla para `*`, es decir, la multiplicación de matrices, primero tendríamos que derivarla. En el caso de la multiplicación matricial general (que es una función $(\mathbb R^{n\times m}, \mathbb R^{m\times k}) \rightarrow \mathbb R^{n\times k}$) ya no estamos tratando con derivadas simples. No vamos a ver aquí los detalles técnicos de estas definiciones, pero debe indicarse que incluso una implementación sencilla es 10 veces más rápida que la implementación vista en la sección previa.
Para implementar una red neuronal completa necesitamos dos reglas más. Una para la no linealidad y otra para concatenar puntos de entrenamiento individuales en un lote:
```julia
σ(x::Real) = 1/(1+exp(-x))
σ(x::AbstractArray) = σ.(x)
function σ(x::TrackedArray)
z = track(σ(x.data))
d = z.data
x.children[z] = Δ -> Δ .* d .* (1 .- d)
z
end
function Base.hcat(xs::TrackedArray...)
y = track(hcat(data.(xs)...))
stops = cumsum([size(x,2) for x in xs])
starts = vcat([1], stops[1:end-1] .+ 1)
for (start,stop,x) in zip(starts,stops,xs)
x.children[y] = function (Δ)
δ = if ndims(x) == 1
Δ[:,start]
else
ds = map(_ -> :, size(x)) |> Base.tail |> Base.tail
Δ[:, start:stop, ds...]
end
δ
end
end
y
end
```
La ejecución de este último script producirá una animación que muestra cómo está aprendiendo la red.
Se puede ver una implementación completa de un AD basado en estas ideas [aquí](https://github.com/JuliaTeachingCTU/Scientific-Programming-in-Julia/blob/master/src/ReverseDiff.jl), y una implementación simple de una Red Neuronal que puede aprender una aproximación a la función `g` [aquí](https://github.com/JuliaTeachingCTU/Scientific-Programming-in-Julia/blob/master/docs/src/lecture_08/reverse-nn.jl).
# Fuentes y Referencias
* [Calculus on Computational Graphs: Backpropagation](https://colah.github.io/posts/2015-08-Backprop/). Blog de Colah. 31 Agosto de 2015.
* [Scientific Programming in Julia](https://juliateachingctu.github.io/Scientific-Programming-in-Julia)
* [ADCME](https://docs.juliahub.com/ADCME)
* [Bulding Computational Graph](https://tomroth.dev/compgraph1/). Blog Tom Roth. Febrero de 2023.
(insert ../menu.md.html here)