Análisis Numérico de la resolución de ecuaciones diferenciales mediante redes neuronales artificiales

TRABAJO FIN DE M´ ASTER An´alisis Num´erico de la resoluci´on de ecuaciones diferenciales mediante redes neuronales artificiales. Realizado por Ana Casado S´anchez Para la obtenci´on del t´ıtulo de M´aster Universitario en Matem´aticas Dirigido por D˜na. Macarena G´omez M´armol Realizado en el departamento de Ecuaciones Diferenciales y An´alisis Num´erico Convocatoria de Julio, curso 2023/24 Agradecimientos En primer lugar, quiero agradecer a mi tutora Maca su atenci´on y su cari˜no a pesar de sus m´ultiples ocupaciones como profesora e investigadora a˜nadidas a sus labores como secretaria del departamento. Adem´as de dirigirme el trabajo, ha mostrado preocupaci´on por mi posible futuro como investigadora (y por m´ı a secas, con sus constantes Ana no te agobies, no pienses tanto en el futuro...), proporcion´andome toda clase de documentos para solicitudes de becas y d´andome la gran oportunidad de tener mi primera comunicaci´on a un congreso, exponiendo mis resultados en el CEDYA / CMA 2024. A continuaci´on, quiero agradecer a mis padres Patrocinio y Juan su amor y su constante compa˜n´ıa. En una ´epoca de incertidumbre como es mi incorporaci´on al mundo laboral tras finalizar la etapa de M´aster, sus consejos, sus bromas y, en general, su apoyo son vitales. Ambos han experimentado grandes p´erdidas este a˜no y a´un as´ı han dejado para m´ı su mejor parte. En relaci´on a estas p´erdidas, quiero darle las gracias a mi abuela To˜ni, cuyo humor y peculiar forma de pensamiento me acompa˜nan cada d´ıa. Asimismo, expreso mi infinita gratitud hacia el resto de miembros de mi familia, que nunca dejan de manifiestar su preocupaci´on por m´ı en la distancia, siendo plenamente conscientes de lo atareada que estoy y deseando siempre lo mejor para m´ı. En especial quiero dar las gracias a mi primo Juan Antonio con sus visitas semanales para cenar y ver pelis cuanto m´as antiguas mejor (seg´un mi padre). Sus bromas y su compa˜n´ıa lo convierten en lo segundo m´as cercano que tengo a un hermano peque˜no. Lo m´as cercano a esta figura fraternal es mi adorable perrito Epsilon, cuyas ocurrencias y lealtad lo convierten en un pilar fundamental de mi vida. Por ´ultimo, destaco a aquellas personas que sin v´ınculos de sangre son tambi´en parte de la misma, como Manolo y Faustino o Maria Jos´e, Javi, Sara y ´ Oscar en calidad de familia de coraz´on y amistades muy preciadas para m´ı como mi querida Ana. Otra parte de mi familia con la que no comparto v´ınculos de sangre (afortunadamente) es mi pareja y compa˜nero de aventuras Dani, quien siempre muestra su amor hacia m´ı y me cuida incluso cuando menos lo merezco, que es cuando m´as lo necesito. Me demuestra d´ıa a d´ıa que mi bienestar es una prioridad en su vida y espero hacerle sentir del mismo modo. Siempre sabe c´omo hacerme feliz y ayudarme a divertirme hasta en los momentos m´as dif´ıciles. Tengo mucha suerte de tener a alguien tan especial como ´el en mi vida y deseo seguir correspondi´endole con mi amor. i Resumen La presente memoria est´a dedicada al estudio de las redes neuronales artificiales aplicadas a la resoluci´on de problemas diferenciales debido a su creciente importancia en el campo de la simulaci´on num´erica. En particular, nos centramos en los modelos de redes neuronales informadas por la F´ısica, conocidos como PINN, para abordar estos problemas. Desde este enfoque, reducimos la resoluci´on de nuestros modelos a la de problemas de optimizaci´on. Asimismo, aportamos estimaciones del error cometido en funci´on de los hiperpar´ametros de la red (n´umero de neuronas, pesos y capas) al aproximar funciones en espacios de Sobolev utilizando funciones de activaci´on ReLU. Adem´as, desde el punto de vista de las aplicaciones, se explora la utilidad de estos m´etodos en diversos problemas relacionados con la teor´ıa de ecuaciones diferenciales. Esto incluye la resoluci´on de ecuaciones y sistemas diferenciales tanto lineales como no lineales, problemas de c´alculo de variaciones, problemas de control y problemas inversos. Para los experimentos num´ericos se emplean las herramientas software Matlab, FreeFem y Python, destacando su m´odulo TensorFlow. Palabras clave: Redes Neuronales Artificiales, Ecuaciones Diferenciales, An´alisis Num´erico, Optimizaci´on ii Abstract The present report is dedicated to the study of artificial neural networks applied to the resolution of differential problems due to their growing importance in the field of numerical simulation. In particular, we focus on Physics-Informed Neural Networks (PINNs) to address these problems. From this approach, we reduce the resolution of our models to optimization problems. Additionally, we provide error estimates based on the network’s hyperparameters (number of neurons, weights, and layers) when approximating functions in Sobolev spaces using ReLU activation functions. Moreover, from the perspective of applications, we explore the utility of these methods in various problems related to differential equation theory. This includes solving both linear and nonlinear differential equations and systems, calculus of variations problems, control problems, and inverse problems. For numerical experiments, the software tools Matlab, FreeFem, and Python are used, with a particular emphasis on its TensorFlow module. Keywords: Neural Networks, Differential Equations, Numerical Analysis, Optimization iii ´ Indice general 1. Introducci´on 1 2. Introducci´on a las redes neuronales artificiales (ANN) 6 2.1. Construcci´on de las ANN . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2. Entrenamiento de las ANN . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1. Dise˜no de la red . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2. Aprendizaje de las ANN. Tipos . . . . . . . . . . . . . . . . . 11 2.3. Funciones de activaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4. Algoritmo de descenso por el gradiente . . . . . . . . . . . . . . . . . 18 2.5. Ventajas de las ANN sobre otras t´ecnicas de Machine Learning . . . . 21 3. Resolviendo problemas de ecuaciones diferenciales con ANN 23 3.1. Un m´etodo de resoluci´on de problemas diferenciales con ANN . . . . 23 3.1.1. Construcci´on de las soluciones de prueba Ψp.......... 25 3.1.2. Optimizaci´on del error . . . . . . . . . . . . . . . . . . . . . . 28 3.2. Primeros experimentos . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4. Estimaciones de error 40 4.1. Cotas superiores sobre el n´umero de elementos de una ANN para aproximarfunciones............................ 40 4.2. Cotas inferiores sobre el n´umero de elementos de una ANN para aproximarfunciones.............................. 49 4.3. Aproximando funciones con redes neuronales artificiales superficiales . 55 4.4. Test de validaci´on de las cotas de errores . . . . . . . . . . . . . . . . 58 5. Algunos modelos implementados con Python 63 5.1. Esquema de resoluci´on de problemas de optimizaci´on con ANN. Optimizadores ................................ 63 5.2. El problema de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.2.1. Ejemplo 1: problema con segundo miembro constante en el cubounidad............................ 66 5.2.2. Ejemplo 2: problema con segundo miembro constante en la bolaunidad ............................ 74 5.2.3. Ejemplo 3: problema con segundo miembro no constante en el cubounidad............................ 78 5.2.4. Ejemplo 4: problema con condiciones de Neumann . . . . . . . 81 5.2.5. Ejemplo 5: problema no lineal . . . . . . . . . . . . . . . . . . 85 5.3. El problema de Laplace como problema de c´alculo de variaciones . . . 91 5.4. Problemadecontrol ........................... 95 5.5. Problema evolutivo: ecuaci´on log´ıstica con difusi´on . . . . . . . . . . 104 5.6. Problema inverso: modelo de Lotka-Volterra . . . . . . . . . . . . . . 112 iv 6. Conclusiones y trabajos futuros 127 7. Bibliograf´ıa 128 v 1. Introducci´on El nacimiento de la teor´ıa de redes neuronales artificiales (Artifitial Neural Networks, ANN) se sit´ua en 1943, con la formulaci´on del primer modelo matem´atico para simular el comportamiento de las neuronas del cerebro humano de la mano de McCulloch y Pitts. Con el desarrollo de este modelo, los autores buscaban comprender el mecanismo mediante el cual nuestras neuronas y sus interconexiones son capaces de procesar la informaci´on recibida mediante est´ımulos y transformarla en pensamientos complejos. Posteriormente, en 1950, con la formulaci´on por Alan Turing del Test que lleva su nombre, se introdujo el concepto de una m´aquina que pudiera aprender y pensar como un hombre. Esto supuso el nacimiento de la Inteligencia Artificial (IA) e indujo al desarrollo de diversas t´ecnicas de estudio en este ´ambito durante esta d´ecada, entre las que destacamos el Aprendizaje Autom´atico (Machine Learning, ML) en el que se desarrollan algoritmos para que las m´aquinas puedan resolver problemas por s´ı mismas a partir de datos. As´ı, las redes neuronales artificiales comenzaron a desarrollarse como una t´ecnica perteneciente a este campo, destacando el nacimiento del perceptr´on o neurona artificial, tal y como la conocemos hoy en d´ıa gracias a Rosenblatt, lo que permiti´o el desarrollo de los conocidos como modelos monocapa conectando perceptrones en paralelo. Sin embargo, la publicaci´on de Perceptrons por parte de Minsky y Papert en 1960, un libro que demostraba las limitaciones de las ANN monocapa junto con la ausencia de algoritmos para trabajar con redes neuronales con conexiones en serie o multicapa, provoc´o que la comunidad cient´ıfica trasladase su inter´es a otras t´ecnicas de Aprendizaje Autom´atico, como los m´etodos de agrupamiento y el reconocimiento de patrones estad´ısticos, dando lugar al per´ıodo conocido como el primer invierno de la IA. No fue hasta la d´ecada de los 80 con la introducci´on de los m´etodos de retropropagaci´on, i.e. la aplicaci´on de m´etodos de descenso por gradiente a ANN, y su capacidad para trabajar con ANN multicapa, que se produjo un resurgimiento en la teor´ıa de redes neuronales artificiales. Por otra parte, durante la segunda mitad del siglo XX, la aparici´on de los computadores y su utilidad como herramientas de c´alculo, tambi´en produjo un gran desarrollo en otra l´ınea de estudio en matem´aticas, el An´alisis Num´erico, que nos permite resolver de forma aproximada modelos complejos, por ejemplo, basados en ecuaciones diferenciales (en la mayor´ıa de los casos no es conocida una soluci´on anal´ıtica). Las ecuaciones diferenciales modelan una gran cantidad de problemas con aplicaciones en m´ultiples campos, como F´ısica, Ingenier´ıa, Biolog´ıa, Ciencias Sociales... De este modo, en la segunda mitad de la d´ecada de los 90, ambas disciplinas 1 se entrelazaron para el estudio de la resoluci´on de problemas diferenciales utilizando ANN. Este es precisamente el objetivo del trabajo realizado: una introducci´on al An´alisis Num´erico de la resoluci´on de ecuaciones diferenciales mediante redes neuronales artificiales. En particular, trataremos diversos problemas utilizando los conocidos como modelos de redes neuronales informadas por la F´ısica, conocidos como PINNs. Desde este enfoque, reduciremos nuestros modelos a la resoluci´on de problemas de optimizaci´on en espacios de Sobolev para minimizar el error cuadr´atico medio cometido respecto a las restricciones de nuestro problema (ecuaci´on y condiciones de contorno y/o iniciales) para ajustar los par´ametros variables de la red neuronal (pesos y b´ıas). En consecuencia, problemas que requieren bucles iterativos adicionales para el M´etodo de los Elementos Finitos, como por ejemplo es el caso de los problemas de ecuaciones diferenciales no lineales o de los evolutivos, se tratan de forma similar al resto de problemas mediante ANN, sin bucles adicionales. Una caracter´ıstica de las PINN es que, a diferencia de los m´etodos num´ericos cl´asicos, no necesitamos una discretizaci´on del espacio espec´ıfica para que los m´etodos converjan, como por ejemplo, ocurre con los m´etodos de Diferencias Finitas y Elementos Finitos. De hecho, es muy frecuente utilizar nubes de puntos aleatorias para realizar experimentos num´ericos. Desde el punto de vista computacional, la resoluci´on de estos problemas conlleva un alto coste debido a que los m´etodos de descenso por gradiente tienen una convergencia muy lenta y costosa. Sin embargo, el ajuste de par´ametros realizado para la resoluci´on de los distintos problemas, nos proporciona una funci´on susceptible a ser evaluada en cualquier punto del dominio, la cual podemos almacenar para su posterior uso sin necesidad de c´alculos adicionales. Incluso podemos partir de modelos previamente entrenados para resolver problemas similares ahorrando c´alculos. En relaci´on al alto coste computacional, hay que mencionar que debido a la ausencia de un hardware lo suficientemente potente para procesar tal cantidad de c´alculos, se produjo un segundo invierno de la IA, lo que provoc´o que el estudio de las ANN y, en particular, de las PINNs se abandonase para retomarse nuevamente a finales de los 2000 con la aparici´on de equipos capaces de procesar grandes vol´umenes de datos y al uso de la GPU para realizar c´alculos. Adem´as, hay librer´ıas, como por ejemplo TensorFlow, disponibles para lenguajes de programaci´on de alto nivel, entre los que destacamos Python, Matlab, C++ y Java. Desde entonces, las redes neuronales han ido tomando una importancia cada vez mayor, gracias a su capacidad para aprender y generalizar a partir de datos. En la actualidad, destacan sus aplicaciones en el reconocimiento de voz y de im´agenes, procesamiento del lenguaje natural, predicci´on y an´alisis financiero, diagn´osticos m´edicos, conducci´on aut´onoma, personalizaci´on de contenido en plataformas digitales, y optimizaci´on en la gesti´on de recursos y energ´ıas. Su versatilidad las convierte en herramientas esenciales para la innovaci´on y la eficiencia en m´ultiples industrias. 2 En el ´ambito de las Ecuaciones Diferenciales y el An´alisis Num´erico tambi´en est´an captando el interes de una gran cantidad de investigadores. Asimismo, un importante objeto de estudio es la obtenci´on de estimaciones acerca de los errores cometidos en las aproximaciones en funci´on de los hiperpar´ametros de la red (n´umero de capas, pesos y neuronas), aunque actualmente es un campo en el que queda una gran labor investigadora por delante. Adem´as, las estimaciones que conocemos ´unicamente son aplicables suponiendo que el algoritmo de retropropagaci´on haya encontrado el m´ınimo global del problema, lo cual en general no ser´a cierto al no estar en condiciones de convexidad adecuadas para poder asegurar la convergencia a un m´ınimo global del problema. En cuanto a las PINNs, comentamos que actualmente, tambi´en se est´an explorando algunas variantes, como son las VPINNs, en las que introducimos el error cometido respecto a la formulaci´on variacional de los problemas a tratar como parte del funcional a minimizar por el m´etodo de descenso. Por su parte, la memoria est´a estructurada como sigue: En el Cap´ıtulo 2, introducimos la teor´ıa correspondiente al campo de las ANN necesaria para poder desarrollar la memoria. En particular, exponemos los elementos fundamentales que conforman una neurona artificial o perceptr´on (entradas, pesos, polarizaciones o b´ıas, funci´on de activaci´on y salida), realizando una analog´ıa con la estructura de las neuronas biol´ogicas. Asimismo, mostramos c´omo construimos redes conectando perceptrones en serie y en paralelo, as´ı como dedicamos una secci´on a las funciones de activaci´on, donde justificamos su importancia para poder aproximar funciones que no sean lineales y presentamos las m´as utilizadas en orden cronol´ogico, resaltando las funciones sigmoide, tangente hiperb´olica y la parte positiva o ReLU (Rectified Linear Unit). Una vez descritos estos elementos, procedemos a definir el entrenamiento de una ANN como un ajuste de sus pesos y b´ıas para producir salidas adecuadas a los datos de entrada para una arquitectura prefijada. As´ı, comentamos tambi´en algunas heur´ısticas utilizadas a la hora de dise˜nar una ANN para su posterior entrenamiento. En adici´on, exponemos los conceptos de aprendizaje supervisado y no supervisado, seg´un conozcamos las salidas esperadas para los datos de entrada o no, respectivamente, as´ı como exponemos los tipos de aprendizaje m´as utilizados en cada una de estas categor´ıas. En este sentido, destacamos las PINNs como parte del aprendizaje no supervisado sujeto a restricciones, ya que son problemas en los que no conocemos la salida esperada pero s´ı propiedades que debe verificar la misma. Es por esto que este tipo de aprendizaje tambi´en suele denominarse aprendizaje supervisado d´ebil y hay fuentes que prefieren considerarlo como un tipo de aprendizaje por s´ı mismo situado a medio camino entre el supervisado y no supervisado. Otro aspecto a comentar en relaci´on al entrenamiento es que, tal y como se ha explicado, se basa en m´etodos de descenso por gradiente, por lo que hemos dedicado una secci´on a recordarlos y enunciar resultados sobre su convergencia. Concluimos este cap´ıtulo incluyendo algunas de las ventajas que poseen las 3 Figura 2.6: Red neuronal multicapa. 2.2. Entrenamiento de las ANN Al inicio de este cap´ıtulo, se coment´o que las redes neuronales pretenden crear una forma de inteligencia artificial capaz de tomar decisiones por s´ı misma. Asimismo, se ha visto que los datos de entrada se procesan a trav´es de la red neuronal con el prop´osito de lograr una salida. De este modo, el objetivo de la red neuronal es aprender a calcular la salida correcta para cada ejemplo de entrada. En esta secci´on, estudiamos c´omo hacerlo. En primer lugar, establecemos que la topolog´ıa de la red y las diferentes funciones de cada neurona (entrada, activaci´on y salida), conocidos como hiperpar´ametros de la red, no pueden cambiar durante el aprendizaje, mientras que los pesos sobre cada una de las conexiones s´ı pueden hacerlo. As´ı pues, es necesario establecer una topolog´ıa para la red adecuada antes de comenzar a entrenarla. 2.2.1. Dise˜no de la red Actualmente, para la mayor´ıa de problemas no se sabe determinar qu´e arquitectura tomar a priori antes de enfrentarse a un problema. De hecho, en el ´ambito que nos interesa, el de las ecuaciones diferenciales, este es el principal problema. Apenas hay resultados que especifiquen la influencia que tiene en el error tomar un n´umero de capas o de neuronas determinado. Este es el objetivo de esta l´ınea de estudio en el an´alisis num´erico. Sin embargo, hay algunas heur´ısticas de selecci´on relacionadas con la dimensionalidad de las capas: 10 Es m´as rentable a˜nadir capas que incrementar la dimensionalidad por capa, por lo que priorizaremos lo primero. Si partimos de una representaci´on del problema con entradas con dimensionalidad baja, iremos a˜nadiendo capas aumentando la dimensionalidad de cada una frente a la anterior paulatinamente hasta alcanzar una dimensionalidad suficientemente alta. Finalmente, en las capas sucesivas, aplicaremos el procedimiento inverso. As´ı, primero aumentamos la cantidad de informaci´on aprendida y, a continuaci´on, para evitar un sobreajuste de los datos, forzaremos a aprender lo realmente clave para nuestro problema. Figura 2.7: Arquitectura de ANN con dimensionalidad creciente en las capas iniciales y decreciente en las posteriores. Si partimos de una representaci´on del problema con entradas con dimensionalidad alta, partimos de un exceso de informaci´on, luego para evitar sobreajustes, iremos reduciendo paulatinamente la dimensionalidad de una capa a la siguiente. Figura 2.8: Arquitectura de ANN con dimensionalidad decreciente. 2.2.2. Aprendizaje de las ANN. Tipos Una vez escogida la topolog´ıa de la red, para conseguir las mejores salidas para cada entrada, los ´unicos par´ametros que podremos modificar son los pesos y 11 polarizaciones. Por tanto, el aprendizaje de una ANN se traduce en una adaptaci´on de sus pesos y b´ıas. Esto nos lleva a la siguiente definici´on: Definici´on 2.7. El aprendizaje oentrenamiento de una red neuronal artificial es el proceso por el cual la red neuronal modifica sus pesos en respuesta a una informaci´on de entrada. Se considerar´a que la red est´a entrenada cuando los valores de los pesos permanezcan estables. Observaci´on 2.8. Los cambios que se producen durante el entrenamiento de una ANN se reducen a la destrucci´on, modificaci´on y creaci´on de conexiones entre las neuronas. La creaci´on de una nueva conexi´on implica que el peso de la misma pasa a tener un valor distinto de cero. De la misma manera, una conexi´on se destruye cuando su peso pasa a ser cero. Para cualquier algoritmo de entrenamiento, siempre comenzaremos dando una asignaci´on de pesos inicial. Esto puede realizarse por varios criterios; por ejemplo uno de ellos es asignar un peso aleatorio perteneciente a un intervalo dado a cada conexi´on. A continuaci´on, entrenaremos la red hasta que se verifique una condici´on de parada establecida, como puede ser el establecimiento de un valor de tolerancia para el error cuadr´atico medio o de un n´umero de iteraciones prefijado. De este modo, el siguiente paso para poder trabajar con ANN es comprender c´omo entrenarlas. Fundamentalmente, se distinguen dos l´ıneas de aprendizaje: Definici´on 2.9. El aprendizaje supervisado se caracteriza porque el proceso de aprendizaje se realiza mediante un entrenamiento controlado por un agente externo (supervisor, maestro) que determina la respuesta que deber´ıa generar la red a partir de una entrada determinada. El supervisor controla la salida de la red y en caso de que ´esta no coincida con la deseada, se proceder´a a modificar los pesos de las conexiones, con el fin de conseguir que la salida obtenida se aproxime a la deseada. En este tipo de aprendizaje se suelen considerar, a su vez, tres formas de llevarlo a cabo, que dan lugar a los siguientes aprendizajes supervisados: 1. Aprendizaje por correcci´on de error: consiste en ajustar los pesos de las conexiones de la red en funci´on de la diferencia entre los valores deseados y los obtenidos a la salida de la red, es decir, en funci´on del error cometido en la salida. 2. Aprendizaje por refuerzo: al no disponer de un ejemplo completo del comportamiento deseado, el supervisor no indica durante el entrenamiento exactamente la salida que se desea que proporcione la red ante una determinada entrada. As´ı, la funci´on del supervisor se reduce a indicar mediante una se˜nal de refuerzo si la salida obtenida en la red se ajusta a la deseada (´exito = +1 o fracaso = -1) y en funci´on de ello se ajustan los pesos 3. Aprendizaje estoc´astico: consiste b´asicamente en realizar cambios aleatorios en los valores de los pesos de las conexiones de la red y evaluar su efecto a partir del objetivo deseado y de distribuciones de probabilidad. 12 Definici´on 2.10. Las redes con aprendizaje no supervisado (tambi´en conocido como autosupervisado) no requieren influencia externa para ajustar los pesos de las conexiones entre sus neuronas. La red no recibe ninguna informaci´on por parte del entorno que le indique si la salida generada en respuesta a una determinada entrada es o no correcta. Estas redes deben encontrar por s´ı mismas caracter´ısticas sobre los datos de entrada que les permitan ajustar los par´ametros de la red adecuadamente para producir sus respectivas salidas. Existen varias posibilidades en cuanto a la interpretaci´on de la salida de estas redes, que dependen de su estructura y del algoritmo de aprendizaje empleado. En cuanto a los algoritmos de aprendizaje no supervisado, en general se suelen considerar: 1. Aprendizaje hebbiano: esta regla de aprendizaje es la base de muchas otras; pretende medir la familiaridad o extraer caracter´ısticas de los datos de entrada. Se fundamenta en la siguiente suposici´on: si dos neuronas NiyNjtoman el mismo estado simult´aneamente (ambas activas o ambas inactivas), el peso de la conexi´on entre ambas se incrementa. En este caso, las entradas y salidas permitidas a la neurona son: {−1,1}o{0,1}(neuronas binarias). 2. Aprendizaje competitivo y comparativo: se orienta a la clusterizaci´on o clasificaci´on de los datos de entrada. La caracter´ıstica principal del aprendizaje competitivo es que, si un patr´on nuevo se determina que pertenece a una clase reconocida previamente, entonces la inclusi´on de este nuevo patr´on a esta clase matizar´a la representaci´on de la misma. Si el patr´on de entrada se determin´o que no pertenece a ninguna de las clases reconocidas anteriormente, entonces la estructura y los pesos de la red neuronal ser´an ajustados para reconocer la nueva clase. 3. Aprendizaje basado en restricciones: el aprendizaje basado en restricciones o aprendizaje supervisado d´ebil es una t´ecnica en la que, a pesar de que no se conocen las salidas de los datos de entrada para el entrenamiento, el modelo se entrena no solo utilizando datos de entrada, sino tambi´en incorporando un conjunto de restricciones o reglas que deben ser satisfechas. Estas restricciones pueden provenir de leyes f´ısicas, reglas l´ogicas, propiedades matem´aticas o cualquier otro conocimiento a priori que describa el comportamiento del sistema bajo estudio. Un ejemplo son las Redes Neuronales Informadas por la F´ısica (Physics- Informed Neural Networks, PINNs) que incorporan leyes f´ısicas (ecuaciones diferenciales) que han de satisfacerse. Dado que el objeto de esta memoria es la resoluci´on de problemas diferenciales, trabajaremos con PINNs. 2.3. Funciones de activaci´on Para la construcci´on de las ANN expuesta en la secci´on 2.1, hemos introducido funciones de activaci´on para cada neurona, una parte fundamental de la estructura 13 del perceptr´on y, por consiguiente, de la topolog´ıa de una red neuronal artificial. Esto se debe a que, en ausencia de estas, por mucho que a˜nadi´esemos nodos y/o capas, las ´unicas operaciones que se realizar´ıan ser´ıan multiplicaciones de matrices, es decir, composiciones de aplicaciones lineales, lo que es a su vez una aplicaci´on lineal. En consecuencia, las ANN solo podr´ıan aproximar funciones lineales y, por tanto, su rango de aplicaciones ser´ıa muy limitado. De este modo, el objetivo de las funciones de activaci´on es romper dicha linealidad. De hecho, en la presente memoria queremos aproximar las soluciones de problemas diferenciales, i.e., buscamos aproximar funciones de diversa ´ındole. Esto ser´ıa imposible sin funciones de activaci´on. A continuaci´on, dedicamos esta secci´on a exponer las funciones de activaci´on m´as importantes cronol´ogicamente. En la d´ecada de los 50, que fue cuando comenz´o a desarrollarse la teor´ıa de redes neuronales artificiales de la mano de McColluch-Pitts, se empleaban Funci´on escal´on unitario o de Heavyside: f(x) = 0,si x < 0, 1,si x≥0 Figura 2.9: Funci´on de Heavyside. Funci´on signo: f(x) = −1,si x < 0, 1,si x≥0 14 Figura 2.10: Funci´on signo. La raz´on de escoger dichas funciones poder separar los datos en dos conjuntos de informaci´on separados. Sin embargo, al ser discontinuas, su uso era muy limitado y adem´as no pose´ıan esquemas de entrenamiento eficientes para entrenar las redes multicapa, por lo que a inicios de los 70 las ANN dejaron de utilizarse. No obstante, a finales de los 80, con la incorporaci´on del algoritmo de backpropagation o descenso por el gradiente (a los que dedicaremos una secci´on), resurgi´o este campo incorporando funciones que no son m´as que regularizaciones de las anteriores para que fueran infinitamente derivables y fueran aptas para dichos algoritmos: Sigmoide: f(x) = σ(x) = 1 1 + e−x Figura 2.11: Funci´on sigmoide. Tangente hiperb´olica: f(x) = ex−e−x ex+e−x 15 Figura 2.12: Funci´on tangente hiperb´olica. Esto permiti´o resolver una gran cantidad de problemas. De hecho, en el cap´ıtulo siguiente, cuando estudiemos la primera formulaci´on de un m´etodo de resoluci´on de problemas de ecuaciones diferenciales utilizando ANN, que es del 98, se utiliza la funci´on sigmoide. A pesar de los avances conseguidos, los c´alculos eran muy lentos computacionalmente, lo que junto con el d´eficit de capacidad de componentes electr´onicos de la ´epoca y los problemas propios del algoritmo de backpropagation (que comentaremos en su secci´on correspondiente), dieron lugar a otro per´ıodo de pausa en el estudio de las redes neuronales. Hace unos 10-15 a˜nos, con los avances tecnol´ogicos correspondientes, es cuando ha tomado m´as fuerza el uso de redes neuronales y por lo que cada vez se investiga m´as en este ´ambito. En este resurgimiento, la funci´on de activaci´on por antonomasia gracias a que es muy r´apida de entrenar, es la parte positiva, a la que los investigadores suelen referirse como ReLU (Rectified Linear Unit) y variantes de esta. ReLU (Rectified Linear Unit): f(x) = 0,si x < 0, x, si x≥0= max(0, x) Figura 2.13: Funci´on ReLU. 16 El problema principal que presenta la ReLU es el conocido como dying ReLU, que consiste en que durante el entrenamiento de la red puede llegar un momento a partir del cual los valores sean negativos y no salgamos de la parte izquierda de la curva. Para paliar esto, dado α > 0, las variantes de la ReLU m´as utilizadas son: Leaky ReLU: f(x) = αx, si x < 0, x, si x≥0= max(αx, x) Figura 2.14: Funci´on Leaky ReLU con α=Leak. ELU (Exponential ReLU): se comenz´o a extender entre la comunidad en 2015, aporta mejores resultados que la ReLU pero es m´as lenta de computar. f(x) = α(ex−1),si x < 0, x, si x≥0 Figura 2.15: Funci´on ELU con α= 1. SELU (Scaled ELU): introducimos otro par´ametro λ > 0 f(x) = λα(ex−1),si x < 0, λx, si x≥0 17 Figura 2.16: Funci´on SELU con α= 1,67, λ= 1,05. Por ´ultimo, hay que comentar que mientras que las funciones signo y de Heavyside actualmente solo se utilizan de forma te´orica, el sigmoide y la tangente hiperb´olica s´ı suelen utilizarse como funciones de activaci´on para la capa de salida. De este modo, se combinan con la ReLU y sus variantes que act´uan como funciones de activaci´on para los nodos de la capa oculta. Con esto, concluimos el estudio de la estructura de las ANN y los elementos que las forman. El siguiente paso es aprender a utilizarlas para resolver problemas. 2.4. Algoritmo de descenso por el gradiente Dedicamos una secci´on al descenso por gradiente o retropropagaci´on (backpropagation), ya que va a ser fundamental en los m´etodos de aproximaci´on para la resoluci´on de ecuaciones diferenciales utilizados en esta memoria. Recordamos que se trata de un algoritmo para resolver un problema de optimizaci´on en un abierto U⊂Rd. En el caso particular de las ANN, se utiliza para actualizar los pesos buscando minimizar el error cometido con respecto a mediciones, en el caso del aprendizaje supervisado, o en las restricciones, en el caso del aprendizaje supervisado d´ebil. Antes de definir este algoritmo, recordamos la idea que induce al mismo. Consideremos un funcional derivable J:U→Ryx0∈U. El problema es encontrar un punto xdonde se alcance el m´ınimo de J. Realmente, en general, lo que vamos a obtener es una sucesi´on {xn} ⊂ U(que suponi´endola acotada) va a converger, para una subsucesi´on, a un punto xtal que ∇J(x) = 0. Puesto que llegaremos a x decreciendo, no podr´a ser un punto de m´aximo, pero s´ı un punto de silla. La construcci´on de la sucesi´on se realizar´a como sigue: dado un punto xn, queremos generar el siguiente elemento de la sucesi´on xn+1 de forma que se verifique J(xn)> J(xn+1). Buscamos xn+1 de la forma xn+1 =xn+τnvn, donde vn∈Rd unitario y τn∈R“peque˜no” para que xn+1 ∈U. Usando el desarrollo de Taylor de 18 primer orden se tiene J(xn+1)∼J(xn) + τn∇J(xn)·vn. La idea es tomar τnyvnminimizando el segundo miembro. Gracias a la desigualdad de Cauchy-Schwarz, la mejor elecci´on de vnviene dada por vn=−∇J(xn) |∇J(xn)|, donde la norma es la norma eucl´ıdea en Rd. El problema es la elecci´on de τn. Seg´un la estrategia seguida para la elecci´on de este par´ametro, existen diversos algoritmos de tipo gradiente. Observemos tambi´en que el c´alculo de ∇J(xn) puede no ser f´acil en algunos casos y entonces utilizaremos aproximaciones. En general, estos algoritmos se suelen escribir de la siguiente forma para una tolerancia prefijada a la que llamaremos tol > 0: Inicializaci´on. Tomar x0∈U. Mientras |∇J(xn)|> tol, se toma xn+1 =xn−ηn∇J(xn) con ηn∈Rsuficientemente peque˜no para que xn+1 ∈U. Como podemos ver, se trata de la idea explicada anteriormente con ηn=τn |∇J(xn)|. En el contexto de las ANN, cada ηnes conocido como la tasa de aprendizaje en la iteraci´on n-´esima. Las dos elecciones m´as cl´asicas de ηnson las siguientes: Algoritmo de descenso por gradiente con paso ´optimo: Consiste en tomar ηnsoluci´on de m´ın xn−η∇J(xn)∈UJ(xn−η∇J(xn)). Esto conlleva la resoluci´on de un problema de m´ınimo en una recta que puede no ser f´acil en la pr´actica salvo en casos muy simples como el correspondiente a J cuadr´atico. Algoritmo de descenso por gradiente con paso fijo: Consiste en tomar ηn=η con η > 0 fijo. Se pueden probar los siguientes resultados (consultar [1]): Teorema 2.11. Supongamos U=Rdy que existe α > 0tal que ∇J(v)− ∇J(u)·v−u≥α|v−u|2,∀u, v ∈Rd.(2.1) Entonces Jtiene un ´unico punto de m´ınimo ¯xy para todo x0∈Rdla sucesi´on {xn} generada por el m´etodo del gradiente con paso ´optimo converge a ¯x. 19 Ejemplo con Ecuaciones en Derivadas Parciales: Por simplificar la exposici´on, nos limitamos al caso de la ecuaci´on de Poisson en un cuadrado: −∆Ψ(x, y) = f(x, y),(x, y)∈[0,1]2, f ∈C0[0,1]2. En la memoria, tratamos m´as tarde con ecuaciones m´as generales, si bien el tratamiento que haremos sobre las restricciones del problema ser´a distinto al aqu´ı expuesto. Recordamos que hasta ahora estamos exponiendo el m´etodo propuesto en [8]. Buscamos las soluciones de prueba de la forma: Ψp(x, y) = A(x, y) + g(x, y)N(x, y, p), donde A, g ∈C2([0,1]2) depender´an de las condiciones de contorno que se hayan impuesto. Resolvemos el problema de optimizaci´on min pX xi,yi∈ˆ D∆Ψp(xi, yi) + f(xi, yi)2. Condiciones de contorno de tipo Dirichlet: Para f0, g0, f1, g1∈C2([0,1]) verificando        f0(0) = g0(0), f1(0) = g0(1), f0(1) = g1(0), f1(1) = g1(1), −(g′′ 0(0) + f′′ 0(0)) = f(0,0),−g′′ 1(0) −f′′ 0(1) = f(0,1), −g′′ 1(1) −f′′ 1(1) = f(1,1),−g′′ 0(1) −f′′ 1(0) = f(1,0), (3.4) consideramos las condiciones de contorno Ψ(0, y) = f0(y),Ψ(1, y) = f1(y), y ∈[0,1] Ψ(x, 0) = g0(x),Ψ(x, 1) = g1(x), x ∈[0,1]. Escribimos las soluciones de prueba como Ψp(x, y) = A(x, y) + g(x, y)N(x, y, p), donde gse anula en la frontera y en ning´un otro punto. A t´ıtulo de ejemplo, tomamos g(x, y) = xy(1 −x)(1 −y). Adem´as, A(x, y) se define para que satisfaga las condiciones de contorno, por ejemplo: A(x, y) = (1 −x)f0(y) + xf1(y) +(1 −y)[g0(x)−(1 −x)g0(0) −xg0(1)] +y[g1(x)−(1 −x)g1(0) −xg1(1)] 26 Condiciones de contorno mixtas: vamos a ver c´omo escribir Ψppara problemas con condiciones de tipo Dirichlet en una parte de la frontera y de Neumann en el resto. Por ejemplo, para f0, g0, f1∈C2([0,1]), g1∈C1([0,1]) verificando f0(0) = g0(0), g0(1) = f1(0), f′ 0(1) = g1(0), f′ 1(1) = g1(1), −(g′′ 0(0) + f′′ 0(0)) = f(0,0),−g′′ 0(1) −f′′ 1(0) = f(1,0), supongamos restricciones del tipo (Ψ(0, y) = f0(y),Ψ(1, y) = f1(y), y ∈[0,1] Ψ(x, 0) = g0(x), ∂yΨ(x, 1) = g1(x), x ∈[0,1] entonces escribimos las soluciones de prueba como Ψp(x, y) = A(x, y) + x(1 −x)yN(x, y, p)−N(x, 1, p)−∂N(x, 1, p) ∂y , donde A(x, y) se define para que satisfaga las condiciones de contorno, por ejemplo: A(x, y) = (1 −x)f0(y) + xf1(y) + g0(x) −(1 −x)g0(0) −xg0(1) + y[g1(x)−(1 −x)g1(0) −xg1(1)]. En el resto de casos, la construcci´on es an´aloga. Observaci´on 3.3. Las condiciones que hemos impuesto a los datos en los problemas anteriores para la ecuaci´on de Poisson son necesarias para la existencia de una soluci´on cl´asica con regularidad C2([0,1]2). Sin embargo recordamos que estas condiciones no son suficientes. En general, para tener una soluci´on de clase C2para la ecuaci´on de Poisson no basta con un segundo miembro continuo. Usando las estimaciones de Schauder, lo usual es tomar un segundo miembro en C0,α((0,1)2) con 0 < α ≤1 lo que asegura que la soluci´on se encuentre en C2,α((0,1)2). Sin embargo incluso pidiendo (adem´as de las condiciones de compatibilidad anteriores) f∈C2,α([0,1]2) y fi, gi∈C2,α([0,1]), no se tiene, en general, que la soluci´on se encuentre en C2,α([0,1]2), debido a la no regularidad del cuadrado por la presencia de las esquinas. Referimos a [10] para estos resultados. Observaci´on 3.4. El problema de minimizaci´on considerado en el algoritmo de redes neuronales puede ser considerado como una discretizaci´on de m´ın  ∆Ψ + f∥2 L2((0,1)2 cuando Ψ verifica las condiciones de contorno correspondientes. El espacio natural donde buscar una soluci´on asociada a esta formulaci´on ser´ıa H2((0,1)2) que aparece por tanto como el espacio natural asociado a la resoluci´on de problemas el´ıpticos usando redes neuronales. Siguiendo [10], esto es lo que se conoce como soluci´on fuerte del problema, la ecuaci´on se verifica en casi todo. 27 En el caso del problema de Dirichlet, para que la soluci´on est´e en H2((0,1)2), basta f∈L2((0,1)2) y que las condiciones de frontera definan una funci´on en H3 2(∂(0,1)2) (lo cual en particular implica que se cumplan las condiciones de compatibilidad en (3.4) que solo hacen intervenir las funciones y no las derivadas), los resultados de regularidad para ecuaciones el´ıpticas garantizan efectivamente la existencia de soluci´on en H2((0,1)2) gracias a la convexidad del cuadrado. En el caso del problema con condiciones mixtas, de nuevo la existencia de soluci´on en H2((0,1)2) no es clara debida a las esquinas (ver [11]) 3.1.2. Optimizaci´on del error La minimizaci´on del error se va a conseguir entrenando la red neuronal hasta obtener los pesos y b´ıas que permitan que para cada entrada xi, la cantidad G(xi,Ψp(xi),∇Ψp(xi),∇2Ψp(xi),...,∇kΨp(xi)) sea lo m´as cercana a cero posible. Para ello, en esta memoria aplicaremos el m´etodo de descenso por el gradiente o retropropagaci´on, aunque son v´alidas otras t´ecnicas como por ejemplo, el BFGS quasi-Newtoniano que puede encontrarse en [9]. Para aplicar dichos m´etodos, en particular necesitaremos las derivadas de la red con respecto a sus par´ametros (pesos y bias). La topolog´ıa escogida para la ANN N(x, p) en [8] es la de una red multicapa con una ´unica capa oculta formada por Hnodos cuyas funciones de activaci´on son el sigmoide. Como se coment´o en la secci´on anterior, es una ANN con un ´unico nodo de salida, al que le asignamos como funci´on de activaci´on la identidad. Consideremos x=   x1 . . . xd   , W =   w11 . . . w1d . . ..... . . wH1. . . wHd   , b =   b1 . . . bH   , v =   v1 . . . vH   , donde x∈ˆ Des un vector de entrada de la red. Wes la matriz de pesos de los nodos de la capa oculta, es decir, (wi1, . . . , wid) representa el vector de pesos correspondiente al nodo i-´esimo, para todo i= 1, . . . , H. bes el vector de bias de los nodos de la capa oculta, es decir, birepresenta el bias del nodo i-´esimo, para todo i= 1, . . . , H. ves el vector de pesos del nodo de la capa de salida. Observaci´on 3.5. En este caso, los par´ametros variables de la red, que hab´ıamos denotado por p, son W,byv, luego nos referiremos a la salida de la ANN por N(x, W, b, v) y a la ANN por N. En lo sucesivo, para futuras redes que aparezcan en el trabajo seguiremos esta misma notaci´on. 28 Para mejor comprensi´on de la arquitectura escogida, la representamos, por ejemplo, para dos entradas y tres nodos en la capa oculta: Figura 3.1: Topolog´ıa de la red para d= 2 y H= 3. Este diagrama ha sido realizado con Python utilizando las librer´ıas networkx y matplotlib. Volviendo al caso general, se tiene que para todo i= 1, . . . , H, la salida correspondiente al nodo i-´esimo de la capa oculta es σi:= σ(zi), donde zi:= d X j=1 wijxj+bi, luego la salida de la red es N(x, W, b, v) = H X i=1 viσi. Por tanto, para todo k∈N, se tiene ∂kN ∂xk j = H X i=1 vi ∂kσi ∂xk j = H X i=1 vi ∂kσi ∂zk i wk ij. Llamando σ(k) i:= ∂kσi ∂zk i , se comprueba que ∂λ1 ∂xλ1 1 ∂λ2 ∂xλ2 2 . . . ∂λd ∂xλd d N= d X i=1 viPiσ(Λ) i,(3.5) 29 donde Pi= d Y k=1 wλk ik ,Λ = d X i=1 λi. De la ecuaci´on (3.5), obtenemos que la derivada de la red con respecto a cualquiera de sus entradas es equivalente a una ANN con la misma distribuci´on de nodos y capas que N, a la que denotaremos por Ng. Para Ng, los nodos de la capa oculta seguir´an teniendo matriz de pesos Wy vector de bias b. Sin embargo, el nodo de salida tendr´a vector de pesos vP =   v1P1 . . . vdPd   . Adem´as, para todo i= 1, . . . , H, el nodo i-´esimo de la capa oculta tiene como funci´on de activaci´on σ(Λ) i. Teniendo esto en cuenta, podemos calcular las derivadas de Ngcon respecto a los par´ametros de la red original: ∂Ng ∂vi =Piσ(Λ) i,∂Ng ∂bi =viPiσ(Λ+1) i, ∂Ng ∂wij =xjviPiσ(Λ+1) i+viλjwλj−1 ij n Y k=1, k=j wλk ik !σ(Λ) i. Una vez calculadas las derivadas de la red, ya se pueden calcular las de las soluciones de prueba y, por tanto, las derivadas del error cometido con respecto a todos los par´ametros de la red. Con esto, ya estamos en disposici´on de aplicar la t´ecnica de minimizaci´on que se prefiera. Asimismo, en el caso de disponer de la posibilidad de programaci´on en paralelo porque el hardware, las derivadas de la red en cada punto de la malla se pueden obtener de forma simult´anea. 3.2. Primeros experimentos Concluimos este cap´ıtulo resolviendo algunos problemas de EDOs utilizando el m´etodo descrito para as´ı esclarecer los algoritmos resultantes. Asimismo, dichos algoritmos ser´an programados utilizando Python sin uso de librer´ıas espec´ıficas resolutoras de problemas con redes neuronales artificiales, ya que su simpleza nos lo permite. Dados T > 0, t0, a ∈Ryf∈C0([t0, t0+T]×R)∩Lipy([t0, t0+T]×R), vamos a aplicar el m´etodo descrito para resolver el problema de Cauchy y′(t) = f(t, y(t)), t ∈[t0, t0+T] y(t0) = a. 30 para ello tomamos una partici´on de este intervalo por puntos t0, t1, ..., tnde [t0, t0+T] (donde puede que t0>0, tn< T ya que por ejemplo podemos tomar los puntos de Tchebycheff) y buscamos yde la forma y(t) = a+ (t−t0) H X i=1 viσ(bi+wit) (3.6) con σ(t) = 1 1 + e−t,∀t∈R, yv1, .., vH,b1, .., bH,w1, ..., wHa escoger de forma que se minimize J(v, b, w) = m X j=0 y′(tj)−f(tj, y(tj)) 2,(3.7) donde y′(t) = H X i=1 viσ(bi+wit)+(t−t0) H X i=1 viwiσ′(bi+wit). Con esta definici´on, las derivadas de Jcon respecto a los par´ametros de la red siguen como queda: ∂viJ(v, b, w)=2 n X j=0 y′(tj)−f(tj, y(tj)) σ(bi+witj)+(tj−t0)wiσ′(bi+witj)−∂yf(tj, y(tj))(tj−t0)σ(bi+wit) ∂biJ(v, b, w)=2vi n X j=0 y′(tj)−f(tj, y(tj)) σ′(bi+witj)+(tj−t0)wiσ′′(bi+witj)−∂yf(tj, y(tj))(tj−t0)σ′(bi+witj) ∂wiJ(v, b, w) = 2vi n X j=0 tjy′(tj)−f(tj, y(tj)) σ′(bi+witj)+(tj−t0)wiσ′′(bi+witj)−∂yf(tj, y(tj))(tj−t0)σ′(bi+witj). Una vez calculadas las derivadas, podemos proceder a la formulaci´on del algoritmo resultante para buscar yseg´un (3.6) minimizando el error cometido. En las pruebas, tomaremos la tasa de aprendizaje fija, aunque como vimos en el primer cap´ıtulo, pueden escogerse otras opciones: ALGORITMO 1. Fijamos el n´umero de puntos de la partici´on del intervalo de tiempos n∈N, una tasa de aprendizaje inicial η0>0, H∈N, una tolerancia tol > 0 y un n´umero m´aximo de iteraciones K. 31 2. Tomamos tj=t0+jT n,0≤j≤n y definimos J:R3H→Rpor (3.7). 3. Consideramos una inicializaci´on de los pesos y bias w0 i, v0 i, b0 i,1≤i≤H. 4. Mientras k < k0,1 3HH X i=1 |∂wiJk|2+|∂biJk|2+|∂viJk|2> tol2,con Jk=J(vk, bk, wk) tomamos ηk=η0, wk+1 i=wk i−ηk∂wiJk, bk+1 i=bk i−ηk∂biJk, vk+1 i=vk i−ηn∂viJk. 5. Devolvemos la yresultante seg´un (3.6). Para la implementaci´on en Python del algoritmo, hemos necesitado las librer´ıas numpy; tambi´en hemos importado matplotlib para hacer representaciones gr´aficas de los resultados obtenidos. Para importar dichas librer´ıas, incluimos las siguientes instrucciones: from matplotlib import pyplot as plt import numpy as np Esencialmente, el c´odigo est´a estructurado en dos m´etodos: El primero, grad desc, realiza el algoritmo de descenso por el gradiente y nos devuelve los par´ametros de la red calculados seg´un el algoritmo as´ı como un historial del valor del funcional Jen cada iteraci´on. El segundo, calcula grad, act´ua como un m´etodo auxiliar del primero y su funci´on es calcular el funcional Jy sus derivadas Jv= (∂v1J, . . . , ∂vHJ), Jw= (∂w1J, . . . , ∂wHJ) y Jb= (∂b1J, . . . , ∂bHJ) en cada paso. Presentamos el c´odigo del primer m´etodo: def grad_desc (n,eta_0 ,H,tol ,K,a,t_0 ,T,f,f_y ): # Inicializamos los pesos y bias como nulos w=np. zeros ((1 ,H)) b=np. zeros ((1 ,H)) v=np. zeros ((1 ,H)) t=np. linspace (t_0 ,T,n) # Inicializamos contador de iteraciones iter=0 32 # Calculamos el funcional J y sus derivadas #con respecto a los pesos y bias J,J_v ,J_b ,J_w= calcula_grad (H,v,b,w,f,f_y ,t,t_0 ,a) # Calculamos el test de parada e1 =np. linalg . norm (J_w ) **2 e2 =np. linalg . norm (J_b ) **2 e3 =np. linalg . norm (J_v ) **2 err_med =( e1+e2 +e3)/(3* H) # Creamos un historial de errores y almacenamos #el error cometido durante la presente iteracion l=[J [0][0]] while iter <=K and err_med > tol **2: # Actualizacion de la tasa de aprendizaje #(en este caso es fija ) eta = eta_0 # Actualizamos los pesos w=w-eta * J_w b=b-eta * J_b v=v-eta * J_v # Calculamos el funcional J y sus derivadas #con respecto a los pesos y bias J,J_v ,J_b ,J_w= calcula_grad (H,v,b,w,f,f_y ,t,t_0 ,a) # Actualiamos el test de parada e1 =np. linalg . norm (J_w ) **2 e2 =np. linalg . norm (J_b ) **2 e3 =np. linalg . norm (J_v ) **2 err_med =( e1+e2 +e3)/(3* H) # Incrementamos el contador de iteraciones iter+=1 # Almacenamos el error cometido #en la presente iteracion l. append (J [0][0]) return w,b,v,l En el caso del segundo m´etodo, dado que el valor del sigmoide y sus derivadas de primer y segundo orden evaluados en wi+bitjpara todos i= 1, . . . , H,j= 0, . . . , n se utilizan varias veces, hemos almacenado estos valores como variables en cada instante de la partici´on. Para ello, nos ha sido de utilidad definir los m´etodos: 33 def sigmoide (t): return 1/(1+ np.exp (-t)) def der_sigmoide (t): return sigmoide (t) *(1 - sigmoide (t)) def der_seg_sigmoide(t): return der_sigmoide (t) *(1 -2* sigmoide (t)) El c´odigo del segundo m´etodo es: def calcula_grad (H,v,b,w,f,f_y ,t,t_0 ,a): # Inicializamos a cero el funcional J y #sus derivadas para hacer una suma iterada t J_w = np. zeros ((1 , H)) J_b = np. zeros ((1 , H)) J_v = np. zeros ((1 , H)) J=0 for tj in t: # Almacenamos el sigmoide y sus derivadas #para ahorrar calculos sigm = sigmoide (b+w* tj) sigmder = der_sigmoide (b+w*tj ) sigmdersec = der_seg_sigmoide (b+w*tj) dift =tj -t_0 # Calculamos y ’(tj) e y (tj) yder_j =np.dot (v ,np. transpose ( sigm ))+ dift *np.dot (np. multiply (v,w) ,np. transpose ( sigmder )) y_j =a+ dift*np.dot (v ,np. transpose ( sigm )) fj=f(tj ,y_j) dif=yder_j -fj f_yj =f_y (tj ,y_j ) # Actualizamos el funcional J y sus derivadas J=J+dif **2 J_v =J_v +dif *( sigm +dift *np . multiply (w ,sigmder )-f_yj * dift*sigm) term =dif *( sigmder + dift *np . multiply (w, sigmdersec ) - f_yj * dift * sigmder ) J_b = J_b + term J_w = J_w +tj*term # Multiplicamos por v el resultado de las sumas anteriores J_b =np. multiply (v,J_b ) J_w =np. multiply (v,J_w ) 34 return J,J_v ,J_b ,J_w Una vez implementado el algoritmo, podemos obtener el valor de ypara los par´ametros obtenidos mediante: def y(w,b,v,a,t_0 ,ti): sigm = sigmoide (b+w* ti) return a+(ti - t_0 )*np.dot (v ,np. transpose ( sigm )) Comenzamos resolviendo un problema muy simple y′(t) = y(t), t ∈[0,1] y(0) = 1.(3.8) Para ello, vamos a tomar H=3 a=1 n=20 eta_0 =0.1 tol =0.0001 K=100 t_0 =0 T=1 def f(t,y): return y def f_y (t,y): return 1 Lo ejecutamos con la instrucci´on w,b,v ,err = grad_desc (n,eta_0 ,H,tol ,K,a ,t_0 ,T,f, f_y ) Comenzamos representando el funcional Ja lo largo de las iteraciones realizadas plt . plot ( np. array ( range(len(err ))),err) 35 De forma similar a la prueba de 1, tenemos que PN k=1 wkxk+b+es combinaci´on lineal de una constante, ρ ai+r0 2rb+r0 2r N X k=1 wkxk!yρ ai−r0 2+r0 2rb+r0 2r N X k=1 wkxk! suponiendo N X k=1 wkxk+b∈[−r, r].(4.3) Como Ω es acotado, podemos tomar ren cada neurona lo suficientemente grande para que se satisfaga (4.3). Finalmente, reemplazando las neuronas y razonando como en el caso anterior, obtenemos una red Ncon 2Uneuronas, profundidad Ly 4Wpesos. Esto concluye la prueba. Por tanto, queda justificado el limitarnos a considerar las ReLU en lugar de otras funciones lineales a trozos. Por simplificar notaci´on, a veces nos referiremos por ANN ReLU a aquellas redes con funciones de activaci´on ReLU. Adem´as, hay que comentar que dada una funci´on φ: [0,1]d→Ry su aproximaci´on por una ANN ˆφ, definiremos el error de aproximaci´on en norma infinito, i.e. ∥φ−ˆφ∥∞= max x∈[0,1]d|φ(x)−ˆφ(x)|. Una vez definido el error, presentamos el siguiente resultado: Proposici´on 4.2. La funci´on φ(x) = x2en el intervalo [0,1] puede ser aproximada con cualquier error ε > 0 por una red neuronal con funci´on de activaci´on ReLU que se anula en 0 y con profundidad, n´umero de pesos y neuronas de orden O(log(1/ε)). Demostraci´on. Consideremos la funci´on diente g: [0,1] →[0,1] definida como sigue g(x) :=      2x, si x < 1 2 2(1 −x),si x < 1 2 Figura 4.1: Funci´on diente. y para todo s∈N, definimos la funci´on sierra por gs(x) := g◦g◦. . . ◦g s (x). 42 Figura 4.2: Funciones sierra para s= 2,3,4. Si gses una funci´on sierra con 2s−1dientes uniformemente distribuidos, se tiene gs(x) :=          2sx−2k 2s,si x∈2k 2s,2k+ 1 2s, k = 0,1,...,2s−1−1, 2s2k 2s−x,si x∈2k−1 2s,2k 2s, k = 1,2,...,2s−1. Este resultado puede encontrarse en [17]. Vamos a ver que φ(x) = x2puede aproximarse con combinaciones lineales de las funciones gs. Para ello, vamos a considerar como φmel spline lineal a trozos que coincide con φen los puntos k/2m,k= 0,...,2m. Ver figura 4.3 para el caso m= 4: Figura 4.3: Recordando la f´ormula de error para el spline lineal, se tiene εm=∥φ−φm∥∞≤∥φ′′∥∞ 81 2m2 =1 22m+2 . 43 Adem´as, observemos que refinar la interpolaci´on de φm−1aφmequivale a ajustarla por una funci´on proporcional a la funci´on sierra: φm−1(x)−φm(x) = gm(x) 22m, de donde aplicando esta igualdad de forma recursiva y observando que g0(x) = x, nos queda φm(x) = x− m X s=1 gs(x) 22s.(4.4) Por otra parte, g(x) se puede expresar como una ANN ReLU ya que se verifica g(x)=2x+−4s−1 2+ + 2(x−1)+. Adem´as, gracias a (4.4), tenemos que la construcci´on de φmrequiere operaciones lineales y composiciones de gdel orden O(m). Por tanto, φmtambi´en puede implementarse usando una ANN con pesos, profundidad y neuronas del orden O(m), lo que concluye la prueba. A partir de esta proposici´on, se tiene el siguiente corolario el cual nos va a permitir aproximar la funci´on producto con precisi´on: Corolario 4.3. Dado M > 0 y ε∈(0,1), existe una red neuronal artificial con funciones de activaci´on ReLU ηcon dos entradas que devuelve una funci´on ˆφ:R2→ Rtal que 1. Para cualquier entrada (x, y)∈R2tal que |x| ≤ My|y| ≤ M, se tiene |ˆφ(x, y)−xy| ≤ ε. 2. Si x= 0 o y= 0, entonces ˆφ(x, y) = 0. 3. El n´umero de capas, el n´umero de pesos y el n´umero de neuronas de ηson como mucho c1log(1/ε) + c2, con c1una constante universal positiva y c2una constante positiva que depende de M. Demostraci´on. En primer lugar, observemos que se tiene xy =1 2((x+y)2−x2−y2).(4.5) Adem´as, tomemos ˆgsq,δ como la funci´on dada por la Proposici´on 4.2 con εreemplazado por δpara δ > 0. A continuaci´on, sin p´erdida de generalidad, suponemos M≥1 y definimos ˆφ(x, y) = M2 8ˆgsq,δ |x+y| 2M−ˆgsq,δ |x| 2M−ˆgsq,δ |y| 2M, con δ= 8ε/3M2. Vamos a comprobar que la ˆφdefinida verifica el enunciado. Claramente, se verifica la afirmaci´on 2 y, usando (4.5), se tiene la afirmaci´on 1. Por su parte, la afirmaci´on 3 se obtiene aplicando la Proposici´on 4.2 y recordando el δ escogido. 44 A continuaci´on, vamos a estudiar el error cometido al aproximar funciones m´as generales. En particular, dado que estamos interesados en aproximar las soluciones de problemas diferenciales, nos va a interesar situarnos en los espacios de Sobolev Wr,∞((0,1)d) = nu∈L∞((0,1)d) : Diu∈L∞((0,1)d)di,∀i≤ro, d, r ∈Z+, dotados de su norma usual, i.e. dada u∈Wr,∞((0,1)d), se verifica ∥u∥Wr,∞((0,1)d)=∥u∥L∞((0,1)d)+ r X i=1 ∥Diu∥L∞((0,1)d)di. En lo sucesivo, denotaremos por Fd,r a la bola unidad en Wr,∞((0,1)d), es decir: Fd,r =u∈Wr,∞((0,1)d) : ∥u∥Wr,∞((0,1)d)≤1.(4.6) Asimismo, introducimos la siguiente definici´on que usaremos m´as adelante: Definici´on 4.4. Diremos que la topolog´ıa o arquitectura de una red neuronal artificial es capaz de expresar cualquier funci´on de Fd,r con error ε > 0si existe una asignaci´on de pesos y b´ıas para la cual esto es cierto. Esto nos lleva a enunciar el teorema principal de esta secci´on. Dicho resultado nos aporta cotas superiores sobre el n´umero de capas, neuronas y pesos necesarios para aproximar funciones de Fd,r mediante ANN con funciones de activaci´on ReLU. Teorema 4.5. Para todos d, r ∈Z+yε∈(0,1), existe una arquitectura de ANN con funciones de activaci´on ReLU y una constante c=c(d, r)verificando: 1. Es capaz de expresar cualquier funci´on de Fd,r con error ε. 2. Como m´aximo tiene n´umero de capas c(log(1/ε) + 1). 3. El n´umero de pesos y neuronas es inferior o igual a cε−d/r(log(1/ε) + 1). Demostraci´on. La prueba de este resultado se divide en dos partes: Primera parte: En primer lugar, aproximaremos φ∈Fd,r por una combinaci´on de sumas y productos φ1de polinomios de Taylor locales y funciones lineales a trozos unidimensionales. Sea Nun entero positivo. Consideremos la partici´on de la unidad en [0,1]d formada por las (N+ 1)dfunciones ϕmdefinidas por ϕm(x) = d Y k=1 ψ3Nxk−mk N,(4.7) 45 donde m= (m1, . . . , md)∈ {0,1, . . . , N}dy ψ(y) =    1,si |y|<1, 2− |y|,si 1 ≤ |y| ≤ 2, 0,si 2 <|y|. Observemos que ∥ψ∥∞= 1,∥ϕm∥∞= 1,∀m(4.8) y sop ϕm⊂x:xk−mk N<1 N,∀k.(4.9) Para cualquier m∈ {0, . . . , N}d, consideramos el polinomio de Taylor de grado r−1 de la funci´on φen x=m/N: Pm(x) = X n∈{0,...,N}d |n|<r Dnφ n!x=m Nx−m Nn,(4.10) siendo n! = Qd k=1 nk! y (x−m/N)n=Qd k=1(xk−mk/N)nk. Aproximamos φpor φ1=X m∈{0,...,N}d ϕmPm.(4.11) Vamos a obtener una cota del error de aproximaci´on cometido. Por (4.8)y(4.9), se tiene |φ(x)−φ1(x)|=X m ϕm(x)(φ(x)−Pm(x)) ≤X |x−m N|∞<1 N |φ(x)−Pm(x)| Adem´as, observando que todo x∈[0,1]dpertenece al soporte de a lo sumo 2d funciones ϕm, tenemos |φ(x)−φ1(x)| ≤ 2dmax |x−m N|∞<1 N |φ(x)−Pm(x)| ≤2ddr r!1 Nr ∥Dnφ(x)∥Wr,∞((0,1)d)≤2ddr r!1 Nr . Por tanto, si tomamos N=&r! 2ddr ε 2−1/r',(4.12) (⌈·⌉ denota la funci´on techo o parte entera por exceso) entonces ∥φ−φ1∥∞≤ε 2.(4.13) Adem´as, por (4.10), los coeficientes de los polinomios Pmest´an uniformemente acotados para todo φ∈Fd,r: Pm(x) = X n:|n|<r am,n x−m Nn,|am,n| ≤ 1.(4.14) 46 Por tanto, la prueba se reduce a construir una arquitectura de ANN capaz de aproximar con error ε/2 cualquier funci´on de la forma (4.11) suponiendo que Nviene dado por (4.12) y que los polinomios Pmson de la forma (4.14). Este es el objetivo de la siguiente parte de la prueba. Segunda parte: Sustituyendo (4.14) en (4.11), tenemos: φ1(x) = X m∈{0,...,N}dX n:|n|<r am,nϕm(x)x−m Nn.(4.15) La expresi´on obtenida es una combinaci´on lineal de no m´as de dn(N+ 1)dt´erminos ϕm(x)(x−m N)n. A su vez, cada uno de estos t´erminos es el producto de a lo sumo d+r−1 funciones lineales a trozos de una variable: dfunciones ψ(3Nxk−3mk) y como mucho n−1 expresiones lineales de la forma xk−mk N. A continuaci´on vamos a aproximar dicho producto utilizando una ANN vali´endonos del Corolario 4.3. As´ı, sea ˆ ×la aproximaci´on del producto obtenida por dicho corolario para M=n+dy un margen de error δa escoger posteriormente. Consideremos la aproximaci´on del producto ϕm(x)(x−m N)ndefinida por la composici´on: ˆφm,n(x) = ˆ ×(ψ(3Nx1−3m1),ˆ ×(ψ(3Nx2−3m2),...,ˆ ×(xk−mk N, . . .). . .)).(4.16) Por la tercera tesis del Corolario 4.3, sabemos que ˆφm,n puede aproximarse por una red neuronal artificial de funciones de activaci´on ReLU con profundidad, n´umero de neuronas y n´umero de pesos no superior a (d+r)c1log(1/δ) para una constante c1 que depende de δyr. Vamos a estimar el error cometido por esta aproximaci´on. Observemos que se tiene |ψ(3Nk−3mk)| ≤ 1,|xk−mk N| ≤ 1,∀k= 1, . . . , d, ∀x∈[0,1]d Por la primera tesis del Corolario 4.3, para a, b ∈Rcon |a| ≤ 1 y |b| ≤ M, se verifica |ˆ ×(a, b)| ≤ |b|+δ. Luego suponiendo δ < 1 y observando que los argumentos de todos los productos aproximados de (4.16) est´an acotados por el Mescogido, tenemos que dicha tesis se verifica para cada uno de dichos productos. Por tanto, se tiene ˆφm,n(x)−ϕm(x)x−m Nn =ˆ ×(ψ(3Nx1−3m1),ˆ ×(ψ(3Nx2−3m2),ˆ ×(ψ(3Nx3−3m3), . . .))) −ψ(3Nx1−3m1)ψ(3Nx2−3m2)ψ(3Nx3−3m3). . .  ≤ˆ ×(ψ(3Nx1−3m1),ˆ ×(ψ(3Nx2−3m2),ˆ ×(ψ(3Nx3−3m3), . . .))) −ψ(3Nx1−3m1)ˆ ×(ψ(3Nx2−3m2),ˆ ×(ψ(3Nx3−3m3), . . .)) 47 +ψ(3Nx1−3m1)ˆ ×(ψ(3Nx2−3m2),ˆ ×(ψ(3Nx3−3m3), . . .)) −ψ((3Nx2−3m2)ˆ ×(ψ(3Nx3−3m3), . . .)+. . . ≤(d+r)δ. Luego hemos llegado a la cota ˆφm,n(x)−ϕm(x)x−m Nn≤(d+r)δ. (4.17) Adem´as, de la segunda afirmaci´on del Corolario 4.3 deducimos ˆφm,n(x) = ϕm(x)x−m Nn, x /∈sop ϕm.(4.18) Ahora definimos la aproximaci´on completa por ˆφ=X m∈{0,...,N}dX n:|n|<r am,n ˆφm,n.(4.19) Estimamos el error cometido al aproximar φ1por ˆφ: ˆφ(x)−φ1(x)= (4,15)X m∈{0,...,N}dX n:|n|<r am,nˆφm,n(x)−ϕm(x)x−m Nn = (4,18)X m:x∈sop ϕmX n:|n|<r am,nˆφm,n(x)−ϕm(x)x−m Nn . Como |an,m| ≤ 1 y x∈sop ϕmpara a lo sumo 2dfunciones ϕm, nos queda: ˆφ(x)−φ1(x)≤2dmax m:x∈sop ϕmX n:|n|<r ˆφm,n(x)−ϕm(x)x−m Nn≤ (4,17) 2ddr(d+r)δ. Luego si tomamos δ=ε 2d+1dr(d+r),(4.20) entonces ∥ˆφ−φ1∥∞≤ε 2y, por (4.13), llegamos a ∥ˆφ−φ∥∞≤ ∥ ˆφ−φ1∥∞+∥φ1−φ∥∞≤ε 2+ε 2≤ε. Esto prueba la primera tesis del teorema. Por otra parte, observemos que por (4.19), ˆφpuede implementarse por una ANN con varias subredes que calculen cada ˆφm,n en paralelo. La salida de la ANN ser´ıa la suma ponderada de las salidas de las subredes por los pesos an,m. La arquitectura de la ANN no depende de la φque estemos aproximando, solo los an,m. Como hemos mencionado anteriormente, cada subred no tiene m´as de c1log(1/δ) capas, pesos y neuronas. Adem´as, no hay m´as de dr(N+1)dde dichas subredes. Por tanto, la red completa no tiene m´as de c1log(1/δ)+1 capas y dr(N+1)d(c1log(1/δ)) pesos y neuronas. Tomando δdado por (4.20) y Npor (4.12), queda probada la segunda tesis del teorema. Esto concluye la prueba. 48 Como hemos visto, este resultado nos aporta cotas superiores para el n´umero de pesos, capas y neuronas necesarios para conseguir redes neuronales con una topolog´ıa capaz de aproximar los elementos de Fd,r con error ε∈(0,1). Sin embargo, en lugar de establecer una arquitectura fija y entrenar la red hasta obtener los pesos y b´ıas que minimicen el error como hemos hecho hasta ahora, podemos adem´as adaptar la arquitectura de la red a la funci´on aproximada. De este modo, cabr´a esperar que podamos encontrar ANN menos complejas (con menos capas y neuronas) que aproximen con el mismo error que las dadas por el Teorema 4.5 a cambio de encontrar la arquitectura apropiada. Asimismo, por el Teorema 4.5, el n´umero de capas est´a acotado superiormente por una cantidad del orden log(1/ε) y el n´umero de pesos y neuronas por un n´umero del orden log(1/ε)/ε para ε∈(0,1). Sin embargo, estas cotas pueden ser muy grandes. As´ı, por ejemplo, en el caso del conjunto F1,1(ver (4.6)), ha sido probado en [12] que estas cotas pueden mejorarse del siguiente modo: Proposici´on 4.6. Existe una constante c > 0 tal que para toda φ∈W1,∞(0,1) con ∥φ∥W1,∞(0,1) ≤1 y todo ε∈(0,1/2), existe una ANN ηde profundidad 6 con funciones de activaci´on ReLU capaz de aproximar φcon error εusando a lo sumo c/(εlog(1/ε)) pesos y neuronas. 4.2. Cotas inferiores sobre el n´umero de elementos de una ANN para aproximar funciones Una vez establecidas cotas superiores sobre el n´umero de capas, pesos y neuronas necesarios para aproximar elementos de Fd,r con error ε∈(0,1), vamos a obtener cotas inferiores en esta secci´on. Para ello, comenzamos enunciando el siguiente teorema probado en [14]: Teorema 4.7. Sean d, r yε > 0. Consideremos Wun entero positivo y una aplicaci´on η:RW→C([0,1]d). Supongamos que existe una aplicaci´on continua M:Fd,r →RWtal que ∥φ−η(M(φ))∥∞≤εpara todo φ∈Fd,r. Entonces se tiene W≥cε−d/r para una constante c=c(r). Vamos a aplicar este teorema tomando ηcomo una ANN ReLU igual que hemos hecho hasta ahora y RWser´a el espacio de pesos de la red. Si suponemos que la red es capaz de aproximar cualquier funci´on de Fd,r con error εescogiendo los pesos correspondientes de forma continua, del teorema se tiene que la red tendr´a al menos cε−d/r pesos. Adem´as, si denotamos por Kal n´umero de conexiones y por U el n´umero de neuronas de la ANN, para que la red est´e bien definida debe verificarse W≤K+U U≤K⇒W≤2K⇒K≥W 2≥cε−d/r 2, luego hemos acotado tambi´en el n´umero de conexiones. 49 Observaci´on 4.8. La selecci´on continua de los pesos es crucial en el Teorema 4.7. De este modo, examinando la prueba de la cota superior obtenida en el Teorema 4.5, observamos que los pesos fueron escogidos de forma continua y se obtuvo que el n´umero de pesos estaba acotado superiormente por una cantidad de O(ε−d/r log(1/ε). Esta cota est´a por encima de la cota inferior obtenida mediante el Teorema 4.7, luego ambos resultados son coherentes. Sin embargo, en el Teorema 4.6, al ser la arquitectura dependiente de la φa aproximar, los pesos no se toman de forma continua. De hecho, observamos que dicho Teorema acota el n´umero de pesos por c/ε log(1/ε), que es menor que la cota dada por 4.7. En consecuencia de la observaci´on anterior, hay que estudiar el caso en el que la dependencia de los pesos seg´un el valor de φno sea continua necesariamente. Para ello, vamos a dar cotas inferiores sobre la arquitectura de la red utilizando la dimen- si´on de Vapnik-Chervonenkis o dimensi´on VC, que es una medida de la complejidad que tiene una familia de funciones capaz de clasificar de forma binaria nubes de puntos de Rd, es decir, de separar dichos puntos en dos conjuntos en funci´on de sus salidas. Introducimos este concepto formalmente mediante las siguientes definiciones dadas en [15]y[16]: Definici´on 4.9. Sea X⊂RdyHuna familia de funciones de Rd⊂Xen R. Al par (X, H)lo denominaremos espacio de aprendizaje. Definici´on 4.10. Dado un espacio de aprendizaje (X, H), diremos que una nube de puntos S⊂Xes desmenuzada por Hsi todo subconjunto de Xes aprendido por H, es decir, si para todo subconjunto A⊂Sexiste h∈Htal que para todo punto de Sse verifique x∈A⇔h(x)>0. Definici´on 4.11. La funci´on crecimiento de un espacio de aprendizaje (X, H) que denotaremos por ΠH:N→N, se define como ΠH(m) = max {|H|S|:S⊂Xcon |S|=m}. Definici´on 4.12. Se define la dimensi´on de Vapnik-Chervonenkis o dimensi´on VC de un espacio de aprendizaje (X, H)como dimV C(X, H) = sup{n∈N| ∃S⊂Xcon |S|=ndesmenuzado por H}. Observaci´on 4.13. Dado S⊂Xcon |S|=m, el cardinal del conjunto de las partes de Sser´a 2m, por lo que claramente se tiene ΠH(m)≤2m. Adem´as, podemos definir la dimensi´on V C de forma equivalente como el mayor mtal que ΠH(m)=2m. Una vez definida la dimensi´on VC, enunciamos dos teoremas cuya prueba est´a en [15] sobre dicha dimensi´on para las clases de funciones aproximadas mediante 50 ANN con funciones de activaci´on polinomiales a trozos en las capas ocultas y una funci´on lineal a trozos en la capa de salida. Dado que el espacio X⊂Rdser´a el conjunto de puntos de entrada de entrenamiento de la red en cuesti´on, en lo sucesivo denotaremos la dimensi´on VC por dimV C(H) para simplificar la notaci´on. Teorema 4.14. Supongamos una red neuronal artificial con una ´unica salida que consta de Wpesos y Uneuronas tal que la funci´on de activaci´on en la capa de salida es lineal a trozos y en el resto de capas son funciones polinomiales con p trozos y grado inferior a l. Entonces, si Hes la clase de funciones capaces de ser aproximadas por dicha red se tiene dimV C(H) = O(W(W+Ul log2p)). Teorema 4.15. Consideremos una ANN de la forma descrita en el Teorema 4.14. Adem´as, supongamos que la profundidad de dicha red es L. Entonces, si Hes la clase de funciones capaces de ser aproximadas por dicha red, se tiene ΠH(m)≤2L(e Ump(l+ 1)L−1)W L,∀m∈N y dimV C(H)≤2WL log24WLUp log 2 + 2WL2log2(l+ 1) + 2L. As´ı, para pylfijos, se verifica dimV C(H) = O(WL log2W+WL2). Utilizando estos resultados, vamos a probar un teorema de cotas inferiores sobre el n´umero de pesos necesarios para aproximar funciones de Fd,r mediante ANN con funciones de activaci´on ReLU. Teorema 4.16. Se tiene: 1. Para cualquier ε∈(0,1), una arquitectura red neuronal artificial con funciones de activaci´on ReLU capaz de aproximar cualquier funci´on f∈Fd,r con error εtiene que tener al menos cε−d/(2r)pesos para una constante c=c(d, r)>0. 2. Sean p≥0yc1>0contantes. Para cualquier ε∈(0,1/2), si una arquitectura de una ANN con funciones de activaci´on ReLU con profundidad L≤c1logp(1/ε)es capaz de aproximar cualquier funci´on f∈Fd,r con error ε, entonces la red debe tener al menos c2ε−d/r pesos para una constante c2= c2(d, r, p, c1). Demostraci´on. Estamos interesados en aplicar los resultados anteriores al caso en el que Hrepresenta la familia de funciones obtenidas aplicando funciones caracter´ıstica de la forma χ{x>a}a una ANN ReLU con arquitectura fija pero pesos variables. En este caso, el Teorema 4.14 nos da dimV C(H)≤c3W2(4.21) 51 En consecuencia, los resultados vistos en este cap´ıtulo nos muestran que las ANN ReLU profundas son m´as eficientes para aproximar funciones derivables que las redes superficiales. As´ı, el Teorema 4.5 hemos comprobado que las funciones del espacio de Sobolev Wr,∞((0,1)d) pueden aproximarse con error ε > 0 por redes ReLU con profundidad de O(log(1/ε) y un n´umero de neuronas de O(ε−d/r log(1/ε). En cambio, por el Teorema 4.24, una funci´on no lineal de clase C2([0,1]d) no puede ser aproximada con error εpor una red ReLU de profundidad Lcon un n´umero de neuronas inferior a cε−1/(2(L−2)) (con cuna constante). En particular, esto implica que, en cuanto al n´umero de neuronas requeridas para aproximar con error εfunciones de Wr,∞((0,1)d), las redes sin cota sobre su profundidad son m´as eficientes que las que tienen una profundidad fija Lal menos en el caso d r<1 2(L−2), que se tendr´a para Lpeque˜no, i.e. redes neuronales superficiales. Para funciones de clase infinito como los polinomios, la eficiencia de la profundidad a´un est´a m´as acentuada, ya que pueden aproximarse con ANN con un n´umero de neuronas de O(log(1/ε)) gracias a la prueba de 4.5. 4.4. Test de validaci´on de las cotas de errores Las cotas superiores de error que hemos obtenido en la secci´on 4.1 est´an basadas en la Proposici´on 4.2, donde se prueba que podemos aproximar la funci´on y=x2en [0,1] con un error menor que ε > 0 usando, a lo sumo, un n´umero de pesos, capas y neuronas de orden log(1/ε). Sin embargo, hay que recordar que en la pr´actica esta aproximaci´on se va a obtener utilizando un algoritmo de descenso el cual podr´ıa converger hacia aproximaciones peores (no podemos asegurar convergencia hacia el m´ınimo global del problema). En esta secci´on, vamos a realizar algunos tests para comprobar si el algoritmo de descenso es, efectivamente, capaz de encontrar la aproximaci´on deseada. En los c´alculos num´ericos se ha utilizado la librer´ıa de Python TensorFlow, cuyo funcionamiento expondremos m´as adelante en el Cap´ıtulo 5. Tomando ε= 10−l, con l= 1,2,3,4,5, hemos resuelto con TensorFlow el problema m´ın p  X x∈ˆ Ωˆφp(x)−x22   , donde ˆ Ω es un subconjunto discreto de [0,1] formado por 1000 puntos uniformemente distribuidos y ˆφpuna red neuronal con pesos y b´ıas dados por py n´umero de pesos, 58 capas y neuronas de orden O(log(1/ε)) = O(l). De hecho, en nuestro caso, hemos tomado l+ 1 capas (locultas m´as la de salida) y 5l+ 1 neuronas (5 neuronas por capa oculta m´as la neurona de la capa de salida). En TensorFlow, la elecci´on de pesos y b´ıas inicial es aleatoria, de forma que al aplicar el mismo algoritmo varias veces no se obtiene siempre la misma aproximaci´on. Debido a ello, para cada lhemos ejecutado el algoritmo 50 veces. En cada ejecuci´on del algoritmo, se han realizado 1500 iteraciones. Para ver si el m´etodo es capaz de darnos la aproximaci´on de la Proposici´on 4.2, hemos incluido un histograma con la frecuencia de errores, que en particular nos va a decir en qu´e porcentaje de ejecuciones obtenemos un error inferior o igual al deseado. Adem´as, indicamos el error obtenido en la mejor ejecuci´on junto con el error medio cometido y la desviaci´on t´ıpica de este entre las 50 ejecuciones. Veamos el resultado obtenido tras ejecutar para los distintos valores de l: Para l= 1, i.e. ε= 0′1, se obtiene: Mejor iteracion: 14 con un error de 5.2365947e-05 El error promedio entre los 50 tests realizados es 0.019358205 La desviaci´ on t´ ıpica del error entre los 50 tests realizados es 0.034907244 Figura 4.4: Histograma de rango de errores en las 50 ejecuciones para l= 1. Para l= 2, i.e. ε= 0′01, se obtiene: Mejor iteracion: 33 con un error de 1.8201674e-05 El error promedio entre los 50 tests realizados es 0.020180117 La desviaci´ on t´ ıpica del error entre los 50 tests realizados es 0.03660553 59 Figura 4.5: Histograma de rango de errores en las 50 ejecuciones para l= 2. Para l= 3, i.e. ε= 0′001, se obtiene: Mejor iteracion: 14 con un error de 9.662816e-06 El error promedio entre los 50 tests realizados es 0.016420849 La desviaci´ on t´ ıpica del error entre los 50 tests realizados es 0.03404927 Figura 4.6: Histograma de rango de errores en las 50 ejecuciones para l= 3. Para l= 4, i.e. ε= 10−4, se obtiene: Mejor iteracion: 21 con un error de 9.730604e-06 El error promedio entre los 50 tests realizados es 0.018242517 La desviaci´ on t´ ıpica del error entre los 50 tests realizados es 0.035436343 60 Figura 4.7: Histograma de rango de errores en las 50 ejecuciones para l= 4. Para l= 5, i.e. ε= 10−5, se obtiene: Mejor iteracion: 48 con un error de 4.789665e-06 El error promedio entre los 50 tests realizados es 0.01786309 La desviaci´ on t´ ıpica del error entre los 50 tests realizados es 0.035607465 Figura 4.8: Histograma de rango de errores en las 50 ejecuciones para l= 5. En vista de los resultados obtenidos, la desviaci´on t´ıpica siempre es elevada, lo que implica una gran dependencia de la inicializaci´on de pesos y b´ıas aleatoria realizada por el algoritmo. Esto explica que aunque en los histogramas se observe una mayor proporci´on de ejecuciones en las que se consigue el error deseado, el promedio suele estar por encima de dicho error. Adem´as, conforme aumenta el valor de l, se observa que la mejor aproximaci´on obtenida es mejor que la esperada, mostrando que las cotas superiores consideradas, 61 probablemente, no son ´optimas, tal y como se coment´o en la secci´on 4.1. As´ı, la mejor iteraci´on en cada caso nos aporta una cota superior sobre el error cometido. Por ´ultimo, hay que comentar que el aumento de ltambi´en conlleva un aumento en la proporci´on de ejecuciones en las que se obtienen errores m´as peque˜nos. Sin embargo, los errores estrictamente inferiores a 10−5se consiguen en un n´umero reducido de ejecuciones incluso para l= 5. 62 5. Algunos modelos implementados con Python En el Cap´ıtulo 3, estudiamos un m´etodo para resolver problemas diferenciales mediante redes neuronales artificiales minimizando el error cuadr´atico cometido al aproximar la soluci´on de la ecuaci´on diferencial correspondiente. A continuaci´on, en este cap´ıtulo vamos a resolver diversos problemas de inter´es mediante una versi´on mejorada de dicho m´etodo utilizando la librer´ıa Tensorflow de Python. 5.1. Esquema de resoluci´on de problemas de optimizaci´on con ANN. Optimizadores En primer lugar, mostramos el esquema de entrenamiento de una ANN para resolver problemas de optimizaci´on mediante el siguiente diagrama de flujo: Figura 5.1: Entrenamiento de una ANN para resolver problemas de optimizaci´on. En el diagrama introducimos dos t´erminos usualmente utilizados en problemas de optimizaci´on: funci´on de p´erdida o residuo y optimizador. Dichos t´erminos representan el funcional a minimizar (o error) y el m´etodo de descenso utilizado, respectivamente. De este modo, como ya hemos visto en los cap´ıtulos anteriores, en cada iteraci´on la red neuronal recibe como entrada los puntos de la discretizaci´on del dominio 63 realizada junto con los par´ametros variables (pesos y b´ıas), y nos da una salida o predicci´on asociada a dichos par´ametros. A continuaci´on, calculamos el error cuadr´atico medio cometido y, utilizando un m´etodo de descenso por gradiente, obtenemos una nueva distribuci´on de pesos y b´ıas que aporten una mejor aproximaci´on de la soluci´on. Sin embargo, vimos que no podemos asegurar convergencia del m´etodo hacia el m´ınimo global de la funci´on de p´erdida, solo hacia m´ınimos locales. Adem´as, el algoritmo de descenso por gradiente cl´asico tiene una convergencia muy lenta y tiene altos costes computacionales. Para tratar de evitar los m´ınimos locales (aunque no podemos asegurarlo) y aumentar el rendimiento del algoritmo, se introducen algunas variantes del algoritmo de la secci´on 3.2, donde recordemos que la actualizaci´on de un par´ametro w(peso o b´ıas) en una etapa kera la siguiente: wk+1 i=wk i−ηk∂wiJk, siendo ηkla tasa de aprendizaje y Jkla funci´on de p´erdida en la etapa k. Presentamos las siguientes variantes en cuanto a la actualizaci´on de los par´ametros en cada etapa: Momentum. Consiste en modificar la actualizaci´on de pesos y b´ıas del algoritmo considerando no ´unicamente el gradiente de la funci´on de p´erdida en el paso actual, sino los anteriores. As´ı, para un par´ametro (peso o b´ıas) wen una etapa k, su actualizaci´on queda como sigue: wk+1 =wk−ηk∂wJk+m∂wJk, donde, al igual que antes ηkes la tasa de aprendizaje y Jkla funci´on de p´erdida. Adem´as, ∂wJkrepresenta el promedio de los gradientes anteriores y m∈[0,1] es un coeficiente a escoger. Stochastic Gradient Descent (SGD). Como comentamos en el Cap´ıtulo 1, las librer´ıas que implementan algoritmos de descenso con ANN tienen implementado como SGD a lo que previamente nos referimos como Mini-Batch. Este m´etodo es utilizado frecuentemente combinado junto con el m´etodo del momentum. Adaptive Gradient Algorithm (Adagrad). Es un optimizador que adapta la tasa de aprendizaje para cada par´ametro individual en funci´on de la suma acumulada de los gradientes cuadrados anteriores. Para ello, introducimos acumuladores de los gradientes de cada par´ametro w(peso o b´ıas) en la etapa k del siguiente modo: G0 w= 0, Gk w=Gk−1 w+∂wJk·(∂wJk)T A partir de Gk w, la actualizaci´on de wen la etapa k, queda como sigue: wk+1 =wk−η pGk w+ε∂wJk. As´ı, en pendientes elevadas, se reduce la tasa ηmientras que en pendientes reducidos esta aumenta. 64 Adadelta. La acumulaci´on continua de gradientes, hace que crezca hasta el punto de que la tasa de aprendizaje sea tan peque˜na que el entrenamiento frene. Para paliar esto, el m´etodo Adadelta surge como una modificaci´on del Adagrad en la que ´unicamente se tienen en cuenta en los acumuladores los n ´ultimos gradientes para nprefijado. Root Mean Square Propagation (RMSProp). Este optimizador es similar a AdaDelta. La diferencia est´a en la forma de asignar el factor de entrenamiento a cada par´ametro interno, de forma que se le da mayor prioridad a los ´ultimos valores del gradiente a˜nadidos al promedio que a los m´as antiguos. Luego el esquema es el siguiente G0 w= 0, Gk w=βGk−1 w+ (1 −β)∂wJk·(∂wJk)T,con β∈[0,1], wk+1 =wk−η pGk w+ε∂wJk. Adaptive moment estimation (Adam). El optimizador Adam combina los m´etodos AdaGrad y RMSProp junto con el momentum, usando as´ı tanto el promedio del hist´orico de gradientes como los acumuladores de gradientes. Este ser´a el escogido para los tests desarrollados en la presente memoria. Una vez afianzado el esquema de resoluci´on para problemas de optimizaci´on, as´ı como presentados los principales optimizadores implementados en las librer´ıas, nos quedan por comentar, para el caso de problemas diferenciales, dos modificaciones al algoritmo de la secci´on 3.2. No vamos a separar las restricciones del problema de la salida red neuronal, sino que vamos a resolver el problema al completo. De este modo, nuestra funci´on de p´erdida, ser´a la suma de los errores cuadr´aticos medios cometidos al predecir el valor de la soluci´on de la ecuaci´on en el interior, las condiciones de contorno en la frontera y/o las condiciones en el tiempo inicial. Utilizaremos m´as de una capa oculta para la resoluci´on de los problemas planteados en la presente memoria, del mismo modo que en el Cap´ıtulo 4para aproximar funciones. Teniendo esto en cuenta, estamos en posici´on de resolver diversos problemas utilizando la librer´ıa Tensorflow, que es muy conocida para la creaci´on de modelos de aprendizaje autom´atico, cuya documentaci´on puede consultarse en [20]. Adem´as, hay que mencionar que TensorFlow est´a disponible para m´ultiples lenguajes, como Matlab, C++, Java, Haskell, R... No obstante, la librer´ıa m´as completa, sencilla de usar y revisada es la de Python, que es con la que vamos a trabajar. Concluimos esta secci´on aportando un enlace con las instrucciones de instalaci´on de Tensorflow en nuestra distribuci´on de Python seg´un el sistema operativo: https://www.tensorflow.org/install/pip?hl=es. As´ı, podremos importar la librer´ıa en los c´odigos realizados en las secciones posteriores. 65 5.2. El problema de Laplace En esta secci´on resolvemos el problema de Laplace en dimensi´on 2 sujeto a diversas condiciones utilizando la librer´ıa TensorFlow de Python. En particular, resolvemos los siguientes ejemplos a cada uno de los cuales dedicamos una subsecci´on: Ejemplo 1: problema de Laplace con segundo miembro constante y condiciones de contorno de tipo Dirichlet homog´eneas en el cubo unidad. Ejemplo 2: problema de Laplace con segundo miembro constante y condiciones de contorno de tipo Dirichlet homog´eneas en la bola unidad (en este caso la soluci´on anal´ıtica es conocida). Ejemplo 3: problema de Laplace con segundo miembro no constante y condiciones de contorno de tipo Dirichlet homog´eneas en el cubo unidad (en este caso la soluci´on anal´ıtica es conocida). Ejemplo 4: problema de Laplace con segundo miembro no constante y condiciones de contorno de tipo Neumann en la bola unidad. Ejemplo 5: problema de Laplace no lineal con segundo miembro constante y condiciones de contorno de tipo Dirichlet homog´eneas en la bola unidad. 5.2.1. Ejemplo 1: problema con segundo miembro constante en el cubo unidad Comenzamos resolviendo el siguiente problema: −∆u= 1 en Ω u= 0 sobre ∂Ω,(5.1) donde Ω = (0,1)×(0,1). Para resolver este problema con redes neuronales, tenemos que escribirlo como un problema de optimizaci´on siguiendo las ideas propuestas en 3. Concretamente, vamos a considerar ˆ Ω, un subconjunto de Ω de Npuntos, y ˆ ∂Ω, un subconjunto de ∂Ω de Mpuntos, que usualmente generaremos de forma aleatoria. En estas condiciones, planteamos el problema de optimizaci´on m´ın p   1 NX x,y∈ˆ Ω (1 + ∆up(x, y))2+k MX x,y∈∂ˆ Ω up(x, y)2   ,(5.2) con k > 0 es una constante de penalizaci´on positiva. Observaci´on 5.1. A diferencia del m´etodo propuesto en el Cap´ıtulo 3, consideramos las restricciones del problema como parte de la funci´on de p´erdida a trav´es del t´ermino k MX x,y∈∂ˆ Ω up(x, y)2, 66 que cuando kes grande fuerza a upa ser nulo en la frontera, de forma similar al bloqueo en el m´etodo de los elementos finitos. Esto resulta m´as eficiente desde el punto de vista computacional y nos evita escoger una funci´on que se anule la frontera como ocurr´ıa en 3. Seguiremos esta t´ecnica en el resto de problemas resueltos en la memoria. Observaci´on 5.2. Tomando |Ω|como la medida de Lebesgue de Ω en R2y|∂Ω| como la medida longitudinal de ∂Ω, observamos que las sumas que aparecen en (5.2) son las aproximaciones de Montecarlo correspondientes de las integrales 1 NX x,y∈∂ˆ Ω up(x, y)2∼1 |Ω|ZΩ (1 + ∆u)2dx, k MX x,y∈∂ˆ Ω up(x, y)2∼k |∂Ω|Z∂Ω u2ds(x). De esta forma, (5.2) se puede considerar una discretizaci´on del problema m´ın u∈H2(Ω) 1 |Ω|ZΩ (1 + ∆u)2dx +k |∂Ω|Z∂Ω u2ds(x).(5.3) Recordemos que el problema (5.1) posee efectivamente una ´unica soluci´on en H2(Ω)∩H1 0(Ω), la cual adem´as es anal´ıtica en Ω, positiva en el interior y es invariante respecto a las rotaciones y simetr´ıas que dejan invariante el cuadrado. Procedemos a la implementaci´on de (5.2) con TensorFlow. En primer lugar, importamos las librer´ıas que vamos a utilizar. Importamos numpy y matplotlib al igual que en la secci´on 3.2, adem´as de TensorFlow: import tensorflow as tf import numpy as np import matplotlib . pyplot as plt A continuaci´on, creamos una clase python cuyos objetos son redes neuronales con la arquitectura escogida para el problema a resolver (esto se incluye en el m´etodo init que es el constructor de los objetos la clase). Asimismo, debemos incluir el m´etodo call que, dada una red de la clase, se encarga de calcular su salida. En particular, para este problema hemos escogido una arquitectura con 5 capas ocultas con 50 neuronas cada una. Adem´as, utilizaremos como funci´on de activaci´on la tangente hiperb´olica, ya que en este problema aparecen derivadas de segundo orden y la ReLU no es ni siquiera derivable. class PINN ( tf. keras . Model ): def __init__ ( self ): super(PINN , self ). __init__ () self . dense1 = tf . keras . layers . Dense (50 , activation = ’tanh ’) self . dense2 = tf . keras . layers . Dense (50 , activation = ’tanh ’) self . dense3 = tf . keras . layers . Dense (50 , activation = ’tanh ’) 67 ERROR L2: 0.0009906976 ERROR H1: 0.016927944 Tal y como cabe esperar, el error en L2(Ω) es inferior al obtenido en H1(Ω). 5.2.2. Ejemplo 2: problema con segundo miembro constante en la bola unidad El m´etodo estudiado en el Cap´ıtulo 3, estaba dise˜nado para ser aplicado en cubos. Sin embargo, el m´etodo que hemos usado para la resoluci´on de (5.1), donde la discretizaci´on del abierto estaba definida a trav´es de una nube de puntos escogida de forma aleatoria en el conjunto y otra en la frontera, es una t´ecnica que podemos extender a cualquier abierto regular Ω. Aunque Ω no sea un cubo, el problema de minimizaci´on (5.2) se puede interpretar como una discretizaci´on de (5.3) y, por tanto, como una aproximaci´on de: −∆u= 1 en Ω u= 0 sobre ∂Ω.(5.4) En esta secci´on consideramos como ejemplo Ω = B((0,0),1). En este caso, la soluci´on anal´ıtica es conocida y viene dada por u(x, y) = 1−x2−y2 4, de forma que podemos representar los resultados obtenidos por nuestra aproximaci´on sin necesidad de usar FreeFem u otro software similar. Vamos a representar gr´aficamente dicha soluci´on anal´ıtica: Figura 5.7: Soluci´on anal´ıtica de (5.4). 74 Vamos a resolver ahora el problema utilizando TensorFlow. Para ello, generamos aleatoriamente puntos en el cubo unidad excluyendo del conjunto de puntos de entrada todos aquellos que no est´en en la bola hasta completar una cantidad prefijada N1, que en este caso es 5000. Los puntos de la frontera se generan sobre la circunferencia utilizando la parametrizaci´on (x, y) = (cos(θ),sin(θ)), θ ∈[0,2π) y tomando θdefinido sobre una malla de [0,2π) con una cantidad de puntos N2 prefijada, que en el test realizado es 500. Por tanto, sustituimos la generaci´on de datos de entrenamiento del problema en el cubo unidad por las siguientes l´ıneas de c´odigo: # Generar datos de entrenamiento en el circulo def generate_circle_points (N, radius =0.5 , center =(0.5 , 0.5) ): X = [] while len (X) < N: x, y = 2*np. random .rand (2) -1 if (x - center [0]) ** 2 + (y - center [1]) ** 2 <= radius ** 2: X. append ([x, y]) return np . array (X) # Generar datos de entrenamiento en la frontera del circulo def generate_boundary_points (N, radius =0.5 , center =(0.5 , 0.5) ): theta = np. linspace (0, 2 * np.pi , N) X_boundary = np. array ([ center [0] + radius * np.cos( theta ), center [1] + radius * np.sin( theta )]).T return X_boundary , np.zeros ((N, 1)) # Generar datos de entrenamiento N1 =5000 N2 = 500 X = generate_circle_points (N1 ,1 ,(0 ,0)) X_boundary , u_boundary = generate_boundary_points (N2 ,1,(0,0)) La conversi´on a elementos de tipo TensorFlow y el proceso de creaci´on y entrenamiento del modelo es id´entico al problema del cubo. Una vez entrenado el modelo, representamos la soluci´on obtenida del siguiente modo: # Generar una cuadricula de puntos en el dominio [0, 1] x [0, 1] 75 num_points = 100 x = np. linspace (-1, 1, num_points ) y = np. linspace (-1, 1, num_points ) X_grid , Y_grid = np . meshgrid (x, y) X_flat = X_grid . flatten () Y_flat = Y_grid . flatten () # Excluimos los puntos fuera del circulo mask = ( X_flat ) ** 2 + ( Y_flat ) ** 2 <= 1 X_test = np. vstack ([ X_flat [ mask ], Y_flat [ mask ]]) .T # Convertir los puntos de prueba a tensores de TensorFlow X_test = tf. convert_to_tensor ( X_test , dtype =tf. float32 ) # Usar el modelo para predecir la solucion en los puntos de prueba u_pred = model ( X_test ) # Aplanar las predicciones u_pred_flat = u_pred . numpy (). flatten () # Dar forma a las predicciones para que coincidan con la cuadricula u_pred_grid = np.full (( num_points , num_points ), np.nan ) u_pred_grid [ mask . reshape ( num_points , num_points )] = u_pred_flat # Graficar la solucion predicha plt . figure ( figsize =(8 , 6)) plt . contourf ( X_grid , Y_grid , u_pred_grid , levels =50 , cmap =’ jet ’) plt . colorbar ( label = ’u(x ,␣y) ’) plt . xlabel (’x’) plt . ylabel (’y’) plt . title (’Predicted ␣ Solution ␣of␣ the ␣PDE ␣( Laplacian ␣=␣1)␣in ␣a␣ Circular ␣ Domain ’) Mostramos la aproximaci´on obtenida tras iterar durante 1500 ´epocas y la evoluci´on del funcional loss: 76 Figura 5.8: Aproximaci´on de la soluci´on de (5.4) tras 1500 iteraciones o ´epocas. Figura 5.9: Evoluci´on de loss tras 1500 iteraciones o ´epocas. En la siguiente figura representamos la diferencia entre la soluci´on proporcionada por el algoritmo y la soluci´on anal´ıtica, mostrando que el m´etodo proporciona una buena aproximaci´on: 77 Figura 5.10: Error cometido al aproximar (5.4) con ANN tras 1500 iteraciones o ´epocas. M´as concretamente, calculando los errores cometidos en L∞(Ω), L2(Ω), H1(Ω) yH2(Ω), se tiene: ERROR Linf: 0.00099805 ERROR L2: 0.00034512096 ERROR H1: 0.0021254956 ERROR H2: 0.014625692 Tal y como cabe esperar, el error cometido en L2(Ω) es inferior al correspondiente a H1(Ω) y este a su vez es menor que el error respecto a la norma de H2(Ω). 5.2.3. Ejemplo 3: problema con segundo miembro no constante en el cubo unidad Mostramos en esta subsecci´on c´omo el m´etodo usado anteriormente puede extenderse sin dificultades al caso en que el segundo miembro no es constante. A t´ıtulo de ejemplo, resolvemos a continuaci´on el siguiente problema: −∆u=f(x, y) en Ω, u= 0 sobre ∂Ω,(5.5) con Ω = (0,1) ×(0,1), f(x, y) = x+y−x2−y2. Al igual que en el caso anterior, se trata de un problema en que la soluci´on anal´ıtica es conocida y viene dada por: u(x, y) = x(1 −x)y(1 −y) 2,∀x, y ∈Ω. 78 Representamos dicha soluci´on anal´ıtica en la siguiente figura: Figura 5.11: Soluci´on anal´ıtica de (5.5). An´alogamente a (5.1), vamos a resolver el siguiente problema de optimizaci´on m´ın p   1 NX x,y∈ˆ Ω (f(x, y)+∆up(x, y))2+k MX x,y∈∂ˆ Ω up(x, y)2   ,(5.6) con k > 0 una constante de penalizaci´on positiva. El ´unico cambio que hay que realizar en el c´odigo con respecto a la resoluci´on de (5.1) es cambiar la funci´on de p´erdida para la ecuaci´on diferencial del siguiente modo: def pde_loss ( model , X): with tf. GradientTape ( persistent = True ) as tape : tape . watch (X) u = model (X) u_x = tape . gradient (u, X)[:, 0] u_y = tape . gradient (u, X)[:, 1] u_xx = tape. gradient (u_x , X)[: , 0] u_yy = tape. gradient (u_y , X)[: , 1] del tape return tf . reduce_mean ( tf. square (- u_xx - u_yy + np. power (X[: ,0] ,2) + np. power (X[: ,1] ,2) -X [: ,0] -X[: ,1]) ) Las representaciones de la aproximaci´on de la soluci´on obtenida y de la evoluci´on del funcional loss tras 7500 ´epocas son: 79 Figura 5.12: Aproximaci´on de la soluci´on de (5.5) con ANN tras 7500 iteraciones o ´epocas. Figura 5.13: Evoluci´on de loss tras 7500 iteraciones o ´epocas. Representamos el error cometido respecto a la soluci´on anal´ıtica en la siguiente figura para comprobar la calidad de la aproximaci´on 80 Figura 5.14: Error cometido al aproximar (5.5) con ANN tras 7500 iteraciones o ´epocas. M´as concretamente, calculando los errores cometidos en L∞(Ω), L2(Ω), H1(Ω) yH2(Ω), se tiene: ERROR Linf: 0.003011947963386774 ERROR L2: 0.0006554645 ERROR H1: 0.0052606706 ERROR H2: 0.03906873 5.2.4. Ejemplo 4: problema con condiciones de Neumann En los ejemplos anteriores hemos estado considerando condiciones de Dirichlet; en esta subsecci´on, vemos c´omo el m´etodo puede adaptarse para resolver problemas con otras condiciones de contorno. Como ejemplo, resolvemos el siguiente problema con condiciones de Neumann: −∆u+u=f(x, y) en Ω ∂νu= 0 sobre ∂Ω,(5.7) donde Ω = B((0,0),1) y f(x, y) = x2+y2. Para aproximar la soluci´on del problema mediante redes neuronales, resolvemos el siguiente problema de optimizaci´on considerando ˆ Ω, un subconjunto de Ω de Npuntos, y ˆ ∂Ω, un subconjunto de ∂Ω de M puntos: m´ın p   1 NX x,y∈ˆ Ω (f(x, y)+∆up(x, y)−up(x, y))2+k MX x,y∈∂ˆ Ω (∂νup(x, y))2   ,(5.8) 81 con k > 0 una constante de penalizaci´on. En este caso, sabemos que la soluci´on es radial pero no la conocemos explicitamente, por lo que para comprobar la fiabilidad del m´etodo, compararemos la aproximaci´on obtenida mediante ANN con la dada por una malla fina y elementos de grado 1 en FreeFem. Adem´as de modificar las funciones de p´erdida en el c´odigo para resolver (5.4), al generar los datos de entrenamiento creamos un nuevo conjunto de puntos X closure que consiste en unir los puntos generados en el c´ırculo y en la frontera. Esto es necesario porque hay que calcular las derivadas en todo el c´ırculo pero solo evaluarlas en la frontera para implementar la condici´on de Neumann y en TensorFlow solo podemos evaluar las derivadas en el conjunto en que han sido calculadas. Por tanto, la generaci´on de datos de entrenamiento queda como sigue: # Generar datos de entrenamiento N1 =5000 N2 = 500 X = generate_circle_points (N1 ,1 ,(0 ,0)) X_boundary , u_boundary = generate_boundary_points (N2 ,1,(0,0)) X_closure = np. append (X, X_boundary , axis =0) # Convertir a tensores de TensorFlow X = tf. convert_to_tensor (X, dtype =tf. float32 ) X_boundary = tf. convert_to_tensor ( X_boundary , dtype =tf. float32 ) u_boundary = tf. convert_to_tensor ( u_boundary , dtype =tf. float32 ) X_closure = tf . convert_to_tensor ( X_closure , dtype = tf . float32 ) As´ı, como el m´etodo append de numpy a˜nade los puntos de X boundary a continuaci´on de los de Xy la variable N2 contiene el n´umero de puntos de la frontera, podemos evaluar las derivadas parciales de los elementos de X closure en los puntos de X boundary de forma simple. Adem´as, como Ω = B((0,0),1), se tiene ∂νup=x∂xup+y∂yu. De este modo, las funciones de p´erdida para este problema son: # Definir la funcion de perdida para el PDE def pde_loss ( model , X): with tf. GradientTape ( persistent = True ) as tape : tape . watch (X) u = model (X) u_x = tape . gradient (u, X)[:, 0] u_y = tape . gradient (u, X)[:, 1] u_xx = tape. gradient (u_x , X)[: , 0] u_yy = tape. gradient (u_y , X)[: , 1] 82 del tape return tf . reduce_mean ( tf. square (- u_xx - u_yy + u - np.power (X [: ,0] ,2) - np. power (X[: ,1] ,2))) # Definir la funcion de perdida para la condicion de frontera def boundary_loss (model , X_boundary , u_boundary , X_closure , N2): with tf. GradientTape ( persistent = True ) as tape : tape . watch ( X_closure ) u = model ( X_closure ) u_x = tape . gradient (u, X_closure )[-N2 , 0] u_y = tape . gradient (u, X_closure )[-N2 , 1] del tape return tf. reduce_mean (tf. square ( u_x * X_boundary [: ,0] + u_y * X_boundary [: ,1] - u_boundary )) En la figura 5.15 representamos la aproximaci´on de (5.7) obtenida tras 1700 iteraciones, en la figura 5.16 la evoluci´on del funcional loss, en la figura 5.17 la aproximaci´on obtenida utilizando Freefem y, por ´ultimo, en la figura 5.18, la diferencia entre las aproximaciones obtenidas por ambos m´etodos. Figura 5.15: Aproximaci´on de la soluci´on de (5.7) con ANN tras 1700 iteraciones o ´epocas. 83 Figura 5.21: Evoluci´on de loss tras 1500 iteraciones o ´epocas. Finalmente, mostramos el error cometido con respecto a la aproximaci´on dada por FreeFem: Figura 5.22: Error cometido al aproximar (5.9) con ANN tras 1500 iteraciones. Se comete un error del orden de 10−3con respecto a la aproximaci´on de Freefem y, efectivamente, no hemos necesitado ning´un m´etodo iterativo adicional. M´as concretamente, calculando los errores cometidos en L∞(Ω), L2(Ω) y H1(Ω), se tiene: ERROR Linf: 0.0021032275781631415 ERROR L2: 0.000941068 ERROR H1: 0.010875117 90 5.3. El problema de Laplace como problema de c´alculo de variaciones Volvemos a resolver (5.1) esta vez como problema de c´alculo de variaciones recordando que su soluci´on coincide con la del problema de optimizaci´on m´ın u∈H1(Ω) (1 |Ω|ZΩ1 2|∇u|2−u2 dx dy +k1 |∂Ω|Z∂Ω u2dx), que podemos resolver directamente utilizando redes neuronales siendo kuna constante de penalizaci´on positiva .Concretamente, para dos subconjuntos ˆ Ω⊂Ω de N puntos y ∂ˆ Ω⊂∂Ω de Mpuntos, nuestro algoritmo resuelve m´ın p   1 NX x,y∈ˆ Ω1 2|∇up(x, y)|2−up(x, y)2 +k1 MX x,y∈∂ˆ Ω u2 p(x, y)   . De esta forma, para la resoluci´on con TensorFlow, el c´odigo de las funciones de p´erdida es el siguiente: # Definir la funcion de perdida para la PDE (Laplaciano = 1) def valor_objetivo (model , X): with tf. GradientTape ( persistent = True ) as tape : tape . watch (X) u = model (X) u_x = tape . gradient (u, X)[:, 0] u_y = tape . gradient (u, X)[:, 1] del tape mod_u =tf. norm (tf. stack ([ u_x , u_y ], axis = -1) , axis = -1) return tf . reduce_mean ( tf. square ( mod_u ) *0.5 - u) # Definir la funcion de perdida para la condicion de frontera def penalizacion ( model , X_boundary ,pen ): u_pred = model ( X_boundary ) return pen *tf. reduce_mean ( tf. square ( u_pred )) Mostramos la aproximaci´on obtenida tras 7500 ´epocas en la figura 5.23 y la evoluci´on del funcional loss en la figura 5.24: 91 Figura 5.23: Aproximaci´on resolviendo (5.1) como problema de c´alculo de variaciones con ANN tras 7500 iteraciones. Figura 5.24: Evoluci´on de loss tras 7500 iteraciones o ´epocas. Finalmente, mostramos el error cometido con respecto a la aproximaci´on que obtuvimos utilizando FreeFem en la subsecci´on (5.2.1) (ver figura 5.5): 92 Figura 5.25: Error cometido en la aproximaci´on de (5.2.1) resolvi´endolo como problema de c´alculo de variaciones respecto a la aproximaci´on en FreeFem representada en 5.5. M´as concretamente, calculando los errores cometidos en L∞(Ω), L2(Ω) y H1(Ω), se tiene: ERROR Linf: 0.010094239192390442 ERROR L2: 0.005254852 ERROR H1: 0.02667785 Dado que en este problema no intervienen derivadas segundas y en el cap´ıtulo anterior nos centramos en la ReLU para dar estimaciones de error, se ha decidido repetir el experimento tomando como funci´on de activaci´on la ReLU en lugar de la tangente hiperb´olica. Mostramos la aproximaci´on obtenida tras 7500 ´epocas en la figura 5.26 y la evoluci´on del funcional loss en la figura 5.27: 93 Figura 5.26: Aproximaci´on resolviendo (5.1) como problema de c´alculo de variaciones con ANN tras 7500 iteraciones. Figura 5.27: Evoluci´on de loss tras 7500 iteraciones o ´epocas. Finalmente, mostramos el error cometido con respecto a la aproximaci´on que obtuvimos utilizando FreeFem en la subsecci´on (5.2.1) (ver figura 5.5): 94 Figura 5.28: Error cometido en la aproximaci´on de (5.2.1) resolvi´endolo como problema de c´alculo de variaciones respecto a la aproximaci´on en FreeFem representada en 5.5. M´as concretamente, calculando los errores cometidos en L∞(Ω), L2(Ω) y H1(Ω), se tiene: ERROR Linf: 0.019560780930233 ERROR L2: 0.0045148563 ERROR H1: 0.077738 Observamos que los errores obtenidos son similares a los obtenidos en el caso anterior. 5.4. Problema de control Teniendo en cuenta la capacidad de las redes neuronales para tratar problemas de optimizaci´on, mostramos en esta secci´on c´omo podemos aplicarlas a la resoluci´on de un problema de control simple. Consideramos el siguiente problema modelo: Dado Ω ⊂RNabierto, acotado y α > 0, queremos resolver m´ın v∈L2(Ω) {J(v)},(5.11) donde el funcional coste es J(v) := 1 2ZΩ |uv−1|2dx +α 2ZΩ |v|2dx 95 y para cada v∈L2(Ω), la funci´on de estado uves la soluci´on de la ecuaci´on de estado (−∆uv=ven Ω uv= 0 sobre ∂Ω. Desde el punto de vista te´orico, se tiene: Teorema 5.4. En las condiciones descritas, el problema (5.11) tiene una ´unica soluci´on v∈L2(Ω). Tomando ula funci´on de estado correspodiente y definiendo el estado adjunto asociado ϕsoluci´on de (−∆ϕ=u−1en Ω ϕ= 0 sobre ∂Ω,(5.12) se tiene ϕ=−αv en Ω. Demostraci´on. Teniendo en cuenta que por el teorema de Lax-Milgram la aplicaci´on v∈L2(Ω) 7→ uv∈H1 0(Ω) es lineal y continua, se tiene que el funcional es estrictamente convexo, continuo y coercitivo sobre el espacio de Hilbert H1 0(Ω). Esto implica la existencia y unicidad de soluci´on de este problema. Adem´as, gracias a la convexidad, sabemos que una funci´on ves soluci´on si y solo si la derivada direccional del funcional coste en vy en cualquier direcci´on w∈L2(Ω) es nula. La derivada del funcional coste viene dada por ⟨DJ(v), w⟩=d dε 1 2ZΩ |uv+εw −1|2dx +α 2ZΩ |v+εw|2dx|ε=0 , donde gracias a la linealidad de la ecuaci´on de estado se verifica uv+εw =u+εz con zsoluci´on de (−∆z=wen Ω z= 0 sobre ∂Ω.(5.13) La condici´on de optimalidad queda por tanto ZΩ(u−1)z+αvwdx = 0,∀w∈L2(Ω). Por otra parte, usando (5.12) y (5.13) se tiene ZΩ (u−1)z dx =ZΩ ∇ϕ· ∇z dx =ZΩ wϕ dx, luego la condici´on queda ZΩ (ϕ+αv)w dx = 0,∀w∈L2(Ω), lo que es equivalente a ϕ+αv = 0 en Ω. 96 Observaci´on 5.5. Por el Teorema 5.4 tenemos que el control ´optimo vy la funci´on de estado asociada uest´an caracterizados como las soluciones del problema de contorno          −∆u=ven Ω −∆v=1 α(1 −u) en Ω u=v= 0 sobre ∂Ω. (5.14) Desde el punto de vista de las redes neuronales, podemos usar la aproximaci´on de (5.11) dada por m´ın u∈H2(Ω),v∈L2(Ω) 1 2|Ω|ZΩ |u−1|2dx +α 2|Ω|ZΩ |v|2dx +k1 |Ω|ZΩ |v+ ∆u|2dx +k2 |∂Ω|Z∂Ω |u|2ds(x), con k1, k2constantes de penalizaci´on positivas a escoger. Considerando dos subconjuntos ˆ Ω⊂Ω de Npuntos y ∂ˆ Ω⊂∂Ω de Mpuntos, nuestro algoritmo resuelve m´ın p   1 2NX x,y∈ˆ Ω |up(x, y)−1|2+α 2NX x,y∈ˆ Ω |vp(x, y)|2+k1 NX x,y∈ˆ Ω |vp(x, y) +∆up(x, y)|2+k2 MX x,y∈∂ˆ Ω |up(x, y)|2   . Para la implementaci´on en TensorFlow, hemos tomado Ω = B((0,0),1) y hemos considerado 5 capas ocultas con 50 neuronas en la primera y ´ultima y 60 en el resto, siguiendo una arquitectura del tipo de la heur´ıstica de la Figura 2.7. Adem´as, hemos elegido k1= 10, k2= 40 y α= 0′05. En el c´odigo utilizado para resolver (5.4), la parte correspondiente a la funci´on de p´erdida y el entrenamiento la modificamos del siguiente modo: def pde_loss ( model , X): with tf. GradientTape ( persistent = True ) as tape : tape . watch (X) u_v = model (X) u = u_v [: , 0] v = u_v [: , 1] u_x = tape . gradient (u, X)[:, 0] u_y = tape . gradient (u, X)[:, 1] u_xx = tape. gradient (u_x , X)[: , 0] u_yy = tape. gradient (u_y , X)[: , 1] del tape return tf. reduce_mean (tf.square ( u_xx + u_yy + v)) 97 def quadrature_loss ( model ): u_v = model (X) u = u_v [: , 0] v = u_v [: , 1] integral_u = tf. reduce_mean (tf.square (u - 1)) integral_v = tf. reduce_mean (tf. square (v)) return integral_u , integral_v for epoch in range ( epochs ): with tf. GradientTape ( persistent = True ) as tape : loss_pde = pde_loss (model , X) loss_boundary = boundary_loss ( model , X_boundary , u_boundary) integral_u , integral_v = quadrature_loss ( model ) loss = 5* loss_pde + 20* loss_boundary + integral_u + 0.05*integral_v #loss = integral_u gradients = tape . gradient (loss , model . trainable_variables) del tape optimizer . apply_gradients (zip( gradients , model . trainable_variables)) Tras ejecutar el c´odigo, las aproximaciones obtenidas para uyvmediante ANN tras 4500 iteraciones y la evoluci´on del funcional loss son los siguientes: Figura 5.29: Aproximaciones de la funci´on de estado uy el control vde (5.11) con ANN tras 4500 ´epocas. 98 Figura 5.30: Evoluci´on de loss durante 4500 ´epocas. Adem´as, resolviendo el problema de contorno que aparece en la Observaci´on 5.5 utilizando FreeFem, obtenemos Figura 5.31: Aproximaciones de la funci´on de estado uy el control vde (5.11) obtenidas resolviendo con FreeFem el problema de contorno de la Observaci´on 5.5. Finalmente, representamos el error cometido Figura 5.32: Errores cometidos al aproximar la funci´on de estado uy el control vde (5.11) con ANN tras 4500 ´epocas respecto a la resoluci´on del problema de contorno. 99 Observamos que, contrariamente a los m´etodos cl´asicos para la resoluci´on de problemas con condiciones iniciales, los valores en un tiempo dado no se obtienen a partir de los valores en tiempos anteriores. El algoritmo trata con todos los tiempos de forma simult´anea. Antes de proceder a la implementaci´on con TensorFlow, supongamos λ= k= 1 (problema adimensionalizado) y procedemos a resolver el problema mediante el siguiente esquema pseudoimpl´ıcito con Matlab para conocer el comportamiento esperado para el sistema generando una animaci´on: un+1 −un h−∂2 xxun+1 =λun1−un k, con h > 0 el paso de tiempo. Mostramos el c´odigo Matlab correspondiente: u0 = @(x) 1 + cos(x); u0c = @(x) [1.1 + cos(x); -sin(x)]; h = 0.01; ih = 1 / h; xd = linspace(0, pi , 20); % Crear el objeto VideoWriter con el perfil ’MPEG -4 ’ vid = VideoWriter (’ logistica .mp4 ’, ’MPEG -4 ’); open ( vid); for j = 0:500 plot(xd , u0(xd)); axis([0 pi -0.2 2.2]) ; frame = getframe(gcf); writeVideo (vid , frame ); % Guardar frames especificos como imagenes JPG if ismember (j , [0 , 50 , 100 , 200 , 300 , 400 , 500]) filename = sprintf (’ frame_ %d. jpg ’, j); saveas(gcf, filename); end fun = @(x, y) [y (2) ; ih * y(1) - (1 + ih) * u0(x) + u0( x)ˆ2]; condcon = @(ya , yb) [ya (2); yb (2) ]; solinit = bvpinit (xd , u0c); sol = bvp4c (fun , condcon , solinit ); u0 = @(x) deval (sol , x, 1) ; u0c = @(x) deval (sol , x); end % Cerrar el archivo de video close(vid ); 106 Representamos algunos instantes de la animaci´on generada: Figura 5.35: Aproximaciones de (5.15) con Matlab en los instantes 0, 1 y 1′5. Figura 5.36: Aproximaciones de (5.15) con Matlab en los instantes 2, 3, 4 y 5. Observamos que, conforme avanza el tiempo, la soluci´on se estabiliza en torno a la soluci´on constante u= 1. Para la implementaci´on con TensorFlow del problema, hemos dividido los intervalos (0,5) y (0, π) en 100 partes iguales, de modo que M1=M2= 100. Adem´as, hemos tomado ˆ Ω (0,5) ×(0, π) de modo que para cada uno de los 100 tiempos del subconjunto ∂ˆ Ω1de (0,5) se consideren todos los 100 puntos del subconjunto ∂ˆ Ω2 de (0, π), i.e., para cada tiempo tenemos todos los espacios. Realmente pod´ıamos haber tomado discretizaciones aleatorias con dos dimensiones igual que hemos hecho en otros ejemplos, pero en los experimentos realizados obtuvimos peores resultados y al no haber escogido muchos puntos no supone un gran coste computacional. La generaci´on de puntos de entrenamiento queda como sigue: N=100 T=np. linspace (0 ,5 ,N) T=T. tolist () X = np. linspace (0, np.pi , N). flatten () X_T = [] X_T_inic=[] for tin T: for xin X: X_T . append ([t ,x]) if t==0: 107 X_T_inic . append ([t ,x]) X_T = np . array ( X_T ) X_T_inic = np. array ( X_T_inic ) X_T_boundary = [] for tin T: X_T_boundary . append ([t ,0]) X_T_boundary . append ([t,np .pi ]) X_T_boundary = np.array ( X_T_boundary ) N2 =np. shape ( X_T_boundary ) [0] u_boundary = np.zeros ((N2 , 1)) X_closure = np. append (X_T , X_T_boundary , axis =0) # Convertir a tensores de TensorFlow X = tf . convert_to_tensor (X , dtype =tf. float32 ) X = tf. reshape (X, (-1, 1)) # Reformatear a (1000 , 1) X_T = tf. convert_to_tensor (X_T , dtype =tf.float32 ) X_T_inic = tf. convert_to_tensor ( X_T_inic , dtype = tf. float32 ) u_inic = tf.cos(X) X_T_boundary = tf. convert_to_tensor ( X_T_boundary , dtype = tf. float32 ) u_boundary = tf. convert_to_tensor ( u_boundary , dtype =tf. float32 ) X_closure = tf . convert_to_tensor ( X_closure , dtype = tf . float32 ) Adem´as, observamos que, al igual que para el problema (5.7), hemos almacenado los puntos de la clausura en un conjunto para implementar el funcional de coste correspondiente a las condiciones de contorno de tipo Neumann del siguiente modo: # Definir la funcion de perdida para la condicion de frontera def boundary_loss ( model , X_closure , u_boundary ): with tf. GradientTape ( persistent = True ) as tape : tape . watch ( X_closure ) u = model ( X_closure ) u_x = tape . gradient (u, X_closure )[: , 1] del tape return tf. reduce_mean (tf. square ( u_x [- N2] - u_boundary )) El c´odigo correspondiente a los residuos correspondientes a la ecuaci´on diferencial y la condicion inicial es: lam =1 k=1 kinv =1/k 108 # Definir la funcion de perdida para la PDE (Laplaciano = 1) def pde_loss ( model , X_T ): with tf. GradientTape ( persistent = True ) as tape : tape . watch ( X_T ) u = model ( X_T ) u_t = tape . gradient (u, X_T )[:, 0] u_x = tape . gradient (u, X_T )[:, 1] u_xx = tape. gradient (u_x , X_T )[: , 0] del tape return tf . reduce_mean ( tf. square ( u_t -u_xx - lam *u*(1 - u* kinv ))) # Definir la funcion de perdida para la condicion inicial def inic_loss ( model , X_T_inic ): u_pred = model ( X_T_inic ) return tf . reduce_mean ( tf. square ( u_pred - u_inic - 1.1) ) Una vez entrenado el modelo, representamos mediante una animaci´on el comportamiento del sistema. Para esto, importamos la librer´ıa animation de matplotlib incluyendo la siguiente l´ınea de c´odigo en la cabecera del programa import matplotlib.animation as animation y ejecutamos # Generar una cuadricula de puntos en el dominio [0, 1] x [0, 1] num_points = 1000 X = np. linspace (0, np.pi , num_points ) fig = plt . figure () # definimos una funcion auxiliar para actualizar la posicion del cono en cada instante def actualizar(i): X_test = [] for xin X: X_test . append ([T[i],x ]) # Convertir los puntos de prueba a tensores de TensorFlow X_test = tf. convert_to_tensor ( X_test , dtype =tf. float32 ) # Usar el modelo para predecir la solucion en los puntos de prueba u_pred = model ( X_test ) # Aplanar las predicciones u_pred_flat = u_pred . numpy (). flatten () 109 fig . clear () # dibujamos la circunferencia completa plt . plot ( X_test [: ,1] , u_pred_flat ) plt . xlim (0,np.pi) plt . ylim ( -0.5 ,2.5) plt . title ( ’Predicted ␣ Solution ␣ of␣the ␣ PDE ’) writervideo = animation . FFMpegWriter (fps =5) ani = animation . FuncAnimation (fig , actualizar , range(len(T)), interval =1, repeat = False ) ani .save (’u.mp4 ’, writer = writervideo ) print(" Animacion ␣ generada ") Mostramos la evoluci´on del funcional coste y algunos instantes de una animaci´on generada tras entrenar el modelo 500 ´epocas: Figura 5.37: Evoluci´on de loss durante 500 ´epocas. Figura 5.38: Aproximaciones de (5.15) con ANN para 500 ´epocas en los instantes 0, 1 y 2. 110 Figura 5.39: Aproximaciones de (5.15) con ANN para 500 ´epocas en los instantes 3, 4 y 5. Observamos que conforme avanza el tiempo, la aproximaci´on de la soluci´on obtenida para utiende a estabilizarse entorno al valor constante u= 1, tal y como cabr´ıa esperar. Para finalizar, cabe destacar que aumentar el n´umero de puntos empeora el resultado en este problema, ya que se produce un sobreajuste u overfitting hacia los puntos del intervalo (0, π/2), de modo que extiende el comportamiento esperado para el sistema en estos puntos al intervalo (π/2, π). Mostramos un ejemplo de esto con M2= 300: Figura 5.40: Aproximaciones de (5.15) con ANN para 1000 ´epocas en los instantes 0, 1 y 2. Figura 5.41: Aproximaciones de (5.15) con ANN para 1000 ´epocas en los instantes 3, 4 y 5. 111 5.6. Problema inverso: modelo de Lotka-Volterra Consideramos el sistema depredador-presa (u′=−αu +βuv, en [a, b] v′=γv −δuv, en [a, b] con T > 0, uyvla cantidad de predadores y presas respectivamente por unidad de tiempo. Nos planteamos el problema inverso correspondiente a encontrar las constantes α, β, γ, δ ≥0 a partir de unas mediciones ui,vide los valores u, v en tiempos t1, ..., tn. A fin de adaptar la resoluci´on de este problema a las redes neuronales podemos pensar en minimizar en u, v, α, β, γ, δ la cantidad 1 n n X i=1 |u(ti)−ui|2+|v(ti)−vi|2+AZb a|u′+αu +βv|2+|v′−γv +δuv|2dt, siendo A > 0 un valor de penalizaci´on alto y [a, b] un intervalo que contenga a todos los puntos ti, por ejemplo se puede tomar a= m´ın i≤i≤nti, b = m´ax 1≤i≤nti. El problema es que los par´ametros α, β, γ, δ no vienen dados a partir de funciones definidas por la red neuronal. Una forma (costosa desde el punto de vista computacional) de arreglar esta dificultad consiste en transformar el sistema en uno de 6 ecuaciones y 6 inc´ognitas dado por        u′=−αu +βuv v′=γv −δuv α′=β′=γ′=δ′= 0. El problema entonces se aproxima por m´ın u,v,α,β,γ,δ (1 n n X i=1 |u(ti)−ui|2+|v(ti)−vi|2 +AZb a|u′+αu +βv|2+|v′−γv +δuv|2+|α′|2+|β′|2+|γ′|2+|δ′|2dt ). En la pr´actica lo que haremos ser´a definir una nube de puntos {si}1≤i≤men [a, b] y resolver m´ın u,v,α,β,γ,δ (1 n n X i=1 |u(ti)−ui|2+|v(ti)−vi|2 +A m m X j=1 |u′(+αu +βv|2+|v′−γv +δuv|2+|α′|2+|β′|2+|γ′|2+|δ′|2t=sj ). 112 Tambi´en se puede a˜nadir la restricci´on de que los coeficientes de la ecuaci´on sean no negativos usando por ejemplo otra penalizaci´on del tipo B m m X j=1 (α−+β−+γ−+δ−)t=sj con B > 0, pero esto complica los c´alculos y quita libertad (el sistema predador presa no deja de ser un modelo matem´atico y quiz´as el sistema en consideraci´on se aproxima mejor por un sistema diferencial donde alguno de los coeficientes es negativo). Observar que con la propuesta que hacemos no obtenemos valores constantes de α, β, γ yδsino funciones. Esto se puede salvar por ejemplo tomando como valor resultante para estas constante las medias en el intervalo [a, b] de las funciones obtenidas. A fin de realizar un experimento num´erico, hemos actuado como sigue: Hemos resuelto con Matlab el problema        u′=−u+ 0′3uv en (0,1′5) v′= 2v−uv en (0,1′5), u(0) = 2, v(0) = 2′8 % Condiciones iniciales u0 = 2; % Valor inicial de u, puedes cambiarlo v0 = 2.8; % Valor inicial de v, puedes cambiarlo y0 = [ u0; v0 ]; % Intervalo de tiempo para la solucion tspan = [0 1.5]; tiempos =0:1.5/200:1.5; % Definir las ecuaciones como una funcion anonima % y resolver el sistema de ODEs usando ode45 fun = @(t ,y)[-y (1) +0.3* y(1) *y (2) ; 2* y(2) -y (1) *y(2) ]; sol = ode45(fun , tspan , y0); y1 =@(t) deval ( sol ,t ,1) ; y2 =@(t) deval ( sol ,t ,2) ; A partir la soluci´on obtenida, hemos extraido su valor en una partici´on uniforme formada por 10 puntos (i.e. n= 10) y los hemos almacenado en ficheros .txt para su posterior lectura en Python. t_medic = linspace(0,1.5,10); u=y1( t_medic ); 113 v=y2( t_medic ); % Guardar los vectores u y v en archivos separados writematrix (u, ’u_vector . txt ’); writematrix (v, ’v_vector . txt ’); Los valores anteriores se han considerado como observaciones en el experimento realizado usando redes neuronales con el objetivo de conocer los valores de los cuatro par´ametros que contiene el sistema de ecuaciones. Mostramos el c´odigo Python correspondiente a la lectura de las datos #leer mediciones with open(’ u_vector . txt ’) as f1: u=[] for line in f1: u. append ( float( line )) with open(’ v_vector . txt ’) as f2: v=[] for line in f2: v. append ( float( line )) Para la implementaci´on del problema de optimizaci´on con TensorFlow, utilizaremos una arquitectura con 4 capas ocultas de 30 neuronas cada una e incluiremos 6 salidas en la capa de salida (las correspondientes a u,vy los cuatro par´ametros) del siguiente modo # Definir la red neuronal class PINN ( tf. keras . Model ): def __init__ ( self ): super(PINN , self ). __init__ () #3n+(k -1)(nˆ2+ n) self . dense1 = tf . keras . layers . Dense (30 , activation = ’tanh ’) self . dense2 = tf . keras . layers . Dense (30 , activation = ’tanh ’) self . dense3 = tf . keras . layers . Dense (30 , activation = ’tanh ’) self . dense4 = tf . keras . layers . Dense (30 , activation = ’tanh ’) self . dense5 = tf . keras . layers . Dense (6 , activation = None) def call (self , inputs ): x = self . dense1 ( inputs ) x = self . dense2 (x) x = self . dense3 (x) x = self . dense4 (x) 114 output = self . dense5 (x) return output A continuaci´on, tomando m= 150, generamos los datos de entrenamiento y creamos el modelo de forma an´aloga al resto de experimentos realizados. Mostramos el c´odigo de los residuos # Definir la funcion de perdida def loss_medic ( model ,T): res = model (T) u=res [: ,0] v=res [: ,1] return tf . reduce_mean ( tf. square (uui)+tf. square (vvi)) def loss_param ( model , X): with tf. GradientTape ( persistent = True ) as tape : tape . watch (X) res = model (X) u=res [: ,0] v=res [: ,1] alpha = res [: ,2] beta =res [: ,3] gamma = res [: ,4] delta = res [: ,5] uder = tape . gradient (u, X) vder = tape . gradient (v, X) alphader = tape . gradient ( alpha , X) betader =tape . gradient ( beta , X) gammader = tape . gradient ( gamma , X) deltader = tape . gradient ( delta , X) del tape return tf . reduce_mean (tf. square ( uder + alpha *u - beta *u*v) +tf. square ( vder - gamma *v+ delta *u*v ) +tf. square ( alphader )+ tf. square ( betader ) +tf. square ( gammader ) +tf. square ( deltader )) Con A= 2 y 10000 iteraciones, la evoluci´on del funcional loss y las funciones obtenidas tras entrenar el modelo son las siguientes: 115 Figura 5.48: Funci´on ual cabo de 104iteraciones. Figura 5.49: Funci´on val cabo de 104iteraciones. Aunque como hemos mencionado el valor del funcional se ha dividido por m´as de 2 ·10−4, los valores correspondientes a los par´ametros est´an muy lejos de los deseados, como ya nos ocurri´o con las ANN. La ventaja es que en este caso podemos iterar m´as porque las iteraciones son menos costosas. Concretamente tenemos µ∼2′5141, ν ∼1′9699, α ∼0′4853, β ∼0′0139, γ ∼1′5135, δ ∼0′7594. El valor del funcional a minimizar en estos puntos es similar al caso ´optimo. Sin embargo, los par´ametros toman valores muy distintos. Podr´ıamos tambi´en pensar que el m´etodo se est´a acercando a otro m´ınimo local, sin embargo, veremos al seguir iterando que esto no parece ser as´ı. Al cabo 122 de otras 10000 iteraciones, el funcional coste desciende a 2′28 ·10−5. Las gr´aficas correspondientes aparecen expuestas Figura 5.50: Funci´on ual cabo de 2 ·104iteraciones. Figura 5.51: Funci´on val cabo de 2 ·104iteraciones. Los valores de las variables son ahora µ∼2′5097, ν ∼1′9917, α ∼0′5169, β ∼0′0333, γ ∼1′8681, δ ∼0′9350. Siguen por tanto lejos de los valores deseados pero se han acercado, en especial los valores de γyδcomienzan a ser similares. Al cabo de otras 10000 iteraciones, el funcional desciende a 7′96 ·10−6. Las gr´aficas son ahora 123 Figura 5.52: Funci´on ual cabo de 3 ·104iteraciones. Figura 5.53: Funci´on val cabo de 3 ·104iteraciones. Los valores de los par´ametros vienen dados por µ∼2′5076, ν ∼1′9980, α ∼0′5441, β ∼0′0491, γ ∼1′9728, δ ∼0′9867. Los valores de αyβsiguen por tanto muy lejos de los deseados aunque las curvas obtenidas por el sistema correspondiente a estos valores est´an muy cerca de los observados. Las aproximaciones de los par´ametros correspondientes a la ecuaci´on para v′s´ı podemos considerarlas como aceptables. Los errores son del orden de una cent´esima. Los valores del funcional siguen descendiendo a medida que iteramos pero muy lentamente. Al cabo de otras 60000 iteraciones tan solo desciende a 3′46 ·10−6. 124 Las gr´aficas son esencialmente las mismas que antes y por tanto no las mostramos. Respecto a los valores de los par´ametros, tenemos µ∼2′5050, ν ∼2′0006, α ∼0′6698, β ∼0′1185, γ ∼2′0143, δ ∼1′0073. Vemos por tanto que aunque el funcional ya pr´acticamente no desciende, los valores de αyβ, aunque siguen lejos de los deseados, s´ı que proporcionan aproximaciones mejores. No est´a, por tanto, claro que el algoritmo se est´e estabilizando en torno a un punto de m´ınimo distinto, sino que puede que est´e convergiendo al correcto pero muy lentamente. El valor de la derivada del funcional coste en la aproximaci´on obtenida es |∇F| ∼ 2′0127 ·10−4. Teniendo en cuenta que la diferencia entre los par´ametros obtenidos y los de partida (en lo cuales la derivada se anula) es del orden de 0′3, por el teorema del valor medio, esto indica la aparici´on de derivadas segundas que toman valores muy peque˜nos, del orden de 0′001, lo cual dificulta la convergencia del m´etodo del gradiente. Finalmente, aportamos el c´odigo Matlab de este experimento: %Generacion de datos de observacion ep =0.01; n =10; al =1; be =0.3; ga =2; de =1; fun =@(t ,y)[y(1) *( - al+ be*y(2) );y (2) *(ga - de*y (1))]; y0=[2.5;2]; sol = ode45(fun ,[0 ,1] , y0 ); tob =0:1/ n :1; uob = deval (sol ,tob ,1) ;vob = deval (sol ,tob ,2) ; %Resolucion del problema inverso al =0; be =0; ga =0; de =0; u0 =0; v0 =0; for i=1:10000 fun =@(t ,y)[y(1) *( - al+ be*y(2) );y (2) *(ga - de*y (1))]; y0 =[ u0 ;v0 ]; es = ode45(fun ,[0 ,1] , y0); u=@(t) deval (es ,t ,1) ;v=@(t) deval (es ,t ,2) ; fuea =@(t,y) [al*y(1)+v(t)*(- be*y (1) +de*y(2) );-ga*y(2)+u (t)*( - be*y (1)+de*y(2) )]; uo=u(tob );vo=v( tob ); sum((uo - uob ) .ˆ2+( vo - vob ) .ˆ2) /(2*( n +1) ) ea0 =[ uob(n+1) -uo(n +1) ;vob (n +1) -vo(n +1) ]; daF=0;dbF=0;dgF=0;ddF=0; for j=n: -1:1 ea = ode45(fuea ,[ tob (j +1) ,tob (j)], ea0 ); phi =@(t) deval (ea ,t ,1) ;psi =@(t) deval (ea ,t ,2) ; tt =tob (j):( tob (j +1) -tob (j)) /100: tob (j+1) ; 125 daF = daF + trapz(tt ,u(tt).* phi (tt)); dbF =dbF - trapz(tt ,u( tt).* v(tt).* phi (tt)); dgF =dgF - trapz(tt ,v( tt).* psi (tt)); ddF = ddF + trapz(tt ,u(tt).*v( tt) .* psi (tt)); ea0 =[ phi(tob (j))+ uob (j) -uo(j); psi (tob (j))+ vob (j)-vo (j)]; end dmF =u0 -uob (1) -phi(tob (1) ); dnF=v0 -vob (1) -psi(tob (1)); u0=max(u0 -ep*dmF ,0);v0=max(v0 -ep*dnF ,0) ;al=max(al -ep* daF ,0); be=max(be -ep*dbF ,0);ga=max(ga -ep*dgF ,0) ;de=max(de -ep* ddF ,0); tt=0:1/100:1; end 126 6. Conclusiones y trabajos futuros Tras la realizaci´on de la presente memoria, concluimos que los m´etodos de resoluci´on de problemas diferenciales bajo condiciones de contorno y/o iniciales utilizando redes neuronales artificiales, se basan en su reducci´on a problemas de optimizaci´on. En particular, utilizamos modelos de redes neuronales informadas por la F´ısica o PINNs, de modo que resolvemos un problema de ajuste de par´ametros para obtener los valores de los pesos y b´ıas de la red que minimizan el error cuadr´atico medio cometido respecto a las restricciones del problema en cuesti´on (ecuaciones diferenciales y las condiciones de contorno y/o iniciales). El gran problema de estos m´etodos es que se resuelven aplicando m´etodos de descenso por gradiente cuya convergencia al m´ınimo global del problema no est´a asegurada debido a la ausencia de condiciones de convexidad adecuadas. Sin embargo, las librer´ıas espec´ıficas para la optimizaci´on con ANN como TensorFlow incorporan herramientas heur´ısticas para evitar m´ınimos locales. Asimismo, estos m´etodos tienen una convergencia lenta y costosa. Respecto al coste computacional, comentamos que, una vez entrenado el modelo podemos almacenarlo como una funci´on susceptible a ser evaluada en todo su dominio. Sin embargo, observamos puntos positivos como la versatilidad de estos modelos para tratar con problemas que requieren bucles iterativos adicionales para los m´etodos de resoluci´on cl´asicos, como los problemas no lineales y los evolutivos. Asimismo, al tratarse de modelos de optimizaci´on, podemos utilizarlos para resolver problemas de control e inversos incluyendo la ecuaci´on diferencial en el funcional coste de forma an´aloga. Desde el punto de vista del An´alisis Num´erico, no hemos encontrado una gran variedad de resultados acerca de las estimaciones de error cometidas por el m´etodo. En aquellos casos en los que se alcanza el m´ınimo global, existen cotas sobre el n´umero de pesos, neuronas y capas a escoger para dise˜nar una arquitectura de red neuronal que nos permita aproximar funciones con un error prefijado, aunque no son optimales. Para trabajos futuros, ser´ıa conveniente estudiar c´omo mejorar dichas cotas, as´ı como dar resultados te´oricos acerca del error cometido teniendo en cuenta la no convergencia del m´etodo e incorporando herramientas estad´ısticas siguiendo la l´ınea de razonamiento usada para los experimentos de la ´ultima secci´on del Cap´ıtulo 4. Para concluir, hay que comentar que la realizaci´on de este trabajo me ha permitido adquirir los conocimientos necesarios para seguir investigando en este tema con un enfoque dirigido a una tesis doctoral, adem´as de formular algoritmos para resolver problemas y realizar experimentos cuyos resultados han sido presentados en el CEDYA / CMA 2024. 127 7. Bibliograf´ıa [1] P.G. Ciarlet,Introduction `a l’analyse num´erique matricielle et `a l’optimisation, Masson, Paris, 1990. [2] D.J. Matich Tema: Redes Neuronales: Conceptos B´asicos y Aplicaciones., apuntes del Departamento de Ingenier´ıa Qu´ımica de la Facultad Regional Rosario, 2001. https://www.frro.utn.edu.ar/repositorio/catedras/quimica/ 5 anio/orientadora1/monograias/matich-redesneuronales.pdf [3] Tema: Redes Neuronales, apuntes del Departamento de Matem´atica Aplicada de la Universidad Polit´ecnica de Madrid para la asignatura Introducci´on al Aprendizaje Autom´atico, 2021 https://dcain.etsin.upm.es/∼carlos/ bookAA/introAA.html [4] M. Hjorth-Jensen Applied Data Analysis and Machine Learning. Solving Differential Equations with Deep Learning, m´odulo de un curso jupyter notebook del Departamento de F´ısica de la Universidad de Oslo. https://compphysics.github.io/MachineLearning/doc/LectureNotes/ build/html/chapter11.html [5] E. Zamora G´ omez Curso de Redes Neuronales Artificiales - Hackeando Tec, curso gratuito de youtube: https://www.youtube.com/playlist?list= PLIyIZGa1sAZo eY8PpuTxfLsja iyytSE [6] E. S´ ANCHEZ KARHUNEN,Deep Learning, seminario para la asignatura Ampliaci´on de Inteligencia Artificial del Departamento de Ciencias de la Computaci´on e Inteligencia Artificial de la Universidad de Sevilla, curso 2021-2022. [7] Y. Lecun, L. Bottou, G. Orr, K.R. M¨ uller,Efficient BackProp. Cap´ıtulo del libro Neural Networks: Tricks of the Trade, Springer, 1998. [8] I.E. Lagaris, A. Likas, D.I. Fotiadis,Artificial Neural Networks For Solving Ordinary And Partial Differential Equations, IEEE Transactions on neural networks, vol. 9, nº5, 1998 . [9] B. Robitaille, B. Marcos, M. Veillette, G. Payre,Quasi-Newton methods for training neural networks, Transactions on Information and Communications Technologies vol 2, 1993. [10] D. Gilbarg, N.S. Trudinger. Elliptic partial differential equations of second order. Springer-Verlag, Berlin 2001. [11] P. Grsvard. Elliptic problems in nonsmooth domains. Classics in Applied Mathematics, SIAM, 2011. [12] D. Yarotsky,Error bounds for approximations with deep ReLU networks, Elsevier Neural Networks, vol. 94, 103-114, 2017 . 128 [13] H.N. Mhaskar,Approximation properties of a multilayered feedforward artificial neural network, Advances in Computational Mathematics, vol. 1, 61 - 80, 1993 [14] R.A. DeVore, R. Howard, C. Micchelli,Optimal nonlinear approximation, Manuscripta Mathematica, vol. 63, 469 - 478, 1989 . [15] M. Anthony, P.L. Bartlett,Neural Network Learning: Theoretical Foundations, Cambridge University Press, New York, ch. 3 y 6, ed. 2009. [16] J. B. D´ ıaz,Las matem´aticas en el pa´ıs de los datos (i): De puntos a m´onadas, La Gaceta de la RSME, 20, 113–142, 2017. [17] M. Telgarsky,Representation Benefits of Deep Feedforward Networks, Cornell University, 2015 . [18] Documentaci´on oficial de numpy: https://numpy.org/doc/stable/ [19] Documentaci´on oficial de matplotlib: https://matplotlib.org/stable/ [20] Documentaci´on oficial de tensorflow: https://www.tensorflow.org/api docs 129