4. Deep Learning para la predicción de series temporales#
4.1. Redes neuronales#
Introducción
Las redes neuronales son sistemas de aprendizaje compuestos por neuronas conectadas en capas que ajustan sus conexiones para aprender. Tras un período de 25 años desde su inicio, las redes neuronales se convirtieron en la norma en el aprendizaje automático.
En un principio, dominaron durante una década, pero luego fueron superadas por máquinas de vectores de soporte. Sin embargo, desde 2010, las redes neuronales profundas se han vuelto populares gracias a mejoras en la tecnología y la disponibilidad de grandes conjuntos de datos, impulsando el campo del aprendizaje automático.
4.2. Gradiente descendente#
El método de gradiente descendente es uno de los mas ampliamente usados para la minimización iterativa de una función de costo diferenciable, \(J(\boldsymbol{\theta}),~\boldsymbol{\theta}\in\mathbb{R}^{l}\) (por ejemplo MSE). Como cualquier otra técnica iterativa, el método parte de una estimación inicial, \(\boldsymbol{\theta}^{(0)}\), y genera una sucesión \(\boldsymbol{\theta}^{(i)},~i=1,2,\dots,\) tal que:
\[ \boldsymbol{\theta}^{(i)}=\boldsymbol{\theta}^{(i-1)}+\mu_{i}\Delta\boldsymbol{\theta}^{(i)},~ i >0,~\mu_{i}>0. \]La diferencia entre cada método radica en la forma que \(\mu_{i}\) y \(\Delta\boldsymbol{\theta}^{(i)}\) son seleccionados. \(\Delta\boldsymbol{\theta}^{(i)}\) es conocido como la dirección de actualización o de búsqueda. La sucesión \(\mu_{i}\) es conocida como el tamaño o longitud de paso en la \(i\)-ésima iteración, estos valores pueden ser constantes o cambiar.
En el método de gradiente descendente, la selección de \(\Delta\boldsymbol{\theta}^{(i)}\) es realizada para garantizar que \(J(\boldsymbol{\theta}^{(i)})<J(\boldsymbol{\theta}^{(i-1)})\), excepto en el minimizador \(\boldsymbol{\theta}_{\star}\).
Fig. 4.1 Función de coste en el espacio de parámetros bidimensional. Fuente [Theodoridis, 2020].#
Suponga que en la iteración \(i-1\) el valor \(\boldsymbol{\theta}^{(i-1)}\) ha sido obtenido
Nótese que seleccionando la dirección tal que \(\nabla^{T}J(\boldsymbol{\theta}^{(i-1)})\Delta\boldsymbol{\theta}^{(i)}<0\), garantizará que \(J(\boldsymbol{\theta}^{(i)})=J(\boldsymbol{\theta}^{(i-1)}+\mu_{i}\Delta\boldsymbol{\theta}^{(i)})<J(\boldsymbol{\theta}^{(i-1)})\). Tal selección de \(\Delta\boldsymbol{\theta}^{(i)}\) y \(\nabla J(\boldsymbol{\theta}^{(i-1)})\) debe formar un ángulo obtuso (\(90^{\circ}<\phi<180^{\circ}\)). Las curvas de nivel asociadas a \(J(\boldsymbol{\theta})\) pueden tomar cualquier forma, la cual va a depender de como está definido \(J(\boldsymbol{\theta})\).
\(J(\boldsymbol{\theta})\) se supone diferenciable, por lo tanto, las curvas de nivel o contornos deben ser suaves y aceptar un plano tangente en cualquier punto. Además, de los cursos de cálculo sabemos que el vector gradiente \(\nabla J(\boldsymbol{\theta})\) es perpendicular al plano tangente (recta tangente) a la correspondiente curva de nivel en el punto \(\boldsymbol{\theta}\).
Nótese que seleccionando la dirección de búsqueda \(\Delta\boldsymbol{\theta}^{(i)}\) que forma un angulo obtuso con el gradiente, se coloca a \(\boldsymbol{\theta}^{(i-1)}+\mu_{i}\Delta\boldsymbol{\theta}^{(i)}\) en un punto sobre el contorno el cual corresponde a un valor menor que \(J(\boldsymbol{\theta})\).
Dos problemas surgen ahora:
Escoger la mejor dirección de búsqueda
Calcular que tan lejos es aceptable un movimiento a traves de esta dirección.
Fig. 4.2 El vector gradiente en un punto \(\boldsymbol{\theta}\) es perpendicular al plano tangente (línea punteada) en la curva de nivel que cruza \(\boldsymbol{\theta}\). La dirección de descenso forma un ángulo obtuso, \(\phi\), con el vector gradiente. Fuente [Theodoridis, 2020]#
Nótese que si \(\mu_{i}\|\Delta\boldsymbol{\theta}^{(i)}\|\) es demasiado grande, entonces el nuevo punto puede ser colocado en un contorno correspondiente a un valor mayor al del actual contorno.
Fig. 4.3 Las correspondientes curvas de nivel para la función de coste, en el plano bidimensional. Nótese que a medida que nos alejamos del valor óptimo, \(\boldsymbol{\theta}_{\star}\), los valores de \(c\) aumentan. Fuente [Theodoridis, 2020].#
Para abordar (1), supongamos que \(\mu_{i}=1\) y buscamos todos los vectores \(\boldsymbol{z}\) con norma Euclidiana unitaria, con inicio (cola) en \(\boldsymbol{\theta}^{(i-1)}\). Entonces, para todas las posibles direcciones, la que entrega el valor más negativo del producto interno, \(\nabla^{T}J(\boldsymbol{\theta}^{(i-1)})z\), es aquella de gradiente negativo
Centrando \(\boldsymbol{\theta}^{(i-1)}\) en la bola con norma Euclideana uno. De todos los vectores con norma unitaria y origen en \(\boldsymbol{\theta}^{(i-1)}\), seleccionamos aquel que apunta en la dirección negativa del gradiente. Por lo tanto, para todos los vectores con norma Euclidiana 1, la dirección de descenso mas pronunciada coincide con la dirección del gradiente descendente, negativo, y la correspondiente actualización recursiva se convierte en
Fig. 4.4 Representación del gradiente negativo, el cual conduce a la máxima disminución de la función de coste.#
La selección de \(\mu_{i}\) debe ser realizada de tal forma que garantice convergencia de la secuencia de minimización. Nótese que el algoritmo puede oscilar en torno al mínimo sin converger, si no seleccionamos la dirección correcta. La selección de \(\mu_{i}\) dependerá de la convergencia a cero del error entre \(\boldsymbol{\theta}^{(i)}\) y el mínimo real en forma de serie geométrica.
Por ejemplo, para el caso de la función de coste del error cuadrático medio, la longitud de paso está dada por: \(0<\mu<2/\lambda_{\max}\), donde \(\lambda_{\max}\) el máximo eigenvalor de la matriz de covarianza \(\Sigma_{x}=\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^{T}]\), donde \(J(\boldsymbol{\theta})=\text{E}[(y-\boldsymbol{\theta}^{T}\boldsymbol{x})^{2}]\) (ver [Theodoridis, 2020]) (
bonus).
4.3. El perceptrón#
Nuestro punto de partida será considerar el problema simple de una tarea de clasificación conformada por dos clases linealmente separables. En otras palabras, dado un conjunto de muestras de entrenamiento, \((y_{n}, \boldsymbol{x}_{n})\), \(n=1,2,\dots,N\), con \(y_{n}\in\{-1,+1\},~\boldsymbol{x}_{n}\in\mathbb{R}^{l}\), suponemos que existe un hiperplano
\[ \boldsymbol{\theta}_{\star}^{T}\boldsymbol{x}=0, \]tal que,
\[\begin{split} \begin{cases} \boldsymbol{\theta}_{\star}^{T}\boldsymbol{x}&>0,\quad\text{si}\quad\boldsymbol{x}\in\omega_{1}\\ \boldsymbol{\theta}_{\star}^{T}\boldsymbol{x}&<0,\quad\text{si}\quad\boldsymbol{x}\in\omega_{2} \end{cases} \end{split}\]en otras palabras, dicho hiperplano clasifica correctamente todos los puntos del conjunto de entrenamiento.
Para simplificar, el término de sesgo del hiperplano ha sido absorbido en \(\boldsymbol{\theta}_{\star}\) después de extender la dimensionalidad del espacio de entrada en uno. El objetivo ahora es desarrollar un algoritmo que calcule iterativamente un hiperplano que clasifique correctamente todos los patrones de ambas clases. Para ello, se adopta una función de costo.
Sea \(\boldsymbol{\theta}\) la estimación del vector de parámetros desconocidos, disponible en la actual iteración. Entonces hay dos posibilidades. La primera es que todos los puntos estén clasificados correctamente; esto significa que se ha obtenido una solución. La otra alternativa es que \(\boldsymbol{\theta}\) clasifique correctamente algunos de los puntos y el resto estén mal clasificados.
Costo perceptrón
Sea \(\mathcal{Y}\) el conjunto de todas las muestras mal clasificadas. La función de costo, perceptrón se define como
donde
Nótese que la función de costo es no negativa. Dado que la suma es sobre los puntos mal clasificados, si \(\boldsymbol{x}_{n}\in\omega_{1}~(\omega_{2}),~\) entonces \(\boldsymbol{\theta}^{T}\boldsymbol{x}_{n}\leq (\geq)~0\), entregando así un producto \(-y_{n}\boldsymbol{\theta}^{T}\boldsymbol{x}_{n}\geq0\).
La función de costo es cero, si no existen puntos mal clasificados, esto es, \(\mathcal{Y}=\emptyset\). La función de costo perceptrón no es diferenciable en todos los puntos, es lineal por tramos. Si reescribimos \(J(\boldsymbol{\theta})\) en una forma ligeramente diferente:
Nótese que esta es una función lineal con respeto a \(\boldsymbol{\theta}\), siempre que el conjunto de puntos mal clasificados permanezca igual. Además, nótese que ligeros cambios del valor \(\boldsymbol{\theta}\) corresponden a cambios de posición del respectivo hiperplano.
Como consecuencia, existirá un punto donde el número de muestras mal clasificadas en \(\mathcal{Y}\), repentinamente cambia; este es el tiempo donde una muestra en el conjunto de entrenamiento cambia su posición relativa con respecto al hiperplano en movimiento, y en consecuencia, el conjunto \(\mathcal{Y}\) es modificado. Después de este cambio, el conjunto, \(J(\boldsymbol{\theta})\), corresponderá a una nueva función lineal.
El algoritmo perceptrón
A partir del método de subgradientes se puede verificar fácilmente que, iniciando desde un punto arbitrario, \(\boldsymbol{\theta}^{(0)}\), el siguiente método iterativo,
converge después de un número finito de pasos. La sucesión de parámetros \(\mu_{i}\) es seleccionada adecuadamente para garantizar convergencia.
Nótese que usando el método de subgradiente se tiene que
Otra versión del algoritmo considera una muestra por iteración en un esquema cíclico, hasta que el algoritmo converge. Denotemos por \(y_{(i)}\), \(\boldsymbol{x}_{i},~i\in\{1,2,\dots,N\}\) los pares de entrenamiento presentados al algoritmo en la iteración \(i\)-ésima. Entonces, la iteración de actualización se convierte en:
Esto es, partiendo de una estimación inicial de forma random, inicializando \(\boldsymbol{\theta}^{(0)}\) con algunos valores pequeños, testeamos cada una de las muestras, \(\boldsymbol{x}_{n},~n=1,2,\dots,N\). Cada vez que una muestra es mal clasificada, se toma acción por medio de la regla perceptrón para una corrección. En otro caso, ninguna acción es requerida. Una vez que todas las muestras han sido consideradas, decimos que una
época (epoch)ha sido completada.
Observación
Si no se obtiene convergencia, todas las muestras son reconsideradas en una segunda época, y así sucesivamente. La versión de este algoritmo es conocida como esquema
pattern-by-pattern. Algunas veces también es referido como el algoritmo online. Nótese que el número total de datos muestrales es fijo, y que el algoritmo los considera en forma cíclica, época por época (epoch-by-epoch).Después de un número finito de épocas, se garantiza que el algoritmo es convergente. Nótese que para obtener dicha convergencia, la sucesión \(\mu_{i}\) debe ser seleccionada apropiadamente. Sin embargo, para el caso del algoritmo perceptrón, la convergencia es garantizada (\(J\) convexa), aun cuando \(\mu_{i}\) es una constante positiva, \(\mu_{i}=\mu>0\), usualmente tomado igual a uno.
La formulación en (4.1) es conocida también como la filosofía de aprendizaje
reward-punishment. Si la actual estimación es exitosa en la predicción de la clase del respectivo patrón, ninguna acción es tomada (reward), en otro caso, el algoritmo es obligado a realizar una actualización (punishment).
Fig. 4.5 El punto \(x\) está mal clasificado por la línea roja. La regla perceptrón gira el hiperplano hacia el punto \(x\), para intentar incluirlo en el lado correcto del nuevo hiperplano y clasificarlo correctamente. El nuevo hiperplano está definido por \(θ^{(i)}\) y se muestra con la línea negra. Fuente [Theodoridis, 2020].#
La Fig. 4.5 ofrece una interpretación geométrica de la regla del perceptrón. Supongamos que la muestra \(\boldsymbol{x}\) está mal clasificada por el hiperplano, \(\boldsymbol{\theta}^{(i-1)}\). Como sabemos, por geometría analítica, \(\boldsymbol{\theta}^{(i-1)}\) corresponde a un vector que es perpendicular al hiperplano que está definido por este vector.
Como \(\boldsymbol{x}\) se encuentra en el lado \((-)\) del hiperplano y está mal clasificado, pertenece a la clase \(\omega_{1}\); asumiendo \(\mu = 1\), la corrección aplicada por el algoritmo es
\[\begin{split} \\[1mm] \boldsymbol{\theta}^{(i)}=\boldsymbol{\theta}^{(i-1)}+\boldsymbol{x}, \end{split}\]y su efecto es girar el hiperplano en dirección a \(\boldsymbol{x}\) para colocarlo en el lado \((+)\) del nuevo hiperplano, que está definido por la estimación actualizada \(\boldsymbol{\theta^{(i)}}\). El algoritmo perceptrón en su modo de funcionamiento patrón por patrón (pattern-by-pattern) se resume en el siguiente algoritmo.
Algorithm 4.1 (Algoritmo perceptrón pattern-by-pattern)
Inicialización
Inicializar \(\boldsymbol{\theta}^{(0)}\); usualmente, de forma
random, pequeñoSeleccionar \(\mu\); usualmente
establecido como 1\(i=1\)
Repeat Cada iteración corresponde a un epoch
counter = 0; Contador del número de actualizaciones porepochFor \(n=1,2,\dots,N\) Do Para cada
epoch, todas las muestras son presentadas una vezIf(\(y_{n}\boldsymbol{x}_{n}^{T}\leq0\)) Then
\(\boldsymbol{\theta}^{(i)}=\boldsymbol{\theta}^{(i-1)}+\mu y_{n}\boldsymbol{x}_{n}\)
\(i=i+1\)
counter = counter + 1
End For
Until counter = 0
Una vez que el algoritmo perceptrón se ha ejecutado y converge, tenemos los pesos, \(\theta_{i},~i = 1,2,\dots,l\), de las sinapsis de la neurona/perceptrón asociada, así como el término de sesgo \(\theta_{0}\). Ahora se pueden utilizar para clasificar patrones desconocidos. Las características \(x_{i}, i = 1, 2,\dots,l\), se aplican a los nodos de entrada.
A su vez, cada característica se multiplica por la sinapsis respectiva (peso), y luego se añade el término de sesgo en su combinación lineal. El resultado de esta operación pasa por una función no lineal, \(f\), conocida como función de activación. Dependiendo de la forma de la no linealidad, se producen diferentes tipos de neuronas. La más clásica conocida como neurona McCulloch-Pitts, la función de activación es la de Heaviside, es decir,
Fig. 4.6 Arquitectura básica de neuronas/perceptrones.#
En la arquitectura básica de neuronas/perceptrones, las características de entrada se aplican a los nodos de entrada y se ponderan por los respectivos pesos que definen las sinapsis. A continuación se añade el término de sesgo en su combinación lineal y el resultado es empujado a través de la no linealidad.
En la neurona McCulloch-Pitts, la salida es 1 para los patrones de la clase \(\omega_{1}\) o 0 para la clase \(\omega_{2}\). La suma y la operación no lineal se unen para simplificar el gráfico.
Para las capas ocultas, la función de activación tangente hiperbólica suele funcionar mejor que la sigmoidea logística. Tanto la función sigmoid como tanh pueden hacer que el modelo sea más susceptible a los problemas durante el entrenamiento, a través del llamado problema de los
gradientes desvanecientes. Los modelos modernos de redes neuronales con arquitecturas comunes, como MLP y CNN, harán uso de la función de activación ReLU, o extensiones.Las redes recurrentes suelen utilizar funciones de activación tanh o sigmoid, o incluso ambas. Por ejemplo, la LSTM suele utilizar la activación sigmoid para las conexiones recurrentes y la activación tanh para la salida.
Fig. 4.8 Selección de función de activación para output layers. (Fuente [Brownlee and Mastery, 2017]).#
Si su problema es de regresión, debe utilizar una función de activación lineal. Si su problema es de clasificación, el modelo predice la probabilidad de pertenencia a una clase, que se puede convertir en una etiqueta de clase mediante redondeo (para sigmoid) o argmax (para softmax).
Activación Lineal y Problemas de Regresión
En las Redes Neuronales Artificiales (ANNs), la función de activación lineal se utiliza comúnmente en la capa de salida para problemas de regresión y series de tiempo debido a las siguientes razones:
Producción de Salidas ContinuasLa regresión y las predicciones de series de tiempo generalmente requieren que el modelo produzca salidas continuas, no discretas. Una función de activación lineal permite que la salida de la red neuronal sea cualquier valor real, lo que es adecuado para este tipo de problemas. Por ejemplo, si estamos prediciendo precios de viviendas o temperaturas futuras, necesitamos que las salidas sean valores continuos, no limitados a un rango específico.
Evitar SaturaciónFunciones de activación no lineales como la sigmoid o tanh tienen salidas acotadas (entre 0 y 1 para sigmoide, y entre -1 y 1 para tanh). Estas funciones pueden saturarse, es decir, para valores grandes de entrada, los gradientes pueden volverse muy pequeños, lo que dificulta el entrenamiento del modelo debido al problema del gradiente desvaneciente. Una función de activación lineal no tiene este problema, ya que la salida es directamente proporcional a la entrada.
Preservación de la Escala de la SalidaEn tareas de regresión y predicción de series de tiempo, a menudo es importante que la escala de las salidas se mantenga consistente con la escala de los valores objetivo. La función de activación lineal, definida como, \(\boldsymbol{f(x)=x}\) no cambia la escala de los valores, permitiendo que la red neuronal modele relaciones directamente proporcionales entre las entradas y las salidas.
4.4. Redes Totalmente Conectadas#
Para resumir de manera más formal el tipo de operaciones que tienen lugar en una red totalmente conectada, centrémonos en, por ejemplo, la capa \(r\) de una red neuronal multicapa y supongamos que está formada por \(k_{r}\) neuronas. El vector de entrada a esta capa está formado por las salidas de los nodos de la capa anterior, que se denomina \(\boldsymbol{y}^{r-1}\).
Sea \(\boldsymbol{\theta}_{j}^{r}\) el vector de los pesos sinápticos, incluido el término de sesgo, asociado a la neurona \(j\) de la capa \(r\), donde \(j = 1,2,\dots, k_{r}\). La dimensión respectiva de este vector es \(k_{r-1} + 1\), donde \(k_{r-1}\) es el número de neuronas de la capa anterior, \(r-1\), y el aumento en 1 representa el término de sesgo. Entonces las operaciones realizadas, antes de la no linealidad, son los productos internos
Colocando todos los valores de salida en un vector \(\boldsymbol{z}^{r}=[z_{1}^{r}, z_{2}^{r},\dots,z_{k_{r}}^{r}]^{T}\), y agrupando todos los vectores sinápticos como filas, una debajo de la otra, en una matriz, podemos escribir colectivamente
El vector de las salidas de la \(r\)-ésima capa oculta, después de empujar cada \(z_{i}^{r}\) a través de la no linealidad \(f\), está finalmente dado por
La siguiente figura describe como es creada la \(j\)-ésima columna de la red neuronal full conectada
Fig. 4.9 \(j\)-ésimo elemento de la capa \(r\) de la red totalmente conectada.#
Observación
La notación anterior significa que \(f\) actúa sobre cada uno de los respectivos componentes del vector, individualmente, y la extensión del vector en uno es para dar cuenta de los términos de sesgo en la práctica estándar. Para redes grandes, con muchas capas y muchos nodos por capa, este tipo de conectividad resulta ser muy costoso en términos del número de parámetros (pesos), que es del orden de \(k_{r}k_{r-1}\).
Por ejemplo, si \(k_{r-1} = 1000\) y \(k_{r} = 1000\), esto equivale a un orden de 1 millón de parámetros. Tenga en cuenta que este número es la contribución de los parámetros de una sola de las capas. Sin embargo, un gran número de parámetros hace que una red sea vulnerable al sobreajuste, cuando se trata de entrenamiento
Se pueden emplear las llamadas técnicas de reparto de pesos, en las que un conjunto de parámetros es compartido entre un número de conexiones, a través de restricciones adecuadamente incorporadas. Las redes neuronales recurrentes y convolucionales pertenecen a esta familia de redes de peso compartido. En una red convolucional, las convoluciones sustituyen a las operaciones de producto interno, lo que permite un reparto de pesos importante que conduce a una reducción sustancial del número de parámetros.
4.5. El Algoritmo De Backpropagation#
Una red neuronal considera una función paramétrica no lineal, \(\hat{y} = f_{\boldsymbol{\theta}}(\boldsymbol{x})\), donde \(\boldsymbol{\theta}\) representa todos los pesos/sesgo presentes en la red. Por lo tanto, el entrenamiento de una red neuronal no parece ser diferente del entrenamiento de cualquier otro modelo de predicción paramétrico.
Todo lo que se necesita es (a) un conjunto de muestras de entrenamiento, (b) una función de pérdida \(\mathcal{L}(y, \hat{y})\), y (c) un esquema iterativo, por ejemplo, el gradiente descendente, para realizar la optimización de la función de coste asociada (pérdida empírica).
La dificultad del entrenamiento de las redes neuronales radica en su estructura multicapa que complica el cálculo de los gradientes, que intervienen en la optimización. Además, la neurona McCulloch-Pitts se basa en la función de activación discontinua de Heaviside, que no es diferenciable.
La
neurona sigmoidea logística: Una posibilidad es adoptar la función sigmoidea logística, es decir,
Nótese que cuanto mayor sea el valor del parámetro \(a\), la gráfica correspondiente se acerca más a la de la función de Heaviside (ver Fig. 4.10).
Fig. 4.10 La función sigmoidea logística para diferentes valores del parámetro \(a\).#
Otra posibilidad sería utilizar la función,
\[ f(z)=a\tanh\left(\frac{cz}{2}\right), \]donde \(c\) y \(a\) son parámetros de control. El gráfico de esta función se muestra en la Fig. 4.11. Nótese que a diferencia de la sigmoidea logística, ésta es una función no simétrica, es decir, \(f(-z)=-f(z)\). Ambas son también conocidas como funciones de reducción, porque limitan la salida a un rango finito de valores.
Fig. 4.11 Función de reducción de la tangente hiperbólica para \(a = 1.7\) y \(c = 4/3\).#
Recordemos que la regla de actualización del algoritmo gradiente descendente, en su versión unidimensional se convierte en
\[\begin{split} \\[1mm] \theta(new)=\theta(old)-\mu\left.\frac{d J}{d\theta}\right|_{\theta(old)}, \end{split}\]y las iteraciones parten de un punto inicial arbitrario, \(\theta^{(0)}\). Si en la iteración actual el algoritmo está digamos, en el punto \(\theta(old) = \theta_{1}\), entonces se moverá hacia el mínimo local, \(\theta_{l}\). Esto se debe a que la derivada del coste en \(\theta_{1}\) es igual a la tangente de \(\phi_{1}\), que es negativa (el ángulo es obtuso) y la actualización, \(\theta(new)\), se moverá a la derecha, hacia el mínimo local, \(\theta_{l}\).
Fig. 4.12 Función no convexa global, con mínimos locales y puntos de silla.#
La elección del tamaño del paso, \(\mu\), es crítica para la convergencia del algoritmo. En problemas reales en espacios multidimensionales, el número de mínimos locales puede ser grande, por lo que el algoritmo puede converger a uno local. Sin embargo, esto no es necesariamente una mala noticia.
Si este mínimo local es lo suficientemente profundo, es decir, si el valor de la función de coste en este punto, por ejemplo, \(J(\theta_{l})\), no es mucho mayor que el alcanzado en el mínimo global, es decir, \(J(\theta_{g})\), la convergencia a dicho mínimo local puede corresponder a una buena solución.
4.6. El Esquema De Backpropagation Para Gradiente descendente#
Habiendo adoptado una función de activación diferenciable, estamos listos para proceder a desarrollar el esquema iterativo de gradiente descendente para la minimización de la función de coste. Formularemos la tarea en un marco general.
Sea \((\boldsymbol{y}_{n}, \boldsymbol{x}_{n}), n = 1, 2,\dots, N\), el conjunto de muestras de entrenamiento. Nótese que hemos asumido output multivariado. Suponemos que la red consta de \(L\) capas, \(L-1\) capas ocultas y una capa de salida. Cada capa consta de \(k_{r}, r = 1, 2,\dots, L\), neuronas. Así, los vectores de salida (objetivo/deseado) son
Para ciertas derivaciones matemáticas, también denotamos el número de nodos de entrada como \(k_{0}\); es decir \(k_{0} = l\), donde \(l\) es la dimensionalidad del espacio de características de entrada.
Sea \(\boldsymbol{\theta}_{j}^{r}\) el vector de los pesos sinápticos asociados a la \(j\)-ésima neurona de la \(r\)-ésima capa, con \(j = 1, 2,\dots, k_{r}\) y \(r = 1, 2,\dots,L\), donde el término de sesgo se incluye en \(\boldsymbol{\theta}_{j}^{r}\) , es decir,
Los pesos sinápticos enlazan la neurona respectiva con todas las neuronas de la capa \(k_{r-1}\) (véase la Fig. 4.13). El paso iterativo básico para el esquema de gradiente descendente se escribe como
El parámetro \(\mu\) es el tamaño de paso definido por el usuario (también puede depender de la iteración) y \(J\) denota la función de coste.
Fig. 4.13 Enlaces y las variables asociadas de la \(j\)-ésima neurona en la \(r\)-ésima capa. \(y_{k}^{r-1}\) es la salida de la \(k\)-ésima neurona de la \((r - 1)\)-ésima capa y \(\theta_{jk}^{r}\) es el peso respectivo que conecta estas dos neuronas.#
Forward/Backward
Las Ecuaciones de actualización (4.3) comprenden el esquema de gradiente descendente para la optimización. Como se ha dicho anteriormente, la dificultad de las redes neuronales feed-forward surge de su estructura multicapa. Para calcular los gradientes en la Ecuación (4.3), para todas las neuronas, en todas las capas, se deben seguir dos pasos en su cálculo
Forward computations: Para un vector de entrada dado \(\boldsymbol{x}_{n}, n = 1, 2,\dots, N\), se utilizan las estimaciones actuales de los parámetros (pesos sinápticos) (\(\boldsymbol{\theta}_{j}^{r}(old)\)) y calcula todas las salidas de todas las neuronas en todas las capas, denotadas como \(y_{nj}^{r}\); en la Fig. 4.13, se ha suprimido el índice \(n\) para no afectar la notación.Backward computations: Utilizando las salidas neuronales calculadas anteriormente junto con los valores objetivo conocidos, \(y_{nk}\), de la capa de salida, se calculan los gradientes de la función de coste. Esto implica \(L\) pasos, es decir, tantos como el número de capas. La secuencia de los pasos algorítmicos se indica a continuación:Calcular el gradiente de la función de coste con respecto a los parámetros de las neuronas de la última capa, es decir, \(\displaystyle{\frac{\partial J}{\partial\boldsymbol{\theta}_{j}^{L}}, j = 1, 2,\dots, k_{L}}\).
For \(r = L-1\) to \(1\), Do
Calcular los gradientes con respecto a los parámetros asociados a las neuronas de la \(r\)-ésima capa, es decir, \(\displaystyle{\frac{\partial J}{\partial\boldsymbol{\theta}_{k}^{r}}, k= 1, 2,\dots, k_{r}}\) basado en todos los gradientes \(\displaystyle{\frac{\partial J}{\partial\boldsymbol{\theta}_{j}^{r+1}}}\), \(j= 1, 2,\dots, k_{r+1}\), con respecto a los parámetros de la capa \(r + 1\) que se han calculado en el paso anterior.
End For
El esquema de cálculo hacia atrás backpropagation es una aplicación directa de la regla de la cadena para las derivadas, y comienza con el paso inicial de calcular las derivadas asociadas a la última capa (de salida), que resulta ser sencillo.
Observación
La última derivada en el algoritmo backpropagation es sencilla porque se calcula directamente a partir de la función de pérdida y las salidas de la red, sin necesidad de propagar errores hacia atrás a través de múltiples capas
A continuación, el algoritmo “fluye” hacia atrás en la jerarquía de capas. Esto se debe a la naturaleza de la red multicapa, donde las salidas, capa tras capa, se forman como funciones de funciones. En efecto, centrémonos en la salida \(y_{k}^{r}\) de la neurona \(k\) en la capa \(r\). Entonces tenemos
\[ y_{k}^{r}=f(\boldsymbol{\theta}_{k}^{r^T}\boldsymbol{y}^{r-1}),\quad k=1,2,\dots, k_{r}, \]donde \(\boldsymbol{y}^{r-1}\) es el vector (ampliado) que comprende todas las salidas de la capa anterior, \(r-1\), y \(f\) denota la no-linealidad.
De acuerdo con lo anterior, la salida de la \(j\)-ésima neurona en la siguiente capa viene dada por
\[\begin{split} y_{j}^{r+1}=f(\boldsymbol{\theta}_{j}^{r+1^T}\boldsymbol{y}^{r})=f\left(\boldsymbol{\theta}_{j}^{r+1^{T}} \begin{bmatrix} 1\\ f(\Theta^{r}\boldsymbol{y}^{r-1}) \end{bmatrix} \right), \end{split}\]donde \(\Theta^{r}:=[\boldsymbol{\theta}_{1}^{r}, \boldsymbol{\theta}_{2}^{r},\dots,\boldsymbol{\theta}_{k_{r}}^{r}]^{T}\) denota la matriz cuyas columnas corresponden al vector de pesos en el layer \(r\).
Nótese que obtuvimos evaluación de “una función interna bajo una función externa”. Claramente, esto continúa a medida que avanzamos en la jerarquía. Esta estructura de evaluación de funciones internas por funciones externas, es el subproducto de la naturaleza multicapa de las redes neuronales, la cual es una operación altamente no lineal, que da lugar a la dificultad de calcular los gradientes, a diferencia de otros modelos, como por ejemplo SVM.
Sin embargo, se puede observar fácilmente que
el cálculo de los gradientes con respecto a los parámetros que definen la capa de salida no plantea ninguna dificultad. En efecto, la salida de la \(j\)-ésima neurona de la última capa (que es en realidad la respectiva estimación de salida actual) se escribe como:
Dado que \(\boldsymbol{y}^{L-1}\) es conocido, después de los cálculos durante el paso adelante, tomando la derivada con respecto a \(\boldsymbol{\theta}_{j}^{L}\) es sencillo; no hay ninguna operación de función sobre función. Por esto es que empezamos por la capa superior y luego nos movemos hacia atrás. Debido a su importancia histórica, se dará la derivación completa del algoritmo backpropagation.
Para la derivación detallada del algoritmo backpropagation, se adopta como ejemplo la función de pérdida del error cuadrático, es decir
(4.4)#\[ J(\boldsymbol{\theta})=\sum_{n=1}^{N}J_{n}(\boldsymbol{\theta})\quad\text{y}\quad J_{n}(\boldsymbol{\theta})=\frac{1}{2}\sum_{k=1}^{k_{L}}(\hat{y}_{nk}-y_{nk})^{2}, \]donde \(\hat{y}_{nk},~k=1,2,\dots,k_{L}\), son las estimaciones proporcionadas en los correspondientes nodos de salida de la red. Las consideraremos como los elementos de un vector correspondiente, \(\hat{\boldsymbol{y}}_{n}\).
4.7. Cálculo de gradientes#
Sea \(z_{nj}^{r}\) la salida del combinador lineal de la \(j\)-ésima neurona en la capa \(r\) en el instante de tiempo \(n\), cuando se aplica el patrón \(\boldsymbol{x}_{n}\) en los nodos de entrada (véase la Fig. 4.13). Entonces, para \(n, j\) fijos, podemos escribir
(4.5)#\[ z_{nj}^{r}=\sum_{m=1}^{k_{r-1}}\theta_{jm}^{r}y_{nm}^{r-1}+\theta_{j0}^{r}=\sum_{m=0}^{k_{r-1}}\theta_{jm}^{r}y_{nm}^{r-1}=\boldsymbol{\theta}_{j}^{r^{T}}\boldsymbol{y}_{n}^{r-1}, \]donde por definición
\[ \boldsymbol{y}_{n}^{r-1}:=[1, y_{n1}^{r-1},\dots, y_{nk_{r-1}}^{r-1}]^{T}, \]y \(y_{n0}^{r}\equiv 1,~\forall~r, n\) y \(\theta_{j}^{r}\) ha sido definido en la Ecuación (4.2).
Para las neuronas de la capa de salida \(r=L,~y_{nm}^{L}=\hat{y}_{nm},~m=1,2,\dots, k_{L}\), y para \(r=1\), tenemos \(y_{nm}^{1}=x_{nm},~m=1,2,\dots, k_{1}\); esto es, \(y_{nm}^{1}\) se fijan iguales a los valores de las características de entrada.
Por lo tanto, podemos escribir ahora
Definamos
Entonces la Ecuación (4.3) puede escribirse como
4.8. Cálculo de \(\delta_{nj}^{r}\)#
Este es el cálculo principal del algoritmo backpropagation. Para el cálculo de los gradientes, \(\delta_{nj}^{r}\), se comienza en la última capa, \(r = L\), y se procede hacia atrás, hacia \(r = 1\); esta “filosofía” justifica el nombre dado al algoritmo.
\(r=L\): Tenemos que
\[ \delta_{nj}^{L}:=\frac{\partial J_{n}}{\partial z_{nj}^{L}}. \]Para la función de pérdida del error al cuadrado,
\[ J_{n}=\frac{1}{2}\sum_{k=1}^{k_{L}}\left(\hat{y}_{nk}-y_{nk}\right)^{2}=\frac{1}{2}\sum_{k=1}^{k_{L}}\left(f(z_{nk}^{L})-y_{nk}\right)^{2}. \]Por lo tanto,
(4.7)#\[\begin{split} \begin{align*} \delta_{nj}^{L}=\frac{\partial}{\partial z_{nj}^{L}}\left(\frac{1}{2}\sum_{k=1}^{k_{L}}\left(f(z_{nk}^{L})-y_{nk}\right)^{2}\right)&=(f(z_{nj}^{L})-y_{nj})f'(z_{nj}^{L})\\ &=(\hat{y}_{nj}-y_{nj})f'(z_{nj}^{L})=e_{nj}f'(z_{nj}^{L}), \end{align*} \end{split}\]donde \(j=1,2,\dots, k_{L}\), \(f'\) denota la derivada de \(f\) y \(e_{nj}\) es el error asociado con el \(j\)-ésima output en el tiempo \(n\). Nótese que para el último layer, el cálculo del gradiente, \(\delta_{nj}^{L}\) es sencillo.
\(r<L\): Debido a la dependencia sucesiva entre las capas, el valor de \(z_{nj}^{r-1}\) influye en todos los valores \(z_{nk}^{r},~k = 1, 2,\dots, k_{r}\), de la capa siguiente. Empleando la regla de la cadena para la diferenciación, obtenemos, para \(r=L, L-1,\dots, 2\):
(4.8)#\[ \delta_{nj}^{r-1}=\frac{\partial J_{n}}{\partial z_{nj}^{r-1}}=\sum_{k=1}^{k_{r}}\frac{\partial J_{n}}{\partial z_{nk}^{r}}\frac{\partial z_{nk}^{r}}{\partial z_{nj}^{r-1}}, \]o
(4.9)#\[ \delta_{nj}^{r-1}=\sum_{k=1}^{k_{r}}\delta_{nk}^{r}\frac{\partial z_{nk}^{r}}{\partial z_{nj}^{r-1}}. \]Además, usando Ecuación (4.5) se tiene,
\[ \frac{\partial z_{nk}^{r}}{\partial z_{nj}^{r-1}}=\frac{\partial}{\partial z_{nj}^{r-1}}\left(\sum_{m=0}^{k_{r-1}}\theta_{km}^{r}y_{nm}^{r-1}\right), \]donde \(\textcolor{blue}{y_{nm}^{r-1}=f(z_{nm}^{r-1})}=f(\boldsymbol{\theta}_{m}^{r-1^{T}}\boldsymbol{y}_{n}^{r-2})\). Entonces,
\[ \frac{\partial z_{nk}^{r}}{\partial z_{nj}^{r-1}}=\theta_{kj}^{r}f'(z_{nj}^{r-1}), \]y combinando las Ecuaciones (4.8) y (4.9), obtenemos la regla recursiva
\[ \delta_{nj}^{r-1}=\left(\sum_{k=1}^{k_{r}}\delta_{nk}^{r}\theta_{kj}^{r}\right)f'(z_{nj}^{r-1}). \]
Manteniendo la misma notación en la Ecuación (4.7), definimos
\[\begin{split} \\[1mm] e_{nj}^{r-1}:=\sum_{k=1}^{k_{r}}\delta_{nk}^{r}\theta_{kj}^{r}, \end{split}\]y finalmente obtenemos,
(4.10)#\[ \delta_{nj}^{r-1}=e_{nj}^{r-1}f'(z_{nj}^{r-1}). \]
El único cálculo que queda es la derivada de \(f\). Para el caso de la función sigmoidea logística se demuestra fácilmente que es igual a
La derivación se ha completado y el esquema backpropagation neural network se resume en el siguiente algoritmo
Algorithm 4.2 (Algoritmo Backpropagation Gradiente descendente)
Inicialización
Inicializar todos los pesos y sesgos sinápticos al azar con valores pequeños, pero no muy pequeños.
Seleccione el tamaño del paso \(\mu\).
Fije \(y_{nj}^{1}=x_{nj},\quad j=1,2,\dots,k_{1}:=l,\quad n=1,2,\dots,N\)
Repeat Cada repetición completa un epoch
For \(n=1,2,\dots,N\) Do
For \(r=1,2,\dots,L\) Do Cálculo Forward
For \(j=1,2,\dots,k_{r}\) Do
Calcule \(z_{nj}^{r}\) a partir de la Ecuación (4.5) Calcule \(y_{nj}^{r}=f(z_{nj}^{r})\)
End For
End For
For \(j = 1, 2,\dots, k_{L}\), Do; Cálculo Backward (output layer)
Calcule \(\delta_{nj}^{L}\) a partir de la Ecuación (4.10)
End For
For \(r=L, L-1,\dots, 2\), Do; Cálculo Backward (hidden layers)
For \(j=1,2,\dots, k_{r}\), Do
Calcule \(\delta_{nj}^{r-1}\) a partir de la Ecuación (4.10)
End For
End For
End For
For \(r=1,2,\dots,L\), Do: Actualice los pesos
For \(j=1,2,\dots,k_{r}\), Do
Calcule \(\Delta\boldsymbol{\theta}_{j}^{r}\) a partir de la Ecuación (4.6)
\(\boldsymbol{\theta}_{j}^{r}=\boldsymbol{\theta}_{j}^{r}+\Delta\boldsymbol{\theta}_{j}^{r}\)
End For
End For
Until Un criterio de parada se cumpla.
El algoritmo de backpropagation puede reivindicar una serie de padres. La popularización del algoritmo se asocia con el artículo clásico [Rumelhart et al., 1986], donde se proporciona la derivación del algoritmo. La idea de backpropagation también aparece en [Bryson Jr et al., 1963] en el contexto del control óptimo.
Existen diferentes variaciones del algoritmo backpropagation, tales como: Gradiende descendente con término de momento, Algoritmo de momentos de Nesterov’s, Algoritmo AdaGrad, RMSProp con momento de Nesterov, Algortimo de estimación de momentos adaptativo los cuales pueden ser utlizados para resolver la tarea de optimización (ver [Theodoridis, 2020]).
4.9. Las capas ocultas#
¿Cuántas capas ocultas?. Si los datos son linealmente separables (lo que a menudo se sabe cuando se empieza a codificar una ANN, SVM puede servir de test), entonces no se necesita ninguna capa oculta. Por supuesto, tampoco se necesita una ANN para resolver los datos, pero está seguirá haciendo su trabajo.
Sobre la configuración de las capas ocultas en las ANNs, existe un consenso dentro de este tema, y es la diferencia de rendimiento al añadir capas ocultas adicionales: las situaciones en las que el rendimiento mejora con una segunda (o tercera, etc.) capa oculta son muy pocas. Una capa oculta es suficiente para la gran mayoría de los problemas.
Entonces, ¿qué pasa con el tamaño de la(s) capa(s) oculta(s), cuántas neuronas?. Existen algunas reglas empíricas; de ellas, la más utilizada es
'The optimal size of the hidden layer is usually between the size of the input and size of the output layers'.Jeff Heaton, the author of Introduction to Neural Networks in Java.
Hay una regla empírica adicional que ayuda en los problemas de aprendizaje supervisado. Normalmente se puede evitar el sobreajuste si se mantiene el número de neuronas por debajo de:
\[ N_{h}=\frac{N_{s}}{(\alpha\cdot(N_{i}+N_{o}))} \]\(N_{i}=\) número de neuronas de entrada
\(N_{o}=\) número de neuronas de salida
\(N_{s}=\) número de muestras en el conjunto de datos de entrenamiento
\(\alpha=\) un factor de escala arbitrario, normalmente 2-10
Un valor de \(\alpha=2\) suele funcionar sin sobreajustar. Se puede pensar en \(\alpha\) como el factor de ramificación efectivo o el número de pesos distintos de cero para cada neurona. Las capas de salida harán que el factor de ramificación “efectivo” sea muy inferior al factor de ramificación medio real de la red. Para profundizar mas en el diseño de redes neuronales, ver el siguiente texto de Martin Hagan.
En resumen, para la mayoría de los problemas, probablemente se podría obtener un rendimiento decente (incluso sin un segundo paso de optimización) estableciendo la configuración de la capa oculta utilizando sólo dos reglas:
el número de capas ocultas es igual a uno
el número de neuronas de esa capa es la media de las neuronas de las capas de entrada y salida.
4.10. Redes Neuronales Recurrentes#
Introducción
Nuestro interés en esta sección se centra en el caso de los datos secuenciales. Es decir, los vectores de entrada no son independientes, sino que aparecen en secuencia. Además, el orden específico en que se producen encierra información importante. Por ejemplo, este tipo de secuencias se dan en el
reconocimiento del habla y en el procesamiento del lenguaje, como la traducción automática, así como también el pronóstico de series de tiempo financieras. Sin duda, la secuencia en la que se producen las palabras es de suma importancia.
Redes Neuronales Recurrentes
Las variables que intervienen en una RNN son:
Vector de estado en el tiempo \(n\), denotado como \(\boldsymbol{h}_{n}\). El símbolo nos recuerda que \(\boldsymbol{h}\) es un vector de variables ocultas (
hidden layer); el vector de estado constituye la memoria del sistema,Vector de entrada en el momento \(n\), denominado \(\boldsymbol{x}_{n}\),
Vector de salida en el momento \(n\), \(\hat{\boldsymbol{y}}_{n}\), y el vector de salida objetivo, \(\boldsymbol{y}_{n}\).
El modelo se describe mediante un conjunto de matrices y vectores de parámetros desconocidos, a saber, \(U, W, V , \boldsymbol{b}\) y \(\boldsymbol{c}\), que deben aprenderse durante el entrenamiento.
Las ecuaciones que describen un modelo RNN son
(4.11)#\[\begin{split} \begin{align*} \boldsymbol{h}_{n}&=f(U\boldsymbol{x}_{n}+W\boldsymbol{h}_{n-1}+\boldsymbol{b})\\ \hat{\boldsymbol{y}}_{n}&=g(V\boldsymbol{h}_{n}+\boldsymbol{c}). \end{align*} \end{split}\]donde las funciones no lineales \(f\) y \(g\) actúan elemento a elemento (element-wise) y se aplican individualmente a cada elemento de sus argumentos vectoriales.
En otras palabras, una vez que se ha observado un nuevo vector de entrada, se actualiza el vector de estado. Su nuevo valor depende de la información más reciente, transmitida por la entrada \(\boldsymbol{x}_{n}\) así como de la historia pasada, ya que esta se ha acumulado en \(\boldsymbol{h}_{n-1}\). La salida depende del vector de estado actualizado, \(\boldsymbol{h}_{n}\). Es decir, depende de la “historia” hasta el instante actual \(n\), tal y como se expresa en \(\boldsymbol{h}_{n}\).
Las opciones típicas para \(f\) son la tangente hiperbólica,
tanh, o las no linealidadesReLU. El valor inicial \(\boldsymbol{h}_{0}\) suele ser igual al vector cero. La no linealidad de salida, \(g\), se elige a menudo para ser la funciónsoftmax.
Fig. 4.14 Arquitectura de una Red Neuronal Recurrente. Fuente [Theodoridis, 2020].#
De las ecuaciones anteriores se deduce que las matrices y vectores de parámetros se comparten en todos los instantes temporales. Durante el entrenamiento, se inicializan mediante números aleatorios. El modelo asociado con la Ecuación (4.11) se muestra en la Fig. 4.14A.
En la Fig. 4.14B, el gráfico se despliega sobre los distintos instantes de tiempo para los que se dispone de observaciones. Por ejemplo, si la secuencia de interés es una frase de 10 palabras, entonces \(N\) se establece igual a 10, mientras que \(\boldsymbol{x}_{n}\) es el vector que codifica las respectivas palabras de entrada.
4.10.1. Backpropagation en tiempo#
Introducción
El entrenamiento de las RNN sigue una lógica similar a la del algoritmo backpropagation para el entrenamiento de redes neuronales de avance. Después de todo, una RNN puede verse como una red feed-forward con \(N\) capas. La capa superior es la del instante de tiempo \(N\) y la primera capa corresponde al instante de tiempo \(n = 1\). Una diferencia radica en que las capas ocultas en una RNN también producen salidas, es decir, \(\hat{\boldsymbol{y}}_{n}\), y se alimentan directamente con entradas. Sin embargo, en lo que respecta al entrenamiento, estas diferencias no afectan al razonamiento principal.
El aprendizaje de las matrices y vectores de parámetros desconocidos se consigue mediante un esquema de gradiente descendente, de acuerdo con Eq. (4.3). Resulta que los gradientes requeridos de la función de coste, con respecto a los parámetros desconocidos, tienen lugar recursivamente, comenzando en el último instante de tiempo, \(N\), y retrocediendo en el tiempo, \(n = N-1, N-2,\dots,1\). Esta es la razón por la que el algoritmo se conoce como bakpropagation a traves del tiempo (BPTT).
La función de coste es la suma a lo largo del tiempo, \(n\), de las correspondientes contribuciones a la función de pérdida, que dependen de los valores respectivos de \(\boldsymbol{h}_{n}, \boldsymbol{x}_{n}\), es decir,
Por ejemplo, tomar \(J_{n}\) como
MSEes la elección más común en series de tiempo conBPTTdebido a su simplicidad y compatibilidad con problemas de predicción continua. Sin embargo, si necesitas una mayor robustez ante outliers, podrías considerarMAEo Huber Loss\[ J_{n}(U, W, V, \boldsymbol{b}, \boldsymbol{c}):=\frac{1}{K}\sum_{k}(y_{nk}-\hat{y}_{nk})^{2}, \]donde la suma es sobre la dimensionalidad de \(\boldsymbol{y}\), y
\[ \hat{\boldsymbol{y}}_{n}=g(\boldsymbol{h}_{n}, V, \boldsymbol{c})~\text{y}~\boldsymbol{h}_{n}=f(\boldsymbol{x}_{n}, \boldsymbol{h}_{n-1}, U, W, \boldsymbol{b}). \]
En el corazón del cálculo de los gradientes de la función de coste con respecto a las diversas matrices y vectores de parámetros se encuentra el cálculo de los gradientes de \(J\) con respecto a los vectores de estado, \(\boldsymbol{h}_{n}\). Una vez calculados estos últimos, el resto de los gradientes, con respecto a las matrices y vectores de parámetros desconocidos, es una tarea sencilla. Para ello, nótese que cada \(h_{n},~n=1, 2,\dots,N-1\), afecta a \(J\) de dos maneras:
Directamente, a traves de \(J_{n}\)
Indirectamente, a través de la cadena que impone la estructura RNN, es decir,
\[ \boldsymbol{h}_{n}\rightarrow\boldsymbol{h}_{n+1}\rightarrow\cdots\rightarrow\boldsymbol{h}_{N}. \]Es decir, \(\boldsymbol{h}_{n}\), además de \(J_{n}\), también afecta a todos los valores de coste posteriores, \(J_{n+1},\dots, J_{N}\).
Nótese que, a través de la cadena, \(\boldsymbol{h}_{n+1}=f(\boldsymbol{x}_{n+1}, \boldsymbol{h}_{n}, U, W, \boldsymbol{b})\), empleando la regla de la cadena para las derivadas, las dependencias anteriores conducen al siguiente cálculo recursivo:
(4.12)#\[ \frac{\partial J}{\partial\boldsymbol{h}_{n}}=\underbrace{{\left(\frac{\partial\boldsymbol{h}_{n+1}}{\partial\boldsymbol{h}_{n}}\right)^{T}}\frac{\partial J}{\partial\boldsymbol{h}_{n+1}}}_{\text{parte recursiva indirecta}}+\underbrace{\left(\frac{\partial\hat{\boldsymbol{y}}_{n}}{\partial\boldsymbol{h}_{n}}\right)^{T}\frac{\partial J}{\partial\hat{\boldsymbol{y}}_{n}}}_{\text{parte directa}}, \]donde, por definición, la derivada de un vector, digamos, \(\boldsymbol{y}\), con respecto a otro vector, digamos, \(\boldsymbol{x}\), se define como la matriz
\[ \left[\frac{\partial\boldsymbol{y}}{\partial\boldsymbol{x}}\right]_{ij}:=\frac{\partial y_{i}}{\partial x_{j}}. \]
Nótese que el gradiente de la función de coste, con respecto a los parámetros ocultos (vector de estado) en la capa “\(n\)”, se da como una función del gradiente respectivo en la capa anterior, es decir, con respecto al vector de estado en el tiempo \(n+1\). Las dos pasadas requeridas por backpropagation en el tiempo se resumen a continuación.
Pasadas de Backpropagation en Tiempo
Paso hacia adelante:
Iniciando en \(n=1\) y utilizando las estimaciones actuales de las matrices y vectores de parámetros implicados, calcular en secuencia,
\[ (\boldsymbol{h}_{1}, \hat{\boldsymbol{y}}_{1})\rightarrow(\boldsymbol{h}_{2}, \hat{\boldsymbol{y}}_{2})\rightarrow\cdots\rightarrow(\boldsymbol{h}_{N}, \hat{\boldsymbol{y}}_{N}). \]Paso hacia atrás:
Empezando en \(n = N\), a calcular en secuencia,
\[ \frac{\partial J}{\partial\boldsymbol{h}_{N}}\rightarrow\frac{\partial J}{\partial\boldsymbol{h}_{N-1}}\rightarrow\cdots\rightarrow\frac{\partial J}{\partial\boldsymbol{h}_{1}}. \]
Nótese que el cálculo del gradiente \(\partial J/\partial\boldsymbol{h}_{N}\) es sencillo, y solo involucra la parte directa en Eq. (4.12).
Para la implementación de la
BPTT, se procede ainicializar aleatoriamente las matrices y vectores desconocidos implicados,
calcular todos los gradientes requeridos, siguiendo los pasos indicados anteriormente, y
realizar las actualizaciones según el esquema de gradiente descendente.
Los pasos (2) y (3) se realizan de forma iterativa hasta que se cumple un criterio de convergencia, de forma análoga al algoritmo estándar de backpropagation.
4.10.2. Desvanecimiento y explosión de gradientes#
Introducción
La tarea de desvanecimiento y explosión de gradientes se ha introducido y discutido en secciones anteriores, en el contexto del algoritmo de backpropagation. Los mismos problemas se presentan en el algoritmo BPTT, dado que, este último es una forma específica del concepto de backpropagation y, como se ha dicho, una RNN puede considerarse como una red multicapa, donde cada instante de tiempo corresponde a una capa diferente. De hecho, en las RNN, el fenómeno de desvanecimiento/explosión de gradiente aparece de una forma bastante “agresiva”, teniendo en cuenta que \(N\) puede alcanzar valores grandes.
La naturaleza multiplicativa de la propagación de gradientes puede verse fácilmente en la Eq. (4.12). Para ayudar a comprender el concepto principal, simplifiquemos el escenario y supongamos que sólo interviene una variante de estado. Entonces los vectores de estado se convierten en escalares, \(h_{n}\) , y la matriz \(W\) en un escalar \(w\). Además, supongamos que las salidas también son escalares. Entonces la recursión en la Eq. (4.12) se simplifica como:
Suponiendo en la Eq. (4.11) que \(f\) es la función \(\tanh(\cdot)\) estándar, teniendo en cuenta que, \(\text{sech}^2(\cdot)=1-\tanh^2(\cdot)\) y \(|\tanh(\cdot)|<1\), sobre su dominio, se ve fácilmente que
Escribiendo la recursión para dos pasos sucesivos, repitiendo el proceso en Eq. (4.13), por ejemplo, para \(\partial h_{n+2}/\partial h_{n+1}\) obtenemos que,
No es difícil ver que la multiplicación de los términos menores que uno puede llevar a valores de desvanecimiento, sobre todo si tenemos en cuenta que, en la práctica, las secuencias pueden ser bastante grandes, por ejemplo, \(N = 100\). Por lo tanto para instantes de tiempo cercanos a \(n = 1\), la contribución al gradiente del primer término del lado derecho en la Eq. (4.14) implicará un gran número de productos de números menores que uno en magnitud. Por otra parte, el valor de \(w\) estará contribuyendo en \(w^{n}\) potencia. Por lo tanto, si su valor es mayor que uno, puede conducir a valores explosivos de los gradientes respectivos.
En varios casos, se puede truncar el algoritmo de backpropagation a unos pocos pasos de tiempo. Otra forma, es sustituir la no linealidad
tanhporReLU. Para el caso del valor explosivo, se puede introducir una técnica que recorte los valores a un umbral predeterminado, una vez que los valores superen ese umbral. Sin embargo, otra técnica que suele emplearse en la práctica es sustituir la formulación RNN estándar descrita anteriormente por una estructura alternativa, que puede hacer mejor frente a estos fenómenos causados por la dependencia a largo plazo de la RNN.
Observation 4.1
RNN profundas: Además de la red RNN básica que comprende una sola capa de estados, se han propuesto extensiones que implican múltiples capas de estados, una encima de otra (ver [Pascanu et al., 2013]).RNN bidireccionales: Como su nombre indica, en las RNN bidireccionales hay dos variables de estado, es decir, una denotada como \(\overset{\rightarrow}{\boldsymbol{h}}\), que se propaga hacia adelante, y otra, \(\overset{\leftarrow}{\boldsymbol{h}}\), que se propaga hacia atrás. De este modo, las salidas dependen tanto del pasado como del futuro (ver [Graves et al., 2013]).
4.10.3. Implementación#
import numpy as np
import matplotlib.pyplot as plt
Funciones de activación y su derivada (
sigmoide)
def sigmoid(x):
return 1 / (1 + np.exp(-x))
def sigmoid_derivative(x):
return x * (1 - x)
Función de activación de salida (
lineal)
def linear(x):
return x
def linear_derivative(x):
return 1
Clase para el modelo
RNNconBPTT
class RNN_BPTT:
def __init__(self, input_size, hidden_size, output_size, learning_rate=0.01, truncation=5):
# Inicialización de los pesos
self.U = np.random.randn(input_size, hidden_size) * 0.1
self.W = np.random.randn(hidden_size, hidden_size) * 0.1
self.V = np.random.randn(hidden_size, output_size) * 0.1
# Inicialización de los sesgos
self.b = np.zeros((1, hidden_size))
self.c = np.zeros((1, output_size))
self.learning_rate = learning_rate
self.truncation = truncation # Truncación de BPTT para evitar largos horizontes de tiempo
# Forward pass
def forward(self, X):
T = len(X)
self.h = np.zeros((T + 1, self.W.shape[0])) # Inicializamos el estado oculto
self.y_hat = np.zeros((T, self.V.shape[1])) # Salida predicha
for t in range(T):
# Estado oculto en el tiempo t
self.h[t] = sigmoid(np.dot(X[t], self.U) + np.dot(self.h[t-1], self.W) + self.b)
# Salida predicha en el tiempo t
self.y_hat[t] = linear(np.dot(self.h[t], self.V) + self.c)
return self.y_hat
# Backward pass (BPTT)
def backward(self, X, y):
T = len(X)
# Gradientes inicializados a cero
dU = np.zeros_like(self.U)
dW = np.zeros_like(self.W)
dV = np.zeros_like(self.V)
db = np.zeros_like(self.b)
dc = np.zeros_like(self.c)
# Inicialización del gradiente del estado oculto
dh_next = np.zeros_like(self.h[0])
for t in reversed(range(T)):
# Gradiente de la pérdida con respecto a la salida (derivada externa)
dy = self.y_hat[t] - y[t]
# Gradientes para V y c
dV += np.outer(self.h[t], dy)
dc += dy
# Gradiente del estado oculto
dh = np.dot(dy, self.V.T) + dh_next
dh_raw = dh * sigmoid_derivative(self.h[t])
# Gradientes para U, W y b
dU += np.outer(X[t], dh_raw)
dW += np.outer(self.h[t-1], dh_raw)
db += dh_raw
# Gradiente para la siguiente iteración (truncación de BPTT)
dh_next = np.dot(dh_raw, self.W.T)
# Actualizamos los pesos con el gradiente descendente
self.U -= self.learning_rate * dU
self.W -= self.learning_rate * dW
self.V -= self.learning_rate * dV
self.b -= self.learning_rate * db
self.c -= self.learning_rate * dc
# Entrenamiento
def train(self, X, y, epochs=1000):
for epoch in range(epochs):
output = self.forward(X)
self.backward(X, y)
if epoch % 100 == 0:
loss = np.mean(np.square(y - output)) # MSE
print(f"Epoch {epoch}, Loss: {loss}")
Dataset ejemplo (
senoidal)
def create_dataset(series, look_back=1):
X, y = [], []
for i in range(len(series) - look_back):
X.append(series[i:(i + look_back)])
y.append(series[i + look_back])
return np.array(X), np.array(y)
Función para graficar la serie de tiempo original y las predicciones usando
BPTT
def plot_predictions(time_series, predictions, look_back):
# Recortamos las predicciones para que coincidan con el tamaño de la serie original
predictions = predictions.flatten() # Convertir a 1D
prediction_shifted = np.zeros_like(time_series)
prediction_shifted[look_back:] = predictions # Desplazamos las predicciones para alinearlas con la serie
# Crear la figura
plt.plot(time_series, label="Serie de tiempo original", color='b')
plt.plot(prediction_shifted, label="Predicciones", color='r', linestyle='--')
# Añadir etiquetas y leyenda
plt.title("Comparación de la serie de tiempo y las predicciones")
plt.xlabel("Tiempo")
plt.ylabel("Valor")
plt.legend()
# Mostrar la figura
plt.show()
Ejemplo de uso
if __name__ == "__main__":
# Generamos una serie de tiempo simple
time_series = np.sin(np.linspace(0, 10, 100)) # Onda senoidal
look_back = 3
# Preparamos los datos
X, y = create_dataset(time_series, look_back)
X = X.reshape(-1, look_back)
y = y.reshape(-1, 1)
# Normalizamos los datos
X = (X - np.min(X)) / (np.max(X) - np.min(X))
y = (y - np.min(y)) / (np.max(y) - np.min(y))
# Definimos la red RNN
input_size = look_back
hidden_size = 10
output_size = 1
learning_rate = 0.01
rnn = RNN_BPTT(input_size, hidden_size, output_size, learning_rate)
rnn.train(X, y, epochs=1000)
predictions = rnn.forward(X)
plot_predictions(time_series, predictions, look_back)
Epoch 0, Loss: 0.38753522909014065
Epoch 100, Loss: 0.009228053210117604
Epoch 200, Loss: 0.0060878161714585256
Epoch 300, Loss: 0.005974296992324269
Epoch 400, Loss: 0.005867766210183003
Epoch 500, Loss: 0.005761767279040523
Epoch 600, Loss: 0.005656052406221255
Epoch 700, Loss: 0.005550395397161551
Epoch 800, Loss: 0.005444579304197365
Epoch 900, Loss: 0.005338393308295797
4.11. Red de memoria a largo plazo (LSTM)#
Observation 4.2
La idea clave de la red LSTM, propuesta en el artículo seminal [Hochreiter and Schmidhuber, 1997], es el llamado estado de celda, que ayuda a superar los problemas asociados con los fenómenos de desvanecimiento/explosión que son causados por las dependencias largo plazo dentro de la red. Las redes LSTM tienen la capacidad incorporada de controlar el flujo de información que entra y sale de la memoria del sistema mediante algoritmos no lineales.
Estas puertas se implementan mediante la no linealidad sigmoidea y un multiplicador. Desde un punto de vista algorítmico, las puertas equivalen a aplicar una ponderación al flujo de información correspondiente. Los pesos se sitúan en el intervalo \([0,1]\) y dependen de los valores de las variables implicadas que activan la no linealidad sigmoidea.
Fig. 4.15 Arquitectura de unidad LSTM. Fuente [Theodoridis, 2020].#
En otras palabras, la ponderación (control) de la información tiene lugar en el contexto. Según este razonamiento, la red tiene la agilidad de olvidar la información que ya ha sido utilizada y ya no es necesaria. La célula/unidad LSTM básica se se muestra en la Fig. 4.15. Se construye en torno a dos conjuntos de variables, apiladas en el vector \(\boldsymbol{s}\) (
cell statelong-term memory), que se conoce como el estado de la célula o unidad, y el vector \(\boldsymbol{h}\) (hidden stateshort-term memory), que se conoce como el vector de variables ocultas. Una red LSTM se construye a partir de la concatenación sucesiva de esta unidad básica. La unidad correspondiente al tiempo \(n\), además del vector de entrada, \(\boldsymbol{x}_{n}\), recibe \(\boldsymbol{s}_{n-1}\) y \(\boldsymbol{h}_{n-1}\) de la etapa anterior y pasa \(\boldsymbol{s}_{n}\) y \(\boldsymbol{h}_{n}\) a la siguiente.
A continuación se resumen las ecuaciones de actualización asociadas,
\[\begin{split} \begin{align*} \boldsymbol{f}&=\sigma(U^{f}\boldsymbol{x}_{n}+W^{f}\boldsymbol{h}_{n-1}+\boldsymbol{b}^{f}),\\ \boldsymbol{i}&=\sigma(U^{i}\boldsymbol{x}_{n}+W^{i}\boldsymbol{h}_{n-1}+\boldsymbol{b}^{i}),\\ \tilde{\boldsymbol{s}}&=\tanh(U^{s}\boldsymbol{x}_{n}+W^{s}\boldsymbol{h}_{n-1}+\boldsymbol{b}^{s}),\\ \boldsymbol{s}_{n}&=\boldsymbol{s}_{n-1}\circ f+\boldsymbol{i}\circ\tilde{\boldsymbol{s}},\\ \boldsymbol{o}&=\sigma(U^{o}\boldsymbol{x}_{n}+W^{o}\boldsymbol{h}_{n-1}+\boldsymbol{b}^{o}),\\ \boldsymbol{h}_{n}&=\boldsymbol{o}\circ\tanh(\boldsymbol{s}_{n}), \end{align*} \end{split}\]donde \(\circ\) denota el producto elemento a elemento entre vectores o matrices (producto de Hadamard), es decir, \((s\circ f)_{i} = s_{i}f_{i}\), y \(\sigma\) denota la función sigmoidea logística.
Puerta de olvido(\(\boldsymbol{f}\))Esta puerta decide qué información se debe olvidar de la celda. Se calcula de la siguiente manera:
\[ \boldsymbol{f} = \sigma(U^{f} \boldsymbol{x}_{n} + W^{f} \boldsymbol{h}_{n-1} + \boldsymbol{b}^{f}) \]\( U^{f} \): pesos de la puerta de olvido aplicados a la entrada actual \( \boldsymbol{x}_{n} \).
\( W^{f} \): pesos de la puerta de olvido aplicados al estado oculto anterior \( \boldsymbol{h}_{n-1} \).
\( \boldsymbol{b}^{f} \): sesgo de la puerta de olvido.
\( \sigma \): función sigmoide.
Puerta de entrada(\(\boldsymbol{i}\))Esta puerta decide qué nueva información se almacenará en la celda:
\[ \boldsymbol{i} = \sigma(U^{i} \boldsymbol{x}_{n} + W^{i} \boldsymbol{h}_{n-1} + \boldsymbol{b}^{i}) \]\( U^{i} \): pesos de la puerta de entrada aplicados a la entrada actual \( \boldsymbol{x}_{n} \).
\( W^{i} \): pesos de la puerta de entrada aplicados al estado oculto anterior \( \boldsymbol{h}_{n-1} \).
\( \boldsymbol{b}^{i} \): sesgo de la puerta de entrada.
Candidato de estado(\( \tilde{\boldsymbol{s}} \))Genera nuevas candidaturas para el estado de la celda:
\[ \tilde{\boldsymbol{s}} = \tanh(U^{s} \boldsymbol{x}_{n} + W^{s} \boldsymbol{h}_{n-1} + \boldsymbol{b}^{s}) \]\( U^{s} \): pesos del candidato de estado aplicados a la entrada actual \( \boldsymbol{x}_{n} \).
\( W^{s} \): pesos del candidato de estado aplicados al estado oculto anterior \( \boldsymbol{h}_{n-1} \).
\( \boldsymbol{b}^{s} \): sesgo del candidato de estado.
\( \tanh \): función tangente hiperbólica.
Estado de la celda(\( \boldsymbol{s}_{n} \))Actualiza el estado de la celda utilizando la puerta de olvido y la puerta de entrada:
\[ \boldsymbol{s}_{n} = \boldsymbol{s}_{n-1} \circ \boldsymbol{f} + \boldsymbol{i} \circ \tilde{\boldsymbol{s}} \]\( \circ \): multiplicación elemento a elemento.
\( \boldsymbol{s}_{n-1} \): estado de la celda en el paso temporal anterior.
Puerta de salida(\( \boldsymbol{o} \))Decide qué parte del estado de la celda se enviará a la salida:
\[ \boldsymbol{o} = \sigma(U^{o} \boldsymbol{x}_{n} + W^{o} \boldsymbol{h}_{n-1} + \boldsymbol{b}^{o}) \]\( U^{o} \): pesos de la puerta de salida aplicados a la entrada actual \( \boldsymbol{x}_{n} \).
\( W^{o} \): pesos de la puerta de salida aplicados al estado oculto anterior \( \boldsymbol{h}_{n-1} \).
\( \boldsymbol{b}^{o} \): sesgo de la puerta de salida.
Estado oculto(\(\boldsymbol{h}_{n}\))Finalmente, se calcula el estado oculto, que es la salida de la celda LSTM:
\[ \boldsymbol{h}_{n} = \boldsymbol{o} \circ \tanh(\boldsymbol{s}_{n}) \]Resumen de salida de una redLSTMEstado de la Celda (\( \boldsymbol{s}_{n} \)): almacena información a largo plazo, afectada por las puertas de olvido y entrada.
Estado Oculto (\( \boldsymbol{h}_{n} \)): representa la salida de la celda LSTM en el tiempo \( n \), que se utiliza como entrada para el siguiente paso temporal y puede usarse como salida de la red.
Observation 4.3
Nótese que el estado de la célula, \(\boldsymbol{s}\), pasa información directa del instante anterior al siguiente. Esta información es controlada primero por la primera puerta, según los elementos en \(f\), que toman valores en el rango \([0, 1]\), dependiendo de la entrada actual y de las variables ocultas que se reciben de la etapa anterior. Esto es lo que decíamos antes, es decir, que la ponderación se ajusta en “contexto”. A continuación, se añade nueva información, es decir, \(\tilde{\boldsymbol{s}}\), a \(\boldsymbol{s}_{n-1}\), que también está controlada por la segunda red de puertas sigmoidales (es decir, \(\boldsymbol{i}\)). Así se garantiza que la información del pasado se transmita al futuro de forma directa, lo cual, ayuda a la red a memorizar información..
Resulta que este tipo de memoria explota mejor las dependencias de largo alcance en los datos, en comparación con la estructura RNN básica. El vector de variables ocultas \(\boldsymbol{h}\) está controlado tanto por el estado de la célula como por los valores actuales de las variables de entrada y de estados anteriores. Todas las matrices y vectores implicados se aprenden en la fase de entrenamiento. Nótese que hay dos líneas asociadas a \(\boldsymbol{h}_{n}\). La de la derecha conduce a la siguiente etapa y la de la parte superior se utiliza para proporcionar la salida, \(\hat{\boldsymbol{y}}_{n}\), en el tiempo \(n\), a través de, digamos, la no linealidad softmax, como en las RNN estándar en Eq. (4.11).
Variantes y Aplicaciones
Además de la estructura LSTM ya comentada, se han propuesto diversas variantes. Un extenso estudio comparativo entre diferentes arquitecturas LSTM y RNN se puede encontrar en [Greff et al., 2016, Jozefowicz et al., 2015]. Las RNNs y las LSTMs se han utilizado con éxito en una amplia gama de aplicaciones, tales como:
el
modelado del lenguaje[Sutskever et al., 2011],traducción de máquinas[Liu et al., 2014],reconocimiento del habla[Graves and Jaitly, 2014],visión artificialpara la generación de descriptores de imágenes [Karpathy and Fei-Fei, 2015],análisis de datos de fMRI para comprender la
dinámica temporal de las redes cerebralesasociadas [Seo et al., 2019].pronostico de
volatilidad de Bitcoin[Pratas et al., 2023]
Por ejemplo, en el procesamiento del lenguaje la entrada suele ser una secuencia de palabras, que se codifican como números (son punteros al diccionario disponible). La salida es la secuencia de palabras que hay que predecir. Durante el entrenamiento, se establece \(\boldsymbol{y}_{n} = \boldsymbol{x}_{n+1}\). Es decir, la red se entrena como un predictor no lineal.
4.12. Series de Tiempo#
Hasta ahora, en este curso hemos descrito métodos estadísticos tradicionales para el análisis de series temporales. En los capítulos anteriores, hemos analizado varios métodos para predecir la serie en un momento futuro a partir de observaciones tomadas en el pasado. Uno de estos métodos para realizar predicciones es el modelo autorregresivo (\(AR\)), que expresa la serie en el tiempo \(t\) como una regresión lineal de \(p\) observaciones anteriores
Aquí, \(\varepsilon_{t}\) es el término de error residual del modelo \(AR\). La idea que subyace al modelo lineal puede generalizarse en el sentido de que el objetivo de la predicción de series temporales es desarrollar una función \(f\) que prediga \(y_{t}\) en función de las observaciones en \(p\) valores previos en el tiempo
En este capítulo, exploraremos tres métodos basados en redes neuronales para desarrollar la función \(f\). Cada método incluye la definición de una arquitectura de red neuronal (en términos del número de capas ocultas, número de neuronas en cada capa oculta, etc.) y luego el algoritmo backpropagation o su variante adecuada para la arquitectura de red utilizada.
Los ejemplos de este capítulo serán implementados utilizando la
API Keraspara el aprendizaje profundo.Kerases unaAPIde alto nivel que permite definir diferentes arquitecturas de redes neuronales y entrenarlas utilizando varios optimizadores basados en gradientes. En el backend, Keras utiliza un marco computacional de bajo nivel implementado enC, C++yFORTRAN. Varios de estos entornos de bajo nivel están disponibles en código abierto.Kerassoporta los tres siguientes entornos: TensorFlow, que fue desarrollado por Google y es el backend por defecto de Keras, CNTK, un marco de código abierto de Microsoft, y Theano, desarrollado originalmente en la Universidad de Montreal, Canadá. Los ejemplos de este libro utilizanTensorFlowcomo backend. Por lo tanto, para ejecutar los ejemplos, necesitarás tantoKerascomoTensorFlowinstalados.
4.13. Perceptrones multicapa#
Los perceptrones multicapa (MLP) son las formas más básicas de las redes neuronales. Un
MLPconsta de tres componentes: una capa de entrada, varias capas ocultas y una capa de salida. Una capa de entrada representa un vector de regresores o características de entrada, por ejemplo, observaciones a partir de \(p\) puntos previos en el tiempo \([y_{t-1}, y_{t-2}, \dots, y_{t-p}]\). Las características de entrada se introducen en una capa oculta que tiene \(n\) neuronas, cada una de las cuales aplica una transformación lineal y una activación no lineal a las características de entrada.
La salida de una neurona es \(g_{i} = h(\boldsymbol{w}_{i}x + b_{i})\), donde \(\boldsymbol{w}_{i}\) y \(b_{i}\) son los pesos y el sesgo de la transformación lineal, y \(h\) es una función de activación no lineal. La función de activación no lineal permite a la red neuronal modelar no linealidades complejas de las relaciones subyacentes entre los regresores y la variable objetivo. Popularmente, \(h\) es la función sigmoid
\[ h(z)=\displaystyle{\frac{1}{1-e^{-z}}}, \]que mapea cualquier número real al intervalo \([0, 1]\).
Debido a esta propiedad, la función sigmoidea se utiliza para generar probabilidades de clase binarias y, por tanto, es comunmente usada en los modelos de clasificación. Otra opción de función de activación no lineal es la función
tanh\[ h(z)=\displaystyle{\frac{1-e^{-z}}{1+e^{-z}}}, \]la cual mapea cualquier número real al intervalo \([-1, 1]\). En algunos casos la función \(h\) es una función lineal o la identidad.
En el caso de una red neuronal de una sola capa oculta, como se muestra en Fig. 4.16, la salida de cada neurona se pasa a la capa de salida, que aplica una transformación lineal y función de activación para generar la predicción de la variable objetivo, que, en el caso de la predicción de series temporales, es el valor predicho de la serie en el tiempo \(t\)
En un
MLP(ver Fig. 4.17), se agrupan varias capas ocultas. La salida de las neuronas de una capa oculta se introduce en la siguiente capa oculta. Las neuronas de esta capa transforman la entrada y la pasan a la siguiente capa oculta. Finalmente, la última capa oculta alimenta la capa de salida
Las capas ocultas del
MLPtambién se denominan capas densas o, a veces, capas totalmente conectadas. El nombre denso tiene su significado en el hecho de que todas las neuronas de una capa densa están conectadas a todas las neuronas de la capa anterior y de la siguiente. Si la capa anterior es la capa de entrada, todas las características de entrada alimentan a cada una de las neuronas de la capa oculta.Debido a las conexiones múltiples entre la capa de entrada y la primera capa densa y entre las capas densas intermedias, un MLP tiene un número enorme de pesos entrenables. Por ejemplo, si el número de características de entrada es \(p\) y hay tres capas densas que tienen número de neuronas \(n_{1}, n_{2}\) y \(n_{3}\) respectivamente, entonces el número de pesos entrenables es
El último elemento de este cálculo es el número de pesos que conectan la tercera capa oculta y la capa de salida. Los MLP profundos tienen varias capas densas y cientos, incluso miles, de neuronas en cada capa. Por lo tanto, el número de pesos entrenables en los MLP profundos es muy grande.
4.14. Entrenamiento de MLP#
Los pesos \(w\) de una red neuronal se calculan mediante un algoritmo de optimización basado en el gradiente, como el gradiente decendiente estocástico, que minimiza iterativamente la función de pérdida o el error (\(L\)) en que incurre la red al hacer predicciones sobre los datos de entrenamiento. El error cuadrático medio (\(MSE\)) y el error absoluto medio (\(MAE\)) se utilizan para tareas de regresión, mientras que la pérdida binaria y categórica logarítmica son funciones de pérdida habituales en los problemas de clasificación.
Para la predicción de series temporales, \(MSE\) y \(MAE\) serían aptos para entrenar los modelos neuronales. Los algoritmos gradiente descendente funcionan moviendo los pesos, en iteraciones \(i\), a lo largo de su trayectoria de gradiente. El gradiente es la derivada parcial de la función de pérdida \(L\) con respecto al peso. La regla de actualización más sencilla para cambiar un peso \(w\) requiere los valores de los pesos, la derivada parcial de \(L\) con respecto a los pesos, y una tasa de aprendizaje \(\alpha\) que controla la rapidez con la que el punto desciende a lo largo del gradiente
Esta regla básica de actualización tiene diversas variantes que influyen en la convergencia del algoritmo. Sin embargo, una entrada crucial para todos los algoritmos basados en el gradiente es la derivada parcial que debe calcularse para todos los pesos de la red. En redes neuronales profundas, algunas de las cuales tienen millones de pesos, el cálculo de la derivada puede ser una tarea computacional gigantesca. En esta dirección es exactamente donde entra en juego el famoso algoritmo backpropagation que resuelve este problema de forma eficiente.
Para entender el algoritmo backpropagation, primero hay que conocer los grafos computacionales y cómo se utilizan para realizar cálculos en una red neuronal. Consideremos una red neuronal simple de una sola capa oculta, que tiene dos unidades ocultas, cada una con una activación sigmoidea. La unidad de salida es una transformación lineal de sus entradas. La red se alimenta con dos variables de entrada, \([y_{1},y_{2}]\). Los pesos se muestran a lo largo de los bordes de la red
La red realiza una serie de sumas, multiplicaciones y un par evaluaciones de funciones sigmoidales para transformar la entrada en una predicción. La transformación de la entrada en una predicción se denomina paso hacia delante de la red neuronal. La Fig. 4.19 muestra cómo se consigue un pase hacia adelante mediante un grafo computacional para un par de entrada \([-1, 2]\). Cada cálculo da como resultado una salida intermedia \(p_{i}\). Los resultados intermedios \(p_{7}\) y \(p_{8}\) son la salida de las neuronas ocultas \(g_{1}\) y \(g_{2}\). Durante el entrenamiento, la función de pérdida \(L\) también se calcula en el paso hacia delante
En este punto, se aplica el algoritmo de backpropagation para calcular las derivadas parciales entre dos nodos conectados por una arista. El recorrido hacia atrás en el grafo para calcular la derivada parcial también se conoce como paso hacia atrás (backward pass). El operador de diferenciación parcial se aplica en cada nodo y las derivadas parciales se asignan a las respectivas aristas que conectan el nodo descendente a lo largo del grafo computacional.
Siguiendo la regla de la cadena la derivada parcial \(\partial_{\omega} L\) se calcula multiplicando las derivadas parciales en todas las aristas que conectan el nodo de peso y el nodo de pérdida. Si existen varios caminos entre un nodo de peso y el nodo de pérdida, las derivadas parciales a lo largo de cada camino se suman para obtener la derivada parcial total de la pérdida con respecto al peso. Esta técnica gráfica de implementar los pasos hacia delante y hacia atrás es la técnica computacional subyacente utilizada en potentes bibliotecas de aprendizaje profundo. El paso hacia atrás se ilustra en Fig. 4.20
Las derivadas parciales de la función de pérdida con respecto a los pesos se obtienen aplicando la regla de la cadena:
Durante el entrenamiento, los pesos se inicializan con números aleatorios comúnmente muestreados a partir de una distribución uniforme con límites superior e inferior en \([-1, 1]\) o una distribución normal que tiene media cero y varianza unitaria. Estos esquemas de inicialización aleatoria tienen algunas variantes que mejoran la convergencia de la optimización. En este caso, vamos a suponer que los pesos son inicializados a partir de una distribución aleatoria uniforme y, por tanto, \(w_{1} = -0.33, w_{2} = -0.33, w_{3} = 0.57, w_{4} = -0.01, w_{5}=0.07\), y \(w_{6} = 0.82\).
Con estos valores, vamos a realizar pasos hacia adelante y hacia atrás sobre el grafo computacional. Actualizamos la figura anterior con los valores calculados durante la pasada hacia adelante en azul y los gradientes calculados durante la pasada hacia atrás en rojo. Para este ejemplo, fijamos el valor real de la variable objetivo como \(y = 1\)
Una vez calculados los gradientes a lo largo de las aristas, las derivadas parciales con respecto a los pesos no son más que una aplicación de la regla de la cadena, de la que ya hemos hablado. Los valores de las derivadas parciales son los siguientes:
El siguiente paso consiste en actualizar los pesos mediante el algoritmo de gradiente descendente. Así, con una tasa de aprendizaje de \(\alpha = 0.01\), el nuevo valor de
El resto de pesos también pueden actualizarse utilizando una regla de actualización similar. El proceso de actualización iterativa de los pesos se repite varias veces. El número de veces que se actualizan los pesos se conoce como número de épocas o pasadas sobre los datos de entrenamiento. Normalmente, un criterio de tolerancia sobre el cambio de la función de pérdida en comparación con la época anterior controla el número de épocas.
Para determinar los pesos de una red neuronal se utiliza el algoritmo backpropagation junto con un optimizador basado en el gradiente. Afortunadamente, existen potentes bibliotecas de aprendizaje profundo, como
Tensorflow, TheanoyCNTKque implementan gráficos computacionales para entrenar redes neuronales de cualquier arquitectura y complejidad. Estas bibliotecas vienen con soporte para ejecutar los cálculos como operaciones matemáticas en matrices multidimensionales y también pueden aprovechar lasGPUpara realizar cálculos más rápidos.
4.15. MLP para la predicción de series temporales#
En esta sección, utilizaremos
MLPpara desarrollar modelos de predicción de series temporales. El conjunto de datos utilizado para estos ejemplos es sobre la contaminación atmosférica medida por la concentración de partículas (PM) de diámetro inferior o igual a 2.5 micrómetros.Hay otras variables, como la presión atmosférica, temperatura del aire, punto de rocío, etc. Se han desarrollado un par de modelos de series temporales, uno sobre la presión atmosférica
PRESy otro sobrepm 2.5. El conjunto de datos se ha descargado del repositorio de aprendizaje automático de la UCI.
import sys
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
import datetime
from sklearn.preprocessing import MinMaxScaler
import tensorflow as tf
from tensorflow import keras
from keras.layers import Dense, Input, Dropout
from keras.optimizers import SGD
from keras.models import Model
from keras.models import load_model
from keras.callbacks import ModelCheckpoint
import os
from sklearn.metrics import r2_score
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings("ignore")
2025-02-07 11:56:31.729300: I tensorflow/core/util/port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
2025-02-07 11:56:31.938400: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:485] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
2025-02-07 11:56:32.017742: E external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:8454] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
2025-02-07 11:56:32.036661: E external/local_xla/xla/stream_executor/cuda/cuda_blas.cc:1452] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
2025-02-07 11:56:32.194974: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F AVX512_VNNI AVX512_BF16 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2025-02-07 11:56:32.960410: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT
devices = tf.config.list_physical_devices()
print("Available devices:")
for device in devices:
print(device)
gpus = tf.config.list_physical_devices('GPU')
if gpus:
print("TensorFlow is using GPU.")
for gpu in gpus:
gpu_details = tf.config.experimental.get_device_details(gpu)
print(f"GPU details: {gpu_details}")
else:
print("TensorFlow is not using GPU.")
Available devices:
PhysicalDevice(name='/physical_device:CPU:0', device_type='CPU')
PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')
TensorFlow is using GPU.
GPU details: {'compute_capability': (8, 9), 'device_name': 'NVIDIA GeForce RTX 4070'}
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1738947393.787247 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947394.071261 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947394.071301 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947394.072500 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
df = pd.read_csv('datasets/PRSA_data_2010.1.1-2014.12.31.csv', usecols=lambda column: column not in ['DEWP', 'TEMP', 'cbwd', 'Iws', 'Is', 'Ir'])
print('Shape of the dataframe:', df.shape)
Shape of the dataframe: (43824, 7)
df.head()
| No | year | month | day | hour | pm2.5 | PRES | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 2010 | 1 | 1 | 0 | NaN | 1021.0 |
| 1 | 2 | 2010 | 1 | 1 | 1 | NaN | 1020.0 |
| 2 | 3 | 2010 | 1 | 1 | 2 | NaN | 1019.0 |
| 3 | 4 | 2010 | 1 | 1 | 3 | NaN | 1019.0 |
| 4 | 5 | 2010 | 1 | 1 | 4 | NaN | 1018.0 |
Para asegurarse de que las filas están en el orden correcto de fecha y hora de las observaciones, se crea una nueva columna datetime a partir de las columnas relacionadas con la fecha y la hora del DataFrame. La nueva columna se compone de objetos
datetime.datetime de Python. El DataFrame se ordena en orden ascendente sobre esta columna
df['datetime'] = df[['year', 'month', 'day', 'hour']].apply(lambda row: datetime.datetime(year=row['year'], month=row['month'], day=row['day'],
hour=row['hour']), axis=1)
df.sort_values('datetime', ascending=True, inplace=True)
df.head()
| No | year | month | day | hour | pm2.5 | PRES | datetime | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 2010 | 1 | 1 | 0 | NaN | 1021.0 | 2010-01-01 00:00:00 |
| 1 | 2 | 2010 | 1 | 1 | 1 | NaN | 1020.0 | 2010-01-01 01:00:00 |
| 2 | 3 | 2010 | 1 | 1 | 2 | NaN | 1019.0 | 2010-01-01 02:00:00 |
| 3 | 4 | 2010 | 1 | 1 | 3 | NaN | 1019.0 | 2010-01-01 03:00:00 |
| 4 | 5 | 2010 | 1 | 1 | 4 | NaN | 1018.0 | 2010-01-01 04:00:00 |
Dibujamos un diagrama de cajas para visualizar la tendencia central y la dispersión de por ejemplo la columna
PRES
plot_fontsize = 18;
plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1)
g1 = sns.boxplot(x=df['PRES'], orient="h", palette="Set2")
g1.set_title('Box plot of Air Pressure', fontsize=plot_fontsize)
g1.set_xlabel('Air Pressure readings in hPa', fontsize=plot_fontsize)
plt.subplot(1, 2, 2)
g2 = sns.boxplot(x=df['pm2.5'], orient="h", palette="Set2")
g2.set_title('Box plot of PM2.5', fontsize=plot_fontsize)
g2.set_xlabel('PM2.5 readings in µg/m³', fontsize=plot_fontsize)
plt.show()
plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1)
g1 = sns.lineplot(data=df['PRES'])
g1.set_title('Time series of Air Pressure', fontsize=plot_fontsize)
g1.set_xlabel('Index', fontsize=plot_fontsize)
g1.set_ylabel('Air Pressure readings in hPa', fontsize=plot_fontsize)
plt.subplot(1, 2, 2)
g2 = sns.lineplot(data=df['pm2.5'])
g2.set_title('Time series of PM2.5', fontsize=plot_fontsize)
g2.set_xlabel('Index', fontsize=plot_fontsize)
g2.set_ylabel('PM2.5 readings in µg/m³', fontsize=plot_fontsize)
plt.show()
Observación
Antes de entrenar el modelo, el conjunto de datos se divide en dos partes: el conjunto de entrenamiento y el conjunto de validación. La red neuronal se entrena en el conjunto de entrenamiento. Esto significa que el cálculo de la función de pérdida, la propagación hacia atrás y los pesos actualizados mediante un algoritmo de gradiente descendente se realizan en el conjunto de entrenamiento.
El conjunto de validación se utiliza para evaluar el modelo y determinar el número de épocas en su entrenamiento. Aumentar el número de épocas reducirá aún más la función de pérdida en el conjunto de entrenamiento, pero no necesariamente tendrá el mismo efecto en el conjunto de validación debido al sobreajuste en el conjunto de entrenamiento, por lo que el número de épocas se controla manteniendo una evaluación y verificación sobre la función de pérdida calculada para el conjunto de validación.
Utilizamos
Kerascon el backendTensorflowpara definir y entrenar el modelo. Todos los pasos implicados en el entrenamiento y validación del modelo se realizan llamando a las funciones apropiadas de laAPIdeKeras.
Los tres primeros años, de 2010 a 2012, se utilizan como
entrenamiento, 2013 se usa como conjunto devalidacióny 2014 se utiliza como conjunto deprueba
split_date_test = datetime.datetime(year=2014, month=1, day=1, hour=0)
df_trainval = df.loc[df['datetime']<split_date_test]
df_test = df.loc[df['datetime']>=split_date_test]
print('Shape of trainval:', df_trainval.shape)
print('Shape of test:', df_test.shape)
Shape of trainval: (35064, 8)
Shape of test: (8760, 8)
Definimos los conjunto de entrenamiento y validación, necesarios para la construcción del modelo de Deep Learning
df_trainval.tail()
| No | year | month | day | hour | pm2.5 | PRES | datetime | |
|---|---|---|---|---|---|---|---|---|
| 35059 | 35060 | 2013 | 12 | 31 | 19 | 22.0 | 1013.0 | 2013-12-31 19:00:00 |
| 35060 | 35061 | 2013 | 12 | 31 | 20 | 18.0 | 1014.0 | 2013-12-31 20:00:00 |
| 35061 | 35062 | 2013 | 12 | 31 | 21 | 23.0 | 1014.0 | 2013-12-31 21:00:00 |
| 35062 | 35063 | 2013 | 12 | 31 | 22 | 20.0 | 1014.0 | 2013-12-31 22:00:00 |
| 35063 | 35064 | 2013 | 12 | 31 | 23 | 23.0 | 1014.0 | 2013-12-31 23:00:00 |
split_date_trainval = datetime.datetime(year=2013, month=1, day=1, hour=0)
df_train = df_trainval.loc[df_trainval['datetime']<split_date_trainval]
df_val = df_trainval.loc[df_trainval['datetime']>=split_date_trainval]
print('Shape of trainval:', df_train.shape)
print('Shape of test:', df_val.shape)
Shape of trainval: (26304, 8)
Shape of test: (8760, 8)
df_train.tail()
| No | year | month | day | hour | pm2.5 | PRES | datetime | |
|---|---|---|---|---|---|---|---|---|
| 26299 | 26300 | 2012 | 12 | 31 | 19 | 104.0 | 1017.0 | 2012-12-31 19:00:00 |
| 26300 | 26301 | 2012 | 12 | 31 | 20 | 131.0 | 1017.0 | 2012-12-31 20:00:00 |
| 26301 | 26302 | 2012 | 12 | 31 | 21 | 113.0 | 1018.0 | 2012-12-31 21:00:00 |
| 26302 | 26303 | 2012 | 12 | 31 | 22 | 45.0 | 1018.0 | 2012-12-31 22:00:00 |
| 26303 | 26304 | 2012 | 12 | 31 | 23 | 39.0 | 1018.0 | 2012-12-31 23:00:00 |
df_train.tail()
| No | year | month | day | hour | pm2.5 | PRES | datetime | |
|---|---|---|---|---|---|---|---|---|
| 26299 | 26300 | 2012 | 12 | 31 | 19 | 104.0 | 1017.0 | 2012-12-31 19:00:00 |
| 26300 | 26301 | 2012 | 12 | 31 | 20 | 131.0 | 1017.0 | 2012-12-31 20:00:00 |
| 26301 | 26302 | 2012 | 12 | 31 | 21 | 113.0 | 1018.0 | 2012-12-31 21:00:00 |
| 26302 | 26303 | 2012 | 12 | 31 | 22 | 45.0 | 1018.0 | 2012-12-31 22:00:00 |
| 26303 | 26304 | 2012 | 12 | 31 | 23 | 39.0 | 1018.0 | 2012-12-31 23:00:00 |
df_val.tail()
| No | year | month | day | hour | pm2.5 | PRES | datetime | |
|---|---|---|---|---|---|---|---|---|
| 35059 | 35060 | 2013 | 12 | 31 | 19 | 22.0 | 1013.0 | 2013-12-31 19:00:00 |
| 35060 | 35061 | 2013 | 12 | 31 | 20 | 18.0 | 1014.0 | 2013-12-31 20:00:00 |
| 35061 | 35062 | 2013 | 12 | 31 | 21 | 23.0 | 1014.0 | 2013-12-31 21:00:00 |
| 35062 | 35063 | 2013 | 12 | 31 | 22 | 20.0 | 1014.0 | 2013-12-31 22:00:00 |
| 35063 | 35064 | 2013 | 12 | 31 | 23 | 23.0 | 1014.0 | 2013-12-31 23:00:00 |
Los algoritmos de gradiente descendente funcionan mejor (por ejemplo, convergen más rápido) si las variables están dentro del intervalo \([-1, 1]\). Muchas fuentes relajan el límite hasta \([-3, 3]\). La variable
PRESes escalada conminmaxpara limitar la variable transformada dentro de \([0,1]\).
scaler = MinMaxScaler(feature_range=(0, 1))
df_train['scaled_PRES'] = scaler.fit_transform(np.array(df_train['PRES']).reshape(-1, 1))
scaler = MinMaxScaler(feature_range=(0, 1))
df_val['scaled_PRES'] = scaler.fit_transform(np.array(df_val['PRES']).reshape(-1, 1))
Restablecemos los índices del conjunto de validación
df_train.reset_index(drop=True, inplace=True)
df_train.head()
| No | year | month | day | hour | pm2.5 | PRES | datetime | scaled_PRES | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 2010 | 1 | 1 | 0 | NaN | 1021.0 | 2010-01-01 00:00:00 | 0.547170 |
| 1 | 2 | 2010 | 1 | 1 | 1 | NaN | 1020.0 | 2010-01-01 01:00:00 | 0.528302 |
| 2 | 3 | 2010 | 1 | 1 | 2 | NaN | 1019.0 | 2010-01-01 02:00:00 | 0.509434 |
| 3 | 4 | 2010 | 1 | 1 | 3 | NaN | 1019.0 | 2010-01-01 03:00:00 | 0.509434 |
| 4 | 5 | 2010 | 1 | 1 | 4 | NaN | 1018.0 | 2010-01-01 04:00:00 | 0.490566 |
df_val.reset_index(drop=True, inplace=True)
df_val.head()
| No | year | month | day | hour | pm2.5 | PRES | datetime | scaled_PRES | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 26305 | 2013 | 1 | 1 | 0 | 35.0 | 1018.0 | 2013-01-01 00:00:00 | 0.490909 |
| 1 | 26306 | 2013 | 1 | 1 | 1 | 31.0 | 1017.0 | 2013-01-01 01:00:00 | 0.472727 |
| 2 | 26307 | 2013 | 1 | 1 | 2 | 32.0 | 1017.0 | 2013-01-01 02:00:00 | 0.472727 |
| 3 | 26308 | 2013 | 1 | 1 | 3 | 21.0 | 1018.0 | 2013-01-01 03:00:00 | 0.490909 |
| 4 | 26309 | 2013 | 1 | 1 | 4 | 16.0 | 1018.0 | 2013-01-01 04:00:00 | 0.490909 |
También se grafican las series temporales de entrenamiento y validación normalizadas para PRES.
plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1)
g1 = sns.lineplot(df_train['scaled_PRES'], color='b')
g1.set_title('Time series of scaled Air Pressure in train set', fontsize=plot_fontsize)
g1.set_xlabel('Index', fontsize=plot_fontsize)
g1.set_ylabel('Scaled Air Pressure readings', fontsize=plot_fontsize);
plt.subplot(1, 2, 2)
g2 = sns.lineplot(df_val['scaled_PRES'], color='r')
g2.set_title('Time series of scaled Air Pressure in validation set', fontsize=plot_fontsize)
g2.set_xlabel('Index', fontsize=plot_fontsize)
g2.set_ylabel('Scaled Air Pressure readings', fontsize=plot_fontsize);
plt.show()
Regresores y Variable Objetivo
Ahora necesitamos generar los regresores (\(X\)) y la variable objetivo (\(y\)) para el entrenamiento y la validación. La matriz bidimensional de regresores y la matriz unidimensional objetivo se crean a partir de la matriz unidimensional original de la columna standardized_PRES en el DataFrame.
Para el modelo de predicción de series temporales, de este ejemplo, se utilizan las observaciones de los últimos siete días para predecir el día siguiente. Esto equivale a un modelo \(AR(7)\). Definimos una función que toma la serie temporal original y el número de pasos temporales en regresores como entrada para generar las matrices \(X\) e \(y\)
def makeXy(ts, nb_timesteps):
"""
Input:
ts: original time series
nb_timesteps: number of time steps in the regressors
Output:
X: 2-D array of regressors
y: 1-D array of target
"""
X = []
y = []
for i in range(nb_timesteps, ts.shape[0]):
if i-nb_timesteps <= 4:
print(i-nb_timesteps, i-1, i)
X.append(list(ts.loc[i-nb_timesteps:i-1])) #Regressors
y.append(ts.loc[i]) #Target
X, y = np.array(X), np.array(y)
return X, y
X_train, y_train = makeXy(df_train['scaled_PRES'], 7)
0 6 7
1 7 8
2 8 9
3 9 10
4 10 11
print('Shape of train arrays:', X_train.shape, y_train.shape)
Shape of train arrays: (26297, 7) (26297,)
X_val, y_val = makeXy(df_val['scaled_PRES'], 7)
0 6 7
1 7 8
2 8 9
3 9 10
4 10 11
print('Shape of validation arrays:', X_val.shape, y_val.shape)
Shape of validation arrays: (8753, 7) (8753,)
Ahora definimos la red
MLPutilizando laAPIfuncional deKeras. En este enfoque una capa puede ser declarada como la entrada de la siguiente capa en el momento de definir la siguiente
input_layer = Input(shape=(7,), dtype='float32')
En este caso,
Inputes una función que se utiliza para crear una capa de entrada en un modelo de red neuronal.shape=(7,)específica la forma de los datos de entrada. En este caso, significa que los datos de entrada tendrán 7 dimensiones.dtype='float32'especifica el tipo de datos de los elementos de la capa de entrada. En este caso, son números de punto flotante de 32 bits.
Las capas densas las definimos en esta caso con activación
tanh. Puede utilizar unGridSearchtal como se hizo en el curso de Machine Learning para encontrar hiperparámetros adecuados minimizando las métricas de regresión.
dense1 = Dense(32, activation='tanh')(input_layer)
dense2 = Dense(16, activation='tanh')(dense1)
dense3 = Dense(16, activation='tanh')(dense2)
I0000 00:00:1738947395.415149 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947395.415201 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947395.415220 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947395.527302 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
I0000 00:00:1738947395.527352 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
2025-02-07 11:56:35.527360: I tensorflow/core/common_runtime/gpu/gpu_device.cc:2112] Could not identify NUMA node of platform GPU id 0, defaulting to 0. Your kernel may not have been built with NUMA support.
I0000 00:00:1738947395.527393 9196 cuda_executor.cc:1001] could not open file to read NUMA node: /sys/bus/pci/devices/0000:01:00.0/numa_node
Your kernel may have been built without NUMA support.
2025-02-07 11:56:35.528327: I tensorflow/core/common_runtime/gpu/gpu_device.cc:2021] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 9558 MB memory: -> device: 0, name: NVIDIA GeForce RTX 4070, pci bus id: 0000:01:00.0, compute capability: 8.9
Dense e input_layer
Dense: Correspone a una capa totalmente conectada (fully connected).Unidades (Neurons): 32 neuronas.
dense1: Esta capa toma como entradainput_layer, que puede ser la capa de entrada del modelo u otra capa anterior.dense2: Esta capa toma como entrada la salida dedense1. Esto significa que los 32 valores de salida dedense1se usan como entrada paradense2. Similarmente, ocurre condense3
Observación
Las múltiples capas ocultas y el gran número de neuronas en cada capa oculta le dan a las redes neuronales la capacidad de modelar la compleja no linealidad de las relaciones subyacentes entre los regresores y el objetivo. Sin embargo, las redes neuronales profundas también pueden sobreajustar los datos de entrenamiento y dar malos resultados en el conjunto de validación o prueba. La función
Dropoutse ha utilizado eficazmente para regularizar las redes neuronales profundas.
En este ejemplo, se añade una capa
Dropoutantes de la capa de salida. Dropout aleatoriamente establece \(p\) fracción de neuronas de entrada a cero antes de pasar a la siguiente capa. La eliminación aleatoria de entradas actúa esencialmente como un tipo de ensamblaje de modelos de agregaciónbootstrap. Al apagar neuronas aleatoriamente,Dropoutestá forzando a la red a no depender excesivamente de ninguna unidad específica, mejorando la generalización.Por ejemplo, el bosque aleatorio utiliza el ensamblaje mediante la construcción de árboles en subconjuntos aleatorios de características de entrada. Utilizamos \(p=0.2\) para descartar el 20% de las características de entrada seleccionadas aleatoriamente.
dropout_layer = Dropout(0.2)(dense3)
Por último, la capa de salida predice la presión atmosférica del día siguiente
output_layer = Dense(1, activation='linear')(dropout_layer)
Las capas de entrada, densa y de salida se empaquetarán ahora dentro de un modelo, que es una clase envolvente para entrenar y hacer predicciones. Como función de pérdida se utiliza el error cuadrático medio (MSE). Los pesos de la red se optimizan mediante el algoritmo
Adam.Adamsignifica Estimación Adaptativa de Momentos y ha sido una opción popular para el entrenamiento de redes neuronales profundas.A diferencia del Gradiente descendente Estocástico, Adam utiliza diferentes tasas de aprendizaje para cada peso y las actualiza por separado a medida que avanza el entrenamiento. La tasa de aprendizaje de un peso se actualiza basándose en medias móviles ponderadas exponencialmente de los gradientes del peso y los gradientes al cuadrado.
ts_model = Model(inputs=input_layer, outputs=output_layer)
ts_model.compile(loss='mean_squared_error', optimizer='adam')
ts_model.summary()
Model: "functional"
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓ ┃ Layer (type) ┃ Output Shape ┃ Param # ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩ │ input_layer (InputLayer) │ (None, 7) │ 0 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense (Dense) │ (None, 32) │ 256 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_1 (Dense) │ (None, 16) │ 528 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_2 (Dense) │ (None, 16) │ 272 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dropout (Dropout) │ (None, 16) │ 0 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_3 (Dense) │ (None, 1) │ 17 │ └─────────────────────────────────┴────────────────────────┴───────────────┘
Total params: 1,073 (4.19 KB)
Trainable params: 1,073 (4.19 KB)
Non-trainable params: 0 (0.00 B)
Observación
En este caso,
Params #es calculado mediante la fórmula:Param # = #Entradas x #Neuronas + #Neuronas. Esto incluye los pesos de cada conexión entre las entradas y las neuronas y un bias para cada neurona. En este casoParams # = 7x32 + 32 = 256.El modelo se entrena llamando a la función
fit()en el objeto modelo y pasándoleX_trainyy_train. El entrenamiento se realiza para un número predefinido de épocas. Además,batch_sizedefine el número de muestras del conjunto de entrenamiento que se utilizarán para una instancia de backpropagation.
El conjunto de datos de validación también se pasa para evaluar el modelo después de cada
epochcompleta. Un objetoModelCheckpointrastrea la función de pérdida en el conjunto de validación y guarda el modelo para la época en la que la función de pérdida ha sido mínima.
save_weights_at = os.path.join('keras_models', 'PRSA_data_Air_Pressure_MLP_weights.{epoch:02d}-{val_loss:.4f}.keras')
save_best = ModelCheckpoint(save_weights_at, monitor='val_loss', verbose=0,
save_best_only=True, save_weights_only=False, mode='min', save_freq='epoch');
Aquí
val_losses el valor de la función de coste para los datos de validación cruzada ylosses el valor de la función de coste para los datos de entrenamiento. Converbose=0, no se imprime ningún mensaje en la consola durante el proceso de guardado del modelo.period=1indica que el modelo se evaluará y potencialmente se guardará después de cada época.
import os
from joblib import dump, load
history_airp = None
if os.path.exists('history_airp.joblib'):
history_airp = load('history_airp.joblib')
print("El archivo 'history_airp.joblib' ya existe. Se ha cargado el historial del entrenamiento.")
else:
history_airp = ts_model.fit(x=X_train, y=y_train, batch_size=16, epochs=20,
verbose=2, callbacks=[save_best], validation_data=(X_val, y_val),
shuffle=True);
dump(history_airp.history, 'history_airp.joblib')
print("El entrenamiento se ha completado y el historial ha sido guardado en 'history_airp.joblib'.")
El archivo 'history_airp.joblib' ya existe. Se ha cargado el historial del entrenamiento.
En este caso, el modo
verbose=2muestra una barra de progreso por cada época. Los modos posibles son:0para no mostrar nada,1para mostrar la barra de progreso,2para mostrar una línea por época. Las muestras se mezclan aleatoriamente antes de cada época (shuffle=True).
Se hacen predicciones para la presión atmosférica a partir del mejor modelo guardado. Las predicciones del modelo sobre la presión atmosférica escalada, se transforman inversamente para obtener predicciones sobre la presión atmosférica original. También se calcula la bondad de ajuste o
R cuadrado.
import os
import re
from tensorflow.keras.models import load_model
model_dir = 'keras_models'
files = os.listdir(model_dir)
pattern = r"PRSA_data_Air_Pressure_MLP_weights\.(\d+)-([\d\.]+)\.keras"
best_val_loss = float('inf')
best_model_file = None
best_model = None
for file in files:
match = re.match(pattern, file)
if match:
epoch = int(match.group(1))
val_loss = float(match.group(2))
if val_loss < best_val_loss:
best_val_loss = val_loss
best_model_file = file
if best_model_file:
best_model_path = os.path.join(model_dir, best_model_file)
print(f"Cargando el mejor modelo: {best_model_file} con val_loss: {best_val_loss}")
best_model = load_model(best_model_path)
else:
print("No se encontraron archivos de modelos que coincidan con el patrón.")
Cargando el mejor modelo: PRSA_data_Air_Pressure_MLP_weights.03-0.0002.keras con val_loss: 0.0002
preds = best_model.predict(X_val)
pred_PRES = scaler.inverse_transform(preds)
pred_PRES = np.squeeze(pred_PRES)
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1738947396.494609 9344 service.cc:146] XLA service 0x7f4528005f80 initialized for platform CUDA (this does not guarantee that XLA will be used). Devices:
I0000 00:00:1738947396.494632 9344 service.cc:154] StreamExecutor device (0): NVIDIA GeForce RTX 4070, Compute Capability 8.9
2025-02-07 11:56:36.504431: I tensorflow/compiler/mlir/tensorflow/utils/dump_mlir_util.cc:268] disabling MLIR crash reproducer, set env var `MLIR_CRASH_REPRODUCER_DIRECTORY` to enable.
2025-02-07 11:56:36.582925: I external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:531] Loaded cuDNN version 8907
1/274 ━━━━━━━━━━━━━━━━━━━━ 1:35 349ms/step
94/274 ━━━━━━━━━━━━━━━━━━━━ 0s 541us/step
195/274 ━━━━━━━━━━━━━━━━━━━━ 0s 520us/step
I0000 00:00:1738947396.799965 9344 device_compiler.h:188] Compiled cluster using XLA! This line is logged at most once for the lifetime of the process.
274/274 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
274/274 ━━━━━━━━━━━━━━━━━━━━ 1s 1ms/step
r2 = r2_score(df_val['PRES'].loc[7:], pred_PRES)
print('R-squared for the validation set:', round(r2,4))
R-squared for the validation set: 0.9942
plt.figure(figsize=(5.5, 5.5))
plt.plot(range(50), df_val['PRES'].loc[7:56], linestyle='-', marker='*', color='r')
plt.plot(range(50), pred_PRES[:50], linestyle='-', marker='.', color='b')
plt.legend(['Actual','Predicted'], loc=2)
plt.title('Actual vs Predicted Air Pressure')
plt.ylabel('Air Pressure')
plt.xlabel('Index');
Para predecir la variable
pm2.5usandoMLPusamos la implementación presentada a continuación
df = pd.read_csv('datasets/PRSA_data_2010.1.1-2014.12.31.csv', usecols=lambda column: column not in ['DEWP', 'TEMP', 'cbwd', 'Iws', 'Is', 'Ir'])
print('Shape of the dataframe:', df.shape)
Shape of the dataframe: (43824, 7)
df.head()
| No | year | month | day | hour | pm2.5 | PRES | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 2010 | 1 | 1 | 0 | NaN | 1021.0 |
| 1 | 2 | 2010 | 1 | 1 | 1 | NaN | 1020.0 |
| 2 | 3 | 2010 | 1 | 1 | 2 | NaN | 1019.0 |
| 3 | 4 | 2010 | 1 | 1 | 3 | NaN | 1019.0 |
| 4 | 5 | 2010 | 1 | 1 | 4 | NaN | 1018.0 |
Las filas con valores NaN en la columna pm2.5 se eliminan
df.dropna(subset=['pm2.5'], axis=0, inplace=True)
df.reset_index(drop=True, inplace=True)
df['datetime'] = df[['year', 'month', 'day', 'hour']].apply(lambda row: datetime.datetime(year=row['year'], month=row['month'], day=row['day'],
hour=row['hour']), axis=1)
df.sort_values('datetime', ascending=True, inplace=True)
plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1)
g1 = sns.boxplot(df['pm2.5'], orient="h", palette="Set2")
g1.set_title('Box plot of pm2.5', fontsize=plot_fontsize);
g1.set_xlabel(xlabel='pm2.5', fontsize=plot_fontsize)
g1.set(yticklabels=[])
g1.tick_params(left=False);
plt.subplot(1, 2, 2)
g2 = sns.lineplot(df['pm2.5'])
g2.set_title('Time series of pm2.5', fontsize=plot_fontsize)
g2.set_xlabel('Index', fontsize=plot_fontsize)
g2.set_ylabel('pm2.5 readings');
plt.show()
plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1)
g1 = sns.lineplot(df['pm2.5'].loc[df['datetime']<=datetime.datetime(year=2010,month=6,day=30)], color='g')
g1.set_title('pm2.5 during 2010', fontsize=plot_fontsize)
g1.set_xlabel('Index', fontsize=plot_fontsize)
g1.set_ylabel('pm2.5 readings', fontsize=plot_fontsize);
plt.subplot(1, 2, 2)
g = sns.lineplot(df['pm2.5'].loc[df['datetime']<=datetime.datetime(year=2010,month=1,day=31)], color='g')
g.set_title('Zoom in on one month: pm2.5 during Jan 2010', fontsize=plot_fontsize)
g.set_xlabel('Index', fontsize=plot_fontsize)
g.set_ylabel('pm2.5 readings', fontsize=plot_fontsize);
plt.show()
scaler = MinMaxScaler(feature_range=(0, 1))
df['scaled_pm2.5'] = scaler.fit_transform(np.array(df['pm2.5']).reshape(-1, 1))
split_date = datetime.datetime(year=2014, month=1, day=1, hour=0)
df_train = df.loc[df['datetime']<split_date]
df_val = df.loc[df['datetime']>=split_date]
print('Shape of train:', df_train.shape)
print('Shape of test:', df_val.shape)
Shape of train: (33096, 9)
Shape of test: (8661, 9)
df_train.head()
| No | year | month | day | hour | pm2.5 | PRES | datetime | scaled_pm2.5 | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 25 | 2010 | 1 | 2 | 0 | 129.0 | 1020.0 | 2010-01-02 00:00:00 | 0.129779 |
| 1 | 26 | 2010 | 1 | 2 | 1 | 148.0 | 1020.0 | 2010-01-02 01:00:00 | 0.148893 |
| 2 | 27 | 2010 | 1 | 2 | 2 | 159.0 | 1021.0 | 2010-01-02 02:00:00 | 0.159960 |
| 3 | 28 | 2010 | 1 | 2 | 3 | 181.0 | 1022.0 | 2010-01-02 03:00:00 | 0.182093 |
| 4 | 29 | 2010 | 1 | 2 | 4 | 138.0 | 1022.0 | 2010-01-02 04:00:00 | 0.138833 |
df_val.head()
| No | year | month | day | hour | pm2.5 | PRES | datetime | scaled_pm2.5 | |
|---|---|---|---|---|---|---|---|---|---|
| 33096 | 35065 | 2014 | 1 | 1 | 0 | 24.0 | 1014.0 | 2014-01-01 00:00:00 | 0.024145 |
| 33097 | 35066 | 2014 | 1 | 1 | 1 | 53.0 | 1013.0 | 2014-01-01 01:00:00 | 0.053320 |
| 33098 | 35067 | 2014 | 1 | 1 | 2 | 65.0 | 1013.0 | 2014-01-01 02:00:00 | 0.065392 |
| 33099 | 35068 | 2014 | 1 | 1 | 3 | 70.0 | 1013.0 | 2014-01-01 03:00:00 | 0.070423 |
| 33100 | 35069 | 2014 | 1 | 1 | 4 | 79.0 | 1012.0 | 2014-01-01 04:00:00 | 0.079477 |
df_val.reset_index(drop=True, inplace=True)
plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1)
g1 = sns.lineplot(df_train['scaled_pm2.5'], color='b')
g1.set_title('Time series of scaled pm2.5 in train set', fontsize=plot_fontsize)
g1.set_xlabel('Index', fontsize=plot_fontsize)
g1.set_ylabel('Scaled pm2.5 readings', fontsize=plot_fontsize);
plt.subplot(1, 2, 2)
g2 = sns.lineplot(df_val['scaled_pm2.5'], color='r')
g2.set_title('Time series of scaled pm2.5 in validation set', fontsize=plot_fontsize)
g2.set_xlabel('Index', fontsize=plot_fontsize)
g2.set_ylabel('Scaled pm2.5 readings', fontsize=plot_fontsize);
plt.show()
X_train, y_train = makeXy(df_train['scaled_pm2.5'], 7)
0 6 7
1 7 8
2 8 9
3 9 10
4 10 11
print('Shape of train arrays:', X_train.shape, y_train.shape)
Shape of train arrays: (33089, 7) (33089,)
X_val, y_val = makeXy(df_val['scaled_pm2.5'], 7)
0 6 7
1 7 8
2 8 9
3 9 10
4 10 11
print('Shape of validation arrays:', X_val.shape, y_val.shape)
Shape of validation arrays: (8654, 7) (8654,)
Definimos la capa de entrada que tiene forma
(None, 7)y de tipofloat32
input_layer = Input(shape=(7,), dtype='float32')
dense1 = Dense(32, activation='relu')(input_layer)
dense2 = Dense(16, activation='relu')(dense1)
dense3 = Dense(16, activation='relu')(dense2)
dropout_layer = Dropout(0.2)(dense3)
output_layer = Dense(1, activation='linear')(dropout_layer)
ts_model = Model(inputs=input_layer, outputs=output_layer)
ts_model.compile(loss='mean_absolute_error', optimizer='adam')
ts_model.summary()
Model: "functional_1"
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓ ┃ Layer (type) ┃ Output Shape ┃ Param # ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩ │ input_layer_1 (InputLayer) │ (None, 7) │ 0 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_4 (Dense) │ (None, 32) │ 256 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_5 (Dense) │ (None, 16) │ 528 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_6 (Dense) │ (None, 16) │ 272 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dropout_1 (Dropout) │ (None, 16) │ 0 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_7 (Dense) │ (None, 1) │ 17 │ └─────────────────────────────────┴────────────────────────┴───────────────┘
Total params: 1,073 (4.19 KB)
Trainable params: 1,073 (4.19 KB)
Non-trainable params: 0 (0.00 B)
save_weights_at = os.path.join('keras_models', 'PRSA_data_PM2.5_MLP_weights.{epoch:02d}-{val_loss:.4f}.keras')
save_best = ModelCheckpoint(save_weights_at, monitor='val_loss', verbose=0,
save_best_only=True, save_weights_only=False, mode='min', save_freq='epoch');
history_pm25 = None
if os.path.exists('history_pm25.joblib'):
history_pm25 = load('history_pm25.joblib')
print("El archivo 'history_pm25.joblib' ya existe. Se ha cargado el historial del entrenamiento.")
else:
history_pm25 = ts_model.fit(x=X_train, y=y_train, batch_size=16, epochs=20,
verbose=2, callbacks=[save_best], validation_data=(X_val, y_val),
shuffle=True);
dump(history_pm25.history, 'history_pm25.joblib')
print("El entrenamiento se ha completado y el historial ha sido guardado en 'history_pm25.joblib'.")
El archivo 'history_pm25.joblib' ya existe. Se ha cargado el historial del entrenamiento.
model_dir = 'keras_models'
files = os.listdir(model_dir)
pattern = r"PRSA_data_PM2.5_MLP_weights\.(\d+)-([\d\.]+)\.keras"
best_val_loss = float('inf')
best_model_file = None
best_model = None
for file in files:
match = re.match(pattern, file)
if match:
epoch = int(match.group(1))
val_loss = float(match.group(2))
if val_loss < best_val_loss:
best_val_loss = val_loss
best_model_file = file
if best_model_file:
best_model_path = os.path.join(model_dir, best_model_file)
print(f"Cargando el mejor modelo: {best_model_file} con val_loss: {best_val_loss}")
best_model = load_model(best_model_path)
else:
print("No se encontraron archivos de modelos que coincidan con el patrón.")
Cargando el mejor modelo: PRSA_data_PM2.5_MLP_weights.11-0.0117.keras con val_loss: 0.0117
preds = best_model.predict(X_val)
pred_pm25 = scaler.inverse_transform(preds)
pred_pm25 = np.squeeze(pred_pm25)
1/271 ━━━━━━━━━━━━━━━━━━━━ 34s 129ms/step
77/271 ━━━━━━━━━━━━━━━━━━━━ 0s 659us/step
169/271 ━━━━━━━━━━━━━━━━━━━━ 0s 597us/step
256/271 ━━━━━━━━━━━━━━━━━━━━ 0s 590us/step
271/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
271/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
mae = mean_absolute_error(df_val['pm2.5'].loc[7:], pred_pm25)
print('MAE for the validation set:', round(mae, 4))
MAE for the validation set: 11.6687
Trazamos los 50 primeros valores reales y predichos de
pm2.5
plt.figure(figsize=(5.5, 5.5))
plt.plot(range(50), df_val['pm2.5'].loc[7:56], linestyle='-', marker='*', color='r')
plt.plot(range(50), pred_pm25[:50], linestyle='-', marker='.', color='b')
plt.legend(['Actual','Predicted'], loc=2)
plt.title('Actual vs Predicted pm2.5')
plt.ylabel('pm2.5')
plt.xlabel('Index');
4.16. LSTM para la predicción de series de tiempo#
Continuaremos usando el conjunto de datos sobre contaminación del aire (pm2.5) usando LSTM. La lectura y el preprocesamiento de los datos se mantienen igual que en los ejemplos de
MLPs. El conjunto de datos original se divide en dos conjuntos: entrenamiento y validación, que se usan para el entrenamiento y validación del modelo respectivamente.
La función
makeXyse utiliza para generar arreglos de regresores y objetivos:X_train, X_val, y_trainyy_val.X_train y X_val, generados por la funciónmakeXy, son arreglos 2D de forma (número de muestras, número de pasos de tiempo). Sin embargo, la entrada a las capas RNN debe ser de forma (número de muestras, número de pasos de tiempo, número de características por paso de tiempo).En este caso, estamos tratando solo con
pm2.5, por lo tanto, el número de características por paso de tiempo es uno. El número de pasos de tiempo es siete y el número de muestras es el mismo que el número de muestras enX_trainyX_val, que se transforman a arreglos 3D
X_train, X_val = X_train.reshape((X_train.shape[0], X_train.shape[1], 1)), X_val.reshape((X_val.shape[0], X_val.shape[1], 1))
print('Shape of 3D arrays:', X_train.shape, X_val.shape)
Shape of 3D arrays: (33089, 7, 1) (8654, 7, 1)
X_trainyX_valse han convertido en matrices 3D y sus nuevas formas se ven en la salida de la instrucción de impresión anterior
from keras.layers import Dense, Input, Dropout
from keras.layers import LSTM
from keras.optimizers import SGD
from keras.models import Model
from keras.models import load_model
from keras.callbacks import ModelCheckpoint
La red neuronal para desarrollar el modelo de predicción de series temporales tiene una capa de entrada, que alimenta a la capa LSTM. La capa LSTM tiene siete pasos temporales, que es el mismo número de observaciones históricas tomadas para hacer la predicción de la presión atmosférica para el día siguiente. Solo el último paso temporal de la LSTM devuelve una salida.
Hay sesenta y cuatro neuronas ocultas en cada paso temporal de la capa LSTM. Por lo tanto, la salida de la LSTM tiene sesenta y cuatro características:
input_layer = Input(shape=(7,1), dtype='float32')
Aquí
shape=(7,1)forma de la entrada a la LSTM.7: Número de pasos de tiempo (timesteps).1: Número de características (features) por cada paso de tiempo.return_sequences=Truese usa cuando se necesita pasar la secuencia completa de salidas a la siguiente capa de la red, que también puede ser una capa recurrente como otra LSTM
lstm_layer1 = LSTM(64, input_shape=(7,1), return_sequences=True)(input_layer)
lstm_layer2 = LSTM(32, input_shape=(7,64), return_sequences=False)(lstm_layer1)
A continuación, la salida de LSTM pasa a una capa de exclusión que elimina aleatoriamente el 20% de la entrada antes de pasar a la capa de salida, que tiene una única neurona oculta con una función de activación lineal
dropout_layer = Dropout(0.2)(lstm_layer2)
output_layer = Dense(1, activation='linear')(dropout_layer)
Finalmente, todas las capas se envuelven en un
keras.models.Modely se entrenan durante veinte epochs para minimizar el MSE utilizando el optimizador Adam
ts_model = Model(inputs=input_layer, outputs=output_layer)
ts_model.compile(loss='mean_absolute_error', optimizer='adam') #SGD(lr=0.001, decay=1e-5))
ts_model.summary()
Model: "functional_2"
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓ ┃ Layer (type) ┃ Output Shape ┃ Param # ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩ │ input_layer_2 (InputLayer) │ (None, 7, 1) │ 0 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ lstm (LSTM) │ (None, 7, 64) │ 16,896 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ lstm_1 (LSTM) │ (None, 32) │ 12,416 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dropout_2 (Dropout) │ (None, 32) │ 0 │ ├─────────────────────────────────┼────────────────────────┼───────────────┤ │ dense_8 (Dense) │ (None, 1) │ 33 │ └─────────────────────────────────┴────────────────────────┴───────────────┘
Total params: 29,345 (114.63 KB)
Trainable params: 29,345 (114.63 KB)
Non-trainable params: 0 (0.00 B)
save_weights_at = os.path.join('keras_models', 'PRSA_data_PM2.5_LSTM_weights.{epoch:02d}-{val_loss:.4f}.keras')
save_best = ModelCheckpoint(save_weights_at, monitor='val_loss', verbose=0,
save_best_only=True, save_weights_only=False, mode='min', save_freq='epoch')
history_pm25_LSTM = None
if os.path.exists('history_pm25_LSTM.joblib'):
history_pm25_LSTM = load('history_pm25_LSTM.joblib')
print("El archivo 'history_pm25_LSTM.joblib' ya existe. Se ha cargado el historial del entrenamiento.")
else:
history_pm25_LSTM = ts_model.fit(x=X_train, y=y_train, batch_size=16, epochs=20,
verbose=2, callbacks=[save_best], validation_data=(X_val, y_val),
shuffle=True);
dump(history_pm25_LSTM.history, 'history_pm25_LSTM.joblib')
print("El entrenamiento se ha completado y el historial ha sido guardado en 'history_pm25_LSTM.joblib'.")
El archivo 'history_pm25_LSTM.joblib' ya existe. Se ha cargado el historial del entrenamiento.
model_dir = 'keras_models'
files = os.listdir(model_dir)
pattern = r"PRSA_data_PM2.5_LSTM_weights\.(\d+)-([\d\.]+)\.keras"
best_val_loss = float('inf')
best_model_file = None
best_model = None
for file in files:
match = re.match(pattern, file)
if match:
epoch = int(match.group(1))
val_loss = float(match.group(2))
if val_loss < best_val_loss:
best_val_loss = val_loss
best_model_file = file
if best_model_file:
best_model_path = os.path.join(model_dir, best_model_file)
print(f"Cargando el mejor modelo: {best_model_file} con val_loss: {best_val_loss}")
best_model = load_model(best_model_path)
else:
print("No se encontraron archivos de modelos que coincidan con el patrón.")
Cargando el mejor modelo: PRSA_data_PM2.5_LSTM_weights.10-0.0116.keras con val_loss: 0.0116
preds = best_model.predict(X_val)
pred_pm25 = scaler.inverse_transform(preds)
pred_pm25 = np.squeeze(pred_pm25)
1/271 ━━━━━━━━━━━━━━━━━━━━ 32s 121ms/step
39/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
73/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
109/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
144/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
181/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
219/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
259/271 ━━━━━━━━━━━━━━━━━━━━ 0s 1ms/step
271/271 ━━━━━━━━━━━━━━━━━━━━ 0s 2ms/step
271/271 ━━━━━━━━━━━━━━━━━━━━ 1s 2ms/step
from sklearn.metrics import mean_absolute_error
mae = mean_absolute_error(df_val['pm2.5'].loc[7:], pred_pm25)
print('MAE for the validation set:', round(mae, 4))
MAE for the validation set: 11.5364
Hemos utilizado
keras.callbacks.ModelCheckpointcomo callback para rastrear elMSEdel modelo en el conjunto de validación y guardar los pesos de la epoch que da el mínimo MSE. ElR-cuadradodel mejor modelo para el conjunto de validación es de 0.9959. Los primeros cincuenta valores reales y predichos se muestran en la siguiente figura
plt.figure(figsize=(5.5, 5.5))
plt.plot(range(50), df_val['pm2.5'].loc[7:56], linestyle='-', marker='*', color='r')
plt.plot(range(50), pred_pm25[:50], linestyle='-', marker='.', color='b')
plt.legend(['Actual','Predicted'], loc=2)
plt.title('Actual vs Predicted pm2.5')
plt.ylabel('pm2.5')
plt.xlabel('Index');
