Resolución numérica de ecuaciones diferenciales con redes neuronales y su implementación en Python

Universidad de Sevilla Facultad de Matem´aticas Trabajo Fin de Estudios Grado en Matem´aticas Resoluci´on num´erica de ecuaciones diferenciales con redes neuronales y su implementaci´on en Python Numerical approximation of differential equations with neural networks and their implementation in Python Carlota Gordillo ´ Alvarez Tutor Anna Doubova Krasotchenko Dpto. Ecuaciones Diferenciales y An´alisis Num´erico. Abstract Differential equations and systems play a fundamental role in the world around us, allowing the description and modelling of diverse phenomena with origins in the Applied Sciences. However, most differential equations cannot be solved exactly, which highlights the need to develop new numerical methods to facilitate their solution. In this context, PINNs (Physics-Informed Neural Networks) are emerging as an innovative solution within the field of deep learning. PINNs are based on combining the learning capabilities of neural networks with the structure of differential equations. Throughout this work, we will explain in detail the operation of neural networks, as well as the theoretical basis for the approximation of differential equations based on the Universal Approximation Theorem, and we will illustrate it with the practical resolution of certain examples of them. We will also show its implementation in “Python” using the PyTorch library. Resumen Las ecuaciones y sistemas diferenciales desempe˜nan un papel fundamental en el mundo que nos rodea, permitiendo la descripci´on y modelizaci´on de diversos fen´omenos con origen en las Ciencias Aplicadas. Sin embargo, la mayor´ıa de las ecuaciones diferenciales no pueden resolverse de forma exacta, lo que pone de manifiesto la necesidad de desarrollar nuevos m´etodos num´ericos que faciliten su soluci´on. En este contexto, las PINNs (Redes Neuronales Informadas por la F´ısica) est´an emergiendo como una soluci´on innovadora dentro del campo del Deep Learning. Las PINNs se basan en combinar las capacidades de aprendizaje de las redes neuronales con la estructura de las ecuaciones diferenciales. A lo largo de este trabajo, explicaremos con detalle el funcionamiento de las redes neuronales, as´ı como los fundamentos te´oricos para la aproximaci´on de ecuaciones diferenciales basados en el Teorema Universal de Aproximaci´on, y lo ilustraremos con la resoluci´on efectiva de ciertos casos de las mismas. Tambi´en, mostraremos su implementaci´on en “Python” utilizando la librer´ıa PyTorch. ´ Indice Introducci´on 5 1. Conceptos previos sobre Redes Neuronales 7 1.1. Concepto de una red neuronal ............................. 8 1.1.1. Ejemplo simple de una red de tres capas ................... 11 1.2. Teorema de Aproximaci´on Universal ......................... 13 1.2.1. Conceptos previos ............................... 14 1.2.2. Principales resultados ............................. 15 1.2.3. Aplicaci´on a Redes Neuronales ........................ 18 1.3. Hiperpar´ametros de la red neuronal .......................... 20 1.3.1. Funci´on de activaci´on ............................. 20 1.3.2. N´umero de capas ocultas y neuronas por capa ................ 23 1.3.3. Inicializaci´on de los par´ametros ........................ 24 1.3.4. Funci´on de p´erdida ............................... 25 1.3.5. Optimizador y learning rate .......................... 26 1.4. PINNs .......................................... 31 1.5. Backpropagation .................................... 33 1.5.1. Ejemplo simple de implementaci´on del algoritmo. .............. 36 2. Redes Neuronales para el ajuste de EDO’s 39 2.1. Caso de una EDO: exponencial ............................ 39 2.1.1. Estructura de la red .............................. 41 2.1.2. N´umero de iteraciones ............................. 44 2.1.3. Funci´on de activaci´on. ............................. 45 2.1.4. N´umero de datos de entrenamiento ...................... 46 2.1.5. Optimizadores ................................. 48 2.2. Caso de un SDO: modelo Presa-Depredador ..................... 50 3. T´ecnicas de aproximaci´on de EDO’s: PINNs 53 3.1. Caso de una EDO ................................... 53 3.2. Caso de un SDO: un modelo de epidemias, SIR ................... 58 4. Resoluci´on de Ecuaciones en Derivadas Parciales (EDP’s). 63 4.1. Ecuaci´on del calor ................................... 63 5. Implementaci´on con Python 69 5.1. PyTorch ......................................... 69 1 5.2. Elecci´on de odeint ................................... 70 5.3. Ejemplo paso a paso .................................. 71 5.3.1. Caso de una EDO ............................... 71 5.3.2. Caso de un sistema de EDOs ......................... 78 5.3.3. Caso de una EDP ............................... 79 Bibliograf´ıa 83 2 ´ Indice de figuras 1.1. Red Neuronal multicapa. La figura se encuentra en [8]. ............... 7 1.2. Estructura de una neurona en una red neuronal. .................. 9 1.3. Estructura de una red neuronal con una capa intermedia .............. 10 1.4. Red neuronal con una ´unica capa oculta del Ejemplo 1.3. ............. 11 1.5. Funci´on sigmoidal discontinua ............................. 19 1.6. Funci´on sigmoide .................................... 21 1.7. Funci´on tangente hiperb´olica. ............................. 21 1.8. Funci´on ReLU. ..................................... 23 1.9. Funci´on Swich para distintos valores de β....................... 23 1.10. Fen´omenos de overfitting yunderfitting........................ 24 1.11. Descripci´on del descenso del gradiente. ........................ 27 1.12. Problemas con la elecci´on del learning rate...................... 27 1.13. Comparaci´on de la convergencia entre las variantes del descenso del gradiente, de izquierda a derecha: Batch gradient descent, SGD y Mini-batch gradient descent, [36]. ........................................... 29 2.1. Ajuste del problema (2.3). ............................... 42 2.2. Resultados con una estructura simple ([1,1,1]) y los hiperpar´ametros referidos en la Tabla 2.1. ....................................... 43 2.3. Resultados con una estructura compleja ([1,50,50,50,1]) y los hiperpar´ametros referidos en la Tabla 2.1. ................................ 43 2.4. Resultados con 1000 iteraciones y los hiperpar´ametros referidos en la Tabla 2.1. . 44 2.5. Resultados con 1000000 Iteraciones y los hiperpar´ametros referidos en la Tabla 2.1. 45 2.6. Resultados con la funci´on de activaci´on ReLU y los hiperpar´ametros referidos en la Tabla 2.1. ....................................... 46 2.7. Resultados con la funci´on de activaci´on sigmoide y los hiperpar´ametros referidos en la Tabla 2.1. ..................................... 46 2.8. Datos de entrenamiento tomados para el ajuste (Ndata = 1) ............ 47 2.9. Resultados con Ndata = 1 y los hiperpar´ametros referidos en la Tabla 2.1. . . . . 47 2.10. Datos de entrenamiento tomados para el ajuste (Ndata = 101) ........... 48 2.11. Resultados con Ndata = 101 y los hiperpar´ametros referidos en la Tabla 2.1. . . . 48 2.12. Resultados con el optimizador SGD y los hiperpar´ametros referidos en la Tabla 2.1. ............................................ 49 2.13. Resultados con el optimizador Adagrad y los hiperpar´ametros referidos en la Tabla 2.1. ............................................ 49 3 2.14. Ajuste del problema (2.6). ............................... 52 3.1. Resoluci´on mediante PINNs del problema (3.5) con los hiperpar´ametros referidos en la Tabla 3.1. ..................................... 55 3.2. Resoluci´on mediante PINNs del problema (3.9) con lo hiperpar´ametros referidos en la Tabla 3.2. ..................................... 57 3.3. Modelo SIR. ....................................... 58 3.4. Resoluci´on mediante PINNs del problema (3.10) con los hiperpar´ametros referidos en la Tabla 3.3 ..................................... 61 4.1. Elecci´on de los puntos de colocaci´on. ......................... 65 4.2. Evoluci´on de error durante el entreamiento de la red para la resoluci´on mediante PINNs del problema (4.2) con los hiperpar´ametros referidos en la Tabla 4.1. . . . 66 4.3. Resoluci´on mediante PINNs del problema (4.2) con los hiperpar´ametros referidos en la Tabla 4.1. ..................................... 67 4 Introducci´on Las ecuaciones y sistemas diferenciales tienen un papel fundamental en la descripci´on de fen´omenos naturales del mundo que nos rodea. Asimismo, tienen gran importancia pr´actica en campos muy diversos, que van desde la F´ısica y la Ingenier´ıa hasta la Biolog´ıa y la Econom´ıa, pues permiten modelar fen´omenos con origen en Ciencias Aplicadas. Como, por ejemplo, el movimiento de los fluidos, la din´amica de poblaciones, o incluso, ayudan analizar una enfermedad. Es bien conocido que la mayor´ıa de las ecuaciones y sistemas diferenciales no se pueden resolver de manera exacta. Tradicionalmente, para su resoluci´on num´erica se han utilizado m´etodos cl´asicos como el M´etodo de las Diferencias Finitas o de Elementos Finitos, ver [37], [28]. Estos m´etodos se basan en la construcci´on de una malla que proporciona el “soporte” a la soluci´on del problema en cuesti´on donde se buscan sus valores aproximados. Cuando nos enfrentamos a problemas con una alta dimensi´on, 2 o m´as, la construcci´on y el uso de una malla fina que garantice una buena aproximaci´on de la soluci´on, supone un reto en los m´etodos num´ericos cl´asicos debido al elevado coste computacional y a la dificultad a la hora de su implementaci´on. En las ´ultimas d´ecadas, debido a su gran desarrollo, la Inteligencia Artificial ha comenzado a jugar un papel fundamental para la resoluci´on de ecuaciones diferenciales. Dentro de este amplio campo, concretamente en la rama del Deep Learning (Aprendizaje Profundo), nacen las Physics-Informed neural networks (PINNs, Redes Neuronales Informadas por la F´ısica) como una poderosa herramienta para su resoluci´on. Las PINNs se basan en combinar la capacidad de aprendizaje de las redes neuronales con la estructura de las ecuaciones diferenciales. Esta t´ecnica permite aproximar sus soluciones a trav´es del entrenamiento de las redes neuronales mediante un proceso de optimizaci´on, incorporando la ecuaci´on diferencial y condiciones adicionales (de contorno, condiciones iniciales), directamente en la funci´on de minimizaci´on. Esto permite obtener una aproximaci´on adecuada del problema planteado. De esta forma, conseguimos disminuir el coste computacional de los m´etodos num´ericos cl´asicos, proporcionando una t´ecnica con muy buenas propiedades de generalizaci´on, pues los ´unicos datos que hacen falta para su resoluci´on son los que intervienen en la definici´on del problema en s´ı. Para poder poner en marcha el funcionamiento de las redes neuronales, hace falta una buena elecci´on de ciertos par´ametros. A pesar de los avances de la inteligencia artificial, no existe de 5 Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 12 Ejemplo 1.3 Sea φ∈ N tal que φ:R2−→ R, es decir, una red neuronal con dos datos de entrada y uno de salida. La capa intermedia tendr´a dos neuronas, por tanto d= 2. La funci´on φdada por (1.2) vendr´a representada por [(A1, b1),(A2, b2)], donde A1∈R2×2,b1∈R2yA2∈R1×2,b2∈R. Para simplificar los c´alculos tomaremos como funci´on de activaci´on f(x) = x2, a pesar de que las funciones de activaci´on deben comunmente cumplir la propiedad de ser sigmoidales, ver Secci´on 1.3.1. Una funci´on sigmoidal es aquella que en su recorrido toma valores entre [0,1]. En secciones posteriores lo explicaremos con m´as detalle. Nuestros datos de entrada ser´an: x="x1 x2#="2 2#. Los datos asociados a la primera capa, (A1, b1): A1="w(1) 1,1w(1) 1,2 w(1) 2,1w(1) 2,2#="1 1 2 0 #, b1="b(1) 1 b(1) 2#="1 2#. Y a la segunda capa, (A2, b2): A2=hw(2) 1,1w(2) 1,2i=h1 1 i, b2=hb(2) 1i=h3i. Hemos denotado w(k) i,j el peso asociado a la i-esima neurona de la capa k-esima de la red neuronal con la j-esima neurona de la siguiente capa y b(k) ies la bias asociada a la i-esima neurona de la k-esima capa. Como podemos observar en la Figura 1.2,x2no est´a conectado con la segunda neurona de la capa intermedia, por ello el elemento w(1) 2,2= 0. Como hemos visto anteriormente, para obtener el dato de salida ytenemos que aplicar la siguiente expresi´on: y=A2f(A1x+b1) + b2. Hag´amoslo por partes: 1. El paso por la primera capa: A1x+b1="w(1) 1,1w(1) 1,2 w(1) 2,1w(1) 2,2#" x1 x2#+"b(1) 1 b(1) 2#="1 1 2 0 #" 2 2#+"1 2#="4 4#+"1 2#="5 6#. 2. Aplicaremos la funci´on de activaci´on, f(x) = x2, que se aplicar´a a cada elemento de la matriz Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 13 f(A1x+b1) = f("5 6#) = "f(5) f(6) #="(5)2 (6)2#="25 36 #. 3. El paso por la ´ultima capa: A2f(A1x+b1) + b2=hw(2) 1,1w(2) 1,2i"25 36 #+hb(2) 1i=h1 1 i"25 36 #+h3i= = [ 61 ] + h3i=h64 i. Finalmente, el dato de salida obtenido es ypred = 64. 1.2. Teorema de Aproximaci´on Universal Dado que el objetivo de este trabajo es usar redes neuronales en la aproximaci´on de funciones, en esta secci´on demostraremos que cualquier funci´on continua de nvariables reales con soporte en el hipercubo unidad puede ser aproximada por una composici´on de funciones de combinaciones lineales finitas. Este resultado muestra que podemos aproximar cualquier funci´on continua a trav´es de una red neuronal con solo una capa oculta y utilizando como funci´on de activaci´on una funci´on no lineal sigmoidal. La aproximaci´on de la funci´on deseada ser´a una combinaci´on lineal de la forma: N X j=1 αjσ(yt jx+θj).(1.4) En el contexto de tratado, esta combinaci´on lineal representa el paso por dicha red neuronal. En este caso, denotamos por xlos datos de entrada, y∈Rnyαson los pesos asociados a las conexiones de la capa de entrada con la capa intermedia y la capa intermedia con la capa de salida respectivamente, y θes la bias asociada a la capa intermedia. Estos tres ´ultimos son datos fijados de antemano o actualizados por la propia red. Por ´ultimo, σrepresenta a la funci´on de activaci´on. Seguiremos la notaci´on y los argumentos del art´ıculo de Cybenko [10]. La demostraci´on de este resultado se basa en el hecho de que las sumas de la forma (1.4) son densas en el espacio de las funciones continuas del hipercubo unidad si σes una funci´on sigmoidal, que definiremos a continuaci´on (ver [10]). Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 14 1.2.1. Conceptos previos En este apartado mencionaremos algunos conceptos de an´alisis funcional imprescindibles para entender los resultados que se nombran a continuaci´on. Las demostraciones de los siguientes resultados los podemos encontrar en [11]. Definici´on 1.4 Sea Vun espacio vectorial real. Se dice que una funci´on f:V−→ Res un funcional sublineal osubnorma si cumple las siguientes propiedades: 1. f(x+y)≤f(x) + f(y)∀x, y ∈V. 2. f(αx) = αf(x)∀x∈Vy∀α∈R, con α > 0. Teorema 1.5 (Teorema de Hahn-Banach) Sea V un espacio vectorial real y M un subespacio de V. Supongamos que p es un funcional sublineal sobre V y que f:M−→ Res una aplicaci´on lineal tal que f(x)≤p(x)para todo x∈M. Entonces existe una aplicaci´on lineal g:V−→ Rtal que g|M=fyg(x)≤p(x)para todo x∈V. Recordemos un corolario muy ´util cuando se trata de probar que un subespacio vectorial es denso. Corolario 1.6 Sea F⊂Eun subespacio vectorial tal que ¯ F=E. Entonces ∃f∈E′,f= 0 tal que ⟨f, x⟩= 0 ∀x∈F, donde ⟨·,·⟩ es el producto escalar en la dualidad E′, E. En consecuencia, para demostrar que un subespacio vectorial F⊂Ees denso, se considera una forma lineal y continua fsobre Etal que f= 0 sobre Fy se prueba que fes id´enticamente nula sobre E. Teorema 1.7 (Teorema de la convergencia dominada de Lebesgue) Sea E⊆Rnmedible- Lebesgue. Sea {fn}n∈Nuna sucesi´on de funciones medibles de E en R,guna funci´on integrable no negativa en E, y f:E−→ R, de modo que: 1. Para n∈N,|fn(x)| ≤ g(x)para casi todo x∈E. 2. Para casi todo x∈E,f(x) = l´ımn→∞ fn(x). Entonces fes integrable y ZE fdµ = l´ım n→∞ ZE fndµ. Sean X, Y dos espacios normados. Consideramos T:X−→ Yaplicaciones lineales y continuas (tambi´en llamados operadores lineales continuos). El conjunto L(X, Y ) de los operadores lineales continuos T:X−→ Yes un espacio vectorial para las operaciones habituales. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 15 En este espacio vectorial la cantidad ∥T∥L(X,Y ):= sup v∈X,v=0 ∥T∥Y ∥v∥X = sup ∥v∥X=1 ∥Tv∥Y es una norma. As´ı, se puede hablar del espacio normado de los operadores lineales continuos de Xen Y, denotado de nuevo por L(X, Y ). Se observa que si Yes un espacio de Banach (normado completo), entonces L(X, Y ) tambi´en lo es. En caso particular Y=R, el espacio de Banach L(X, Y ) se denota X′y se denomina dual topol´ogico (o simplemente dual) de X. Si X=Hes un espacio de Hilbert (espacio normado completo, cuya norma es inducida por un producto escalar), entonces Hpuede identificarse algebraica y topol´ogicamente con su dual. Es lo que dice el teorema de representaci´on de Riesz que se presenta a continuaci´on. Teorema 1.8 (Teorema de representaci´on de Riesz) Sea H un espacio de Hilbert, con producto escalar (. , .)Hy norma ∥·∥Hasociada. Dada f∈H′, existe un ´unico uf∈Htal que f(v) = (uf, v)H,∀v∈H. Adem´as, la aplicaci´on f∈H′7→ uf∈Hes un isomorfismo isom´etrico, es decir, es una aplicaci´on lineal, continua, biyectiva y conserva la norma ∥f∥H′=∥uf∥H,∀f∈H′. 1.2.2. Principales resultados Sea In= [0,1]nel hipercubo unidad n-dimensional yC(In) el espacio de las funciones continuas en In. Dada f∈C(In), una funci´on continua en In, denotaremos a la norma del supremo de fcomo ||f||. Se denotar´a M(In) el conjunto de las medidas de Borel finitas. El principal resultado que vamos a presentar en este apartado pretende establecer bajo qu´e condiciones una suma de la forma G(x) = N X j=1 αjσ(yt jx+θj),(1.5) es densa en C(In). Definici´on 1.9 Diremos que una funci´on σes discriminadora si para una medida µ∈ M(In) la igualdad ZIn σ(ytx+θ)dµ(x)=0,∀y∈Rn, θ ∈R implica que µ= 0. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 16 Definici´on 1.10 Decimos que una funci´on σ:R7→ Res sigmoidal si σ(t)→     1 si t→+∞, 0 si t→ −∞. (1.6) Teorema 1.11 Sea σ:R7→ Rcualquier funci´on continua y discriminadora. Entonces las funciones de la forma G(x) = N X j=1 αjσ(yt jx+θj), son densas en C(In)respecto de la norma del supremo. En otras palabras, dada cualquier f∈ C(In)yϵ > 0, existen α, y, θ ∈Rntales que G(x) = G(x;α, y, θ), de la forma (1.5)satisface |G(x)−f(x)|< ϵ, ∀x∈In. Demostraci´on. Sea S⊂C(In) el conjunto de funciones de la forma G(x) como hemos visto en (1.5). Est´a claro que Ses un subespacio vectorial de C(In). Por definici´on, Ses denso en C(In) si ¯ S=C(In), donde ¯ Ses la clausura de S. Vamos a razonar por reducci´on al absurdo. Supongamos que ¯ S⫋C(In), es decir, que la clausura de S no est´a totalmente contenida en C(In). Entonces, la clausura de S, que denotaremos por R=¯ Ses un subespacio propio cerrado de C(In). Por el Teorema de Hanh-Banach, (Corolario 1.6) sabemos que hay un funcional lineal acotado llamado L, que verifica la propiedad de que L= 0, pero L(R) = L(S) = 0. Por el Teorema de Representaci´on de Riesz, este funcional lineal acotado, L, es de la forma L(h) = ZIn h(x)dµ(x) = 0 ∀h∈C(In) para alg´un µ∈ M(In). En particular, como σ(ytx+θ) est´a en R∀y, θ, entonces se tiene que ZIn σ(ytx+θ)dµ(x)=0 ∀y, θ. Luego, como tenemos como hip´otesis que σes discriminadora y, por definici´on sabemos que esta condici´on implica que µ= 0, luego L= 0. Llegamos a una contradici´on, ya que se deber´ıa verificar que L= 0. Por tanto, concluimos que ¯ S=C(In), es decir, el subespacio Sdebe ser denso en C(In). Esto termina la demostraci´on. ■ Este teorema afirma que las sumas de la forma (1.5) son densas en C(In) siendo σuna funci´on continua y discriminadora. Para aplicar este resultado a nuestro marco de trabajo, demostraremos ahora que cualquier funci´on continua sigmoidal σ, de la forma descrita anteriormente (1.6), es discriminadora. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 17 En el contexto de las redes neuronales, las funciones de activaci´on continuas sigmoidales son normalmente tomadas con un crecimiento mon´otono, pero la monoton´ıa no es requerida en estos resultados. Lema 1.12 Cualquier funci´on medible sigmoidal acotada σes discriminadora. En particular, cualquier funci´on sigmoidal continua es discriminadora. Demostraci´on. Sea σuna funci´on sigmoidal continua, es decir, una funci´on continua que cumple la propiedad (1.6). Para demostrar que σes discriminadora notaremos que para cualquier x, y ∈ Rn,θ, φ ∈Rtenemos que σ(λ(ytx+θ) + φ)                −→ 1para ytx+θ > 0cuando λ −→ +∞, −→ 0para ytx+θ < 0cuando λ+−→ ∞, =σ(φ)para ytx+θ= 0 ∀λ. Luego, las funciones σλ(x) = σλ(λ(ytx+θ)+σ) convergen puntualmente y de forma acotada a la funci´on γ(x) =                1para ytx+θ > 0, 0para ytx+θ < 0, σ(φ)para ytx+θ= 0, cuando λ−→ +∞. Sea Πy,θel hiperplano definido por {x|ytx+θ= 0}yHy,θel semiespacio abierto definido por {x|ytx+θ > 0}. Entonces, por el Teorema de la convergencia dominada de Lebesgue, tenemos que 0 = ZIn σλ(x)dµ(x) = ZIn γ(x)dµ(x) = σ(φ)µ(Πy,θ) + µ(Hy,θ)∀φ, θ, y. Ahora probaremos que la medida de todos los semiplanos es 0 implica que la medida µdebe ser, entonces, 0. Fijamos y. Para una funci´on medible acotada, h, definimos el funcional lineal, F, de la siguiente forma F(h) = ZIn h(ytx)dµ(x). Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 18 Se observa que Fes un funcional acotado en L∞(R) ya que µes una medida finita con signo. Sea hla funci´on indicador del intervalo [θ, ∞), es decir, h(x) =        1si x ≥θ, 0si x < θ. En este caso, h(ytx) =        1si ytx−θ≥0, 0si ytx−θ < 0. En particular, hes sigmoidal. Por tanto, nos queda F(h) = ZIn h(ytx)dµ(x) = µ(Πy,−θ) + µ(Hy,−θ)=0. De forma similar obtenemos que F(h) = 0 si hes la funci´on indicador en el intervalo (θ, ∞). Por linealidad, F(h) = 0 para la funci´on indicador en cualquier intervalo y, por tanto, para cualquier funci´on simple (i.e, que sea suma de funciones indicadoras de intervalos). Como las funciones simples son densas en L∞(R) tenemos entonces que F= 0. En particular, sea s(u) = sen (m·u) y c(u) = cos (m·u), dos funciones medibles acotadas, tenemos que F(s+ic) = ZIn cos (mtx) + isin (imtx)dµ(x) = ZIn exp(imtx)dµ(x)=0 ∀m. Luego, la transformada de Fourier de µes 0 y esto implica que µdebe ser 0. Por tanto, σes discriminadora. ■ 1.2.3. Aplicaci´on a Redes Neuronales En este apartado, aplicaremos los resultados probados en el apartado anterior a la teor´ıa de redes neuronales. Demostraremos que a trav´es de una red neuronal con una sola capa interna y una funci´on de activaci´on sigmoidal cualquiera podemos aproximar funciones continuas combinando el Teorema 1.11 y el Lema 1.12 vistos en la secci´on anterior. La precisi´on de esta aproximaci´on va a depender del n´umero de neuronas y par´ametros de la propia red. Teorema 1.13 Dada σcualquier funci´on continua sigmoidal. Entonces las sumas finitas de la forma G(x) = N X j=1 αjσ(yt jx+θj), Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 19 son densas en C(In). En otras palabras, dada cualquier f∈C(In)yϵ > 0, entonces hay una suma, G(x), de la forma (1.5), tal que |G(x)−f(x)|< ϵ ∀x∈In. Demostraci´on. Basta aplicar el Teorema 1.11 y el Lema 1.12, pues cualquier funci´on sigmoidal es discriminadora. ■ Hasta ahora hemos considerado que la funci´on de activaci´on, σ, es una funci´on continua sigmoidal. Pero, tambi´en, la funci´on de activaci´on puede ser discontinua. Como veremos en secciones posteriores, en la pr´actica este tipo es el m´as usado. Consideraremos la funci´on sigmoidal discontinua de la forma (ver Figura 1.5) σ(x) =        1si x ≥0, 0si x < 0. Figura 1.5: Funci´on sigmoidal discontinua En el resultado siguiente asumiremos que σes una funci´on medible acotada sigmoidal. En estas condiciones, tenemos el siguiente resultado an´alogo al Teorema 1.13: Teorema 1.14 Sea σuna funci´on medible, acotada y sigmoidal. Entonces la suma finita de la forma G(x) = N X j=1 αjσ(yt jx+θj), es densa en L1(In). Es decir, dada una funci´on cualquiera f∈L1(In)yϵ > 0, existe una suma, G(x)tal que Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 20 ∥G(x)−f∥L1=ZIn|G(x)−f(x)|dx < ϵ. La prueba sigue la prueba de los Teoremas 1.11 y1.13 sustituyendo las funciones continuas C(In) por funciones integrables, L1(In), y usando el hecho de que L∞(In) es el dual de L1(In). La noci´on de ser discriminadora cambia de acuerdo a los siguientes cambios: para h∈L∞(In) la condici´on de ZIn σ(ytx+θ)h(x)dx = 0 ∀y, θ implica que h(x) = 0 ∀x. Ya hemos visto en el Lema 1.12 que las funciones sigmoidales son discriminadoras. En conclusi´on, hemos demostrado que cualquier superposici´on finita fija de una funci´on discriminadora puede aproximar uniformemente cualquier funci´on continua de nvariables reales con soporte el hipercubo unidad [0,1]n. Tanto las funciones sigmoidales continuas como discontinuas son, com´unmente usadas en el marco de las redes neuronales, son discriminadoras. La uni´on de todos estos resultados demuestra que, a trav´es de una red neuronal con una sola capa interna y una funci´on de activaci´on sigmoidal puede aproximar uniformemente una funci´on continua. 1.3. Hiperpar´ametros de la red neuronal 1.3.1. Funci´on de activaci´on Definiremos a la funci´on de activaci´on f:R−→ Rcomo una funci´on continua que a priori podr´ıa ser tanto lineal como no lineal. Esta funci´on es la encargada de determinar el valor de salida de cada neurona de la red. En este caso, el uso de funciones de activaci´on no lineales juega un papel fundamental debido a la importancia que tiene a la hora de aproximar funciones no lineales. Este tipo de funciones inducen la no linealidad, es decir, si trabaj´aramos con funciones de activaci´on lineales, la red neuronal ser´ıa equivalente a un modelo lineal lo que dificultar´ıa la minimizaci´on del error de la funci´on aproximada y la deseada. Y por consiguiente, la aproximaci´on de dicha funci´on adecuadamente. Tambi´en estas funciones facilitan la derivaci´on recursiva, factor muy importante ya que permite la retropropagaci´on del error y el ajuste de los pesos de la red durante el entrenamiento, a trav´es del algoritmo backpropagation. Por el teorema de Aproximaci´on Universal sabemos que las funciones de activaci´on que requerimos en este caso son aquellas que cumplan la propiedad de ser sigmoidal (1.6). Entre las m´as comunes destacamos las funciones que siguen. 1. Funci´on sigmoide. Es una de las primeras funciones de activaci´on utilizadas en redes neuronales. Suelen ser usadas en problemas de clasificaci´on. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 21 Se define como σ:R−→ [0,1], σ(x) = 1 1 + e−x, (ver Figura 1.6) y su derivada es: dσ dx =e−x (1 + e−x)2=σ(x)(1 −σ(x)). Podemos observar que la derivada puede escribirse en t´erminos de la propia funci´on, esto Figura 1.6: Funci´on sigmoide facilitar´a el c´alculo de las derivadas y a su vez, el aprendizaje de la red. 2. Tangente hiperb´olica. Esta funci´on es similar a la sigmoide, la diferencia es que su recorrido est´a entre [-1,1]. Se define como f:R−→ [−1,1], f(x) = tanh(x) = ex−e−x ex+e−x, Figura 1.7: Funci´on tangente hiperb´olica. y su derivada es: Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 28 imento del par´ametro de entrenamiento (decay rate). Esta t´ecnica se basa en utilizar un learning rate variable donde en las primeras iteraciones sea grande y conforme se vaya avanzando en las iteraciones disminuya. Los algoritmos m´as comunes son el decaimiento inverso (inverse decay) y el decaimiento exponencial (exponential decay). El learning rate de la iteraci´on t,ηt, puede expresarse como sigue: ηt=η0 1 + k·t(decaimiento inverso), ηt=η0exp (−k·t) (decaimiento exponencial), donde η0es el learning rate incial, kes un par´ametro que controla el decaimiento y tes el n´umero de iteraciones realizadas, es decir, n´umero de veces que los datos han pasado a trav´es de la red. Una vez explicado el concepto de learning rate, nos centraremos en el proceso de actualizaci´on de los par´ametros, es decir el tipo de optimizador utilizado. Hay tres variantes del descenso del gradiente: Batch gradient descent, descenso del gradiente estoc´astico (SGD) y mini-batch gradient descent (ver en [33]). Se diferencian entre ellas por la cantidad de datos que usan a la hora de calcular el gradiente de la funci´on objetivo. Batch gradient descent. Calcula el gradiente de la funci´on objetivo para todos los datos en cada actualizaci´on. Se define como θ=θ−η∇θJ(θ), donde θdenota los par´ametros a actualizar y ηes el learning rate. Este m´etodo tiene garantizada la convergencia cuando trabajamos con funciones de p´erdida convexas. Por tanto, se alcanzar´a el valor ´optimo de θen un n´umero finito de iteraciones. Al trabajar con toda la muestra, este m´etodo ser´a computacionalmente redundante para grandes conjuntos de datos. Descenso del gradiente estoc´astico (Stochastic gradient descent, SGD). A diferencia del m´etodo anterior, actualiza el par´ametro para una ´unica muestra del conjunto de datos, (x(i), y(i)). Se puede expresar como θ=θ−η∇θJ(θ;x(i);y(i)), donde θdenota los par´ametros a actualizar y ηes el learning rate. En la pr´actica, la muestra se elige de forma aleatoria. Esto hace que las actualizaciones sean m´as frecuentes y as´ı, se alcance m´as r´apidamente la convergencia respecto del m´etodo anterior. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 29 Mini-batch gradient descent. . Este m´etodo lleva a cabo actualizaciones para muestras de tama˜no n+ 1 del conjunto de datos tratado: θ=θ−η∇θJ(θ;x(i:i+n);y(i:i+n)), donde θson los par´ametros a actualizar, ηel learning rate y (x(i:i+n);y(i:i+n)) son el conjunto de puntos de iai+nque conforman la muestra. Normalmente, los grupos tienen un tama˜no entre 50 y 256 muestras. Este es el algoritmo elegido a la hora de entrenar una red neuronal, ya que al escoger conjuntos de datos pueden ser calculados matricialmente de forma m´as r´apida y eficaz. Respecto a los anteriores, este m´etodo reduce la varianza de la actualizaci´on de los par´ametros, dot´andolo de una convergencia m´as estable. Figura 1.13: Comparaci´on de la convergencia entre las variantes del descenso del gradiente, de izquierda a derecha: Batch gradient descent, SGD y Mini-batch gradient descent, [36]. En la Figura 1.13, comparamos la evoluci´on del error tras la aplicaci´on de los m´etodos de optimizaci´on utilizados estudiados anteriormente. Podemos observar como el m´etodo que converge m´as r´apidamente es el mini-batch gradient descent y el que presenta m´as fluctuaciones es el SGD, debido al tama˜no de la muestra utilizada. Como hemos mencionado anteriormente, la elecci´on del learning rate es muy importante para conseguir que el algoritmo converja. Adem´as, conviente tener en cuenta que en nuestro caso, mantener constante el par´ametro de aprendizaje no favorecer´a la convergencia ya que nuestros datos normalmente no ser´an homog´eneos. Ahora, destacaremos algunos optimizadores com´unmente usados en el entrenamiento de una red neuronal. Entre ellos, diferenciamos los que mantienen el learning rate constante, como el Momentum, y los que lo actualizan a lo largo del entrenamiento, como Adagrad yAdam (v´ease [5]). Momentum.Es un m´etodo que ayuda a acelerar la convergencia y disminuir las oscilaciones del SGD. Consiste en vt=γvt−1+η∇θJ(θ), t ≥1, Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 30 θ=θ−vt, donde vtyvt−1son los vectores de actualizaci´on en la iteraci´on tyt−1 respectivamente, γ∈Res una constante dada, ηes el learning rate yθes el vector de par´ametros de la red. La constante γse suele tomar igual a 0.9 o a un valor similar. Adagrad.Es un algoritmo basado en el descenso del gradiente que adapta el learning rate de todos los par´ametros en cada iteraci´on. En la primera iteraci´on, se lleva a cabo una actualizaci´on de cada uno de par´ametros θi usando el learning rate η. Posteriormente, Adagrad usa diferentes learning rate para cada par´ametro θiy tiempo t. Por simplicidad, denotaremos a gt,i al gradiente de la funci´on objetivo del par´ametro θiy el tiempo t: gt,i =∇θtJ(θt,i).(1.10) En cada iteraci´on, Adagrad modifica cada θien el tiempo tcomo sigue θt+1,i =θt,i −η pGt,ii +ϵgt,i,(1.11) donde Gtes una matriz diagonal, que en la posici´on (i, i) tiene la suma de cuadrados de los gradientes correspondientes a θien el tiempo t. Para evitar dividir por cero se le a˜nade el t´ermino ϵ. Como Gtes una matriz diagonal, podemos reescribir la ecuaci´on (1.11) de forma matricial como θt+1 =θt−η √Gt+ϵ⊙gt, donde ⊙es la multiplicaci´on matriz-vector elemento a elemento. Uno de los principales beneficios de este algoritmo es la adaptaci´on del learning rate, normalmente se toma como inciaci´on η= 0.01. Mientras que, su principal inconveniente es la acumulaci´on de los gradientes al cuadrado en su denominador. Esto puede causar que el learning rate se reduzca, e incluso que se convierta en un valor infitesimal y el algoritmo no se actualice. Adaptative Moment Estimation (Adam).Al igual que el anterior es otro m´etodo que calcula el learning rate para cada par´ametro. Para ello, introduciremos los siguientes par´ametros mt=β1mt−1+ (1 −β1)gt, vt=β2vt−1+ (1 −β2)g2 t, Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 31 donde tyt−1 corresponden a la iteraci´on actual y la anterior respectivamente y gtes gradiente de la funci´on objetivo del par´ametro θiy el tiempo t, notaci´on usada anteriormente en (1.10). mtyvtson los estimadores de el primer momento (la media) y el segundo momento (la varianza) de los gradientes respectivamente, de ah´ı el nombre del algoritmo. Para la inicializaci´on tomamos mt, vt= 0. Los hiperpar´ametros β1, β2∈[0,1), controlan las tasas de decaimiento exponencial de mtyvt. Antes de la actualizaci´on de los par´ametros, calculamos los estimadores ˆmty ˆvtde mty vt, de la siguiente forma ˆmt=mt 1−βt 1 , ˆvt=mt 1−βt 2 , donde βt 1, βt 2denotan respectivamente a β1, β2elevado a t, siendo tla iteraci´on en la que nos encontramos. La regla de actualizaci´on quedar´ıa como sigue: θt+1 =θt−η √ˆvt+ϵˆmt, En [5], se sugiere que tomemos como datos β1= 0.9, β2= 0.999 y ϵ= 10−8. Podemos encontrar informaci´on m´as detallada de los optimizadores mencionados anteriormente en [33]. 1.4. PINNs Physics-Informed Neural Networks (PINNs) es una t´ecnica eficiente y directa de Machine Learning para resolver Ecuaciones Diferenciales Ordinarias (EDOs) y Ecuaciones en Derivadas Parciales (EDPs). PINNs aproxima estos problemas a trav´es del entrenamiento de redes neuronales para minimizar la funci´on de p´erdida, que incluir´a tanto las condiciones de contorno e iniciales como la propia ecuaci´on diferencial. A pesar de ser un algoritmo de aprendizaje supervisado, ya que podemos suponer que conocemos la soluci´on del problema a resolver, ya sea por que se pueda calcular expl´ıcitamente o se ha podido aproximar mediante m´etodos num´ericos como por ejemplo el m´etodo cl´asico de Runge Kutta en el caso de una EDO y el m´etodo de diferencias finitas o Elementos finitos en el caso de una EDP. Este algoritmo se plante´o como un problema de aprendizaje no supervisado que no necesite datos conocidos de la propia soluci´on, convirtiendo as´ı el problema diferencial en un problema de optimizaci´on de la funci´on de p´erdida, con el fin de resolver problemas o ecuaciones de los cuales se desconozcan las soluciones y sean inviables de calcular. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 32 Este m´etodo se puede aplicar a las EDPs o EDOs de la forma N[u;λ] = 0 donde N[u;λ] es un operador diferencial, en general no lineal, con posibles par´ametros λ, y u=u(x, t) denota la funci´on u:Rn+1 7→ Rnen el caso de una EDP para x∈Ω⊂Rny t∈I⊂Ry de orden uno, y en el caso de un sistema de EDO‘s, u=u(x) denotar´a la funci´on u:Rn7→ Rncon x∈Ω⊂Rn, La metodolog´ıa de PINN consiste en determinar los par´ametros θde la red neuronal, los pesos y las bias de cada capa, minimizado la funci´on de p´erdida L(θ), es decir resolviendo el problema θ=arg m´ın θL(θ), donde L(θ) = wFLF(θ) + wBLB(θ) + wdLdata(θ).(1.12) Los tres t´erminos de Ldescriben el error de las condiciones de contorno, LB, el error respecto de la EDO o EDP, LF(θ), y la validaci´on de los valores conocidos de la soluci´on, Ldata. Las costantes wF,wB,wdata son los pesos asociados a cada t´ermino, respectivamente. Generalmente, el que m´as valor tendr´a ser´a wF, ya que el principal objetivo es aproximar la ecuaci´on adecuadamente. Para ello, podemos utilizar los errores vistos en la Secci´on 1.3.4. Para explicar el procedimiento que se sigue, consideraremos un problema de contorno para una EDP escrito como:                F(x, t, ux, ut, ..) = 0, x ∈Ω, t ∈[t0, T ], u(x, t) = uΓ(t), x ∈∂Ω, u(x, t0) = u0(x), x ∈Ω, donde t0yTson el tiempo inicial y final respectivamente, u0la ecuaci´on en el tiempo inicial, t0, y uΓla ecuaci´on para los valores de xsobre la frontera Γ. Las variables x∈Rdyt∈Rdenotan las variables espacial y temporal, respectivamente. Para su resoluci´on asumiremos que tenemos un conjunto de Ndata datos conocidos de la soluci´on ui data en diferentes tiempos ti data y puntos espaciales, xi data, es decir, {ti data, xi data, ui data}Ndata i=1 que constuir´an los datos de entrenamiento incluyendo aqu´ı la condici´on inicial. Denotaremos ˆuθ(x, t) a la predicci´on de la soluci´on de la EDP obtenida por la red neuronal. La funci´on correspondiente a Ldata donde estar´an incluidas las condiciones iniciales del problema u0(usando el error cuadr´atico medio, visto en la Secci´on 1.3.4): Ldata(θ) = 1 Ndata Ndata X i=1 ∥ˆuθ(xdata i, tdata i)−ui data∥2, Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 33 donde ∥·∥es la norma eucl´ıdea. De manera similar, definimos LBla funci´on de p´erdida correspondiente a las condiciones de contorno: LB(θ) = 1 Nb Nb X i=1 ∥ˆuθ(xb i, tb i)−ui b∥, al igual que antes, tenemos un conjunto de Nbpuntos conocidos correspondientes a las condiciones de contorno, {ti b, xi b, ui b}Nb i=1. Por ´ultimo, la funci´on de p´erdida referida LFse obtiene a trav´es de la propia ecuaci´on LF(θ) = 1 Nc Nc X i=1 ∥F(ˆuθ(xc i, tc i))∥, que debe ser evaluada en un conjunto de Ncdatos, llamados puntos de colocaci´on, que son puntos elegidos de manera uniforme en el dominio. Una ventaja de la aproximaci´on a trav´es de redes neuronales es la posibilidad de evaluar la propia ecuaci´on diferencial en los puntos de colocaci´on usando la diferenciaci´on autom´atica. Para la diferenciaci´on autom´atica se necesita calcular derivadas con respecto a los par´ametros de la red, por tanto es necesario implementar un procedimiento de optimizaci´on, vistos anteriormente en la Secci´on 1.3.5. En el caso de una EDO, wB= 0 ya que no posee condiciones de contorno, pero s´ı las condiciones iniciales. 1.5. Backpropagation En una red neuronal con una sola capa intermedia, el proceso de entrenamiento es relativamente directo debido a que la funci´on de p´erdida (o error), estudiada en la Secci´on 1.3.4, es vista como una funci´on de los par´ametros, en este caso, de los pesos y las bias, lo que facilita el c´alculo del gradiente. En el caso de las redes neuronales multicapa, el problema es que la funci´on de p´erdida es una complicada funci´on compuesta de los par´ametros a trav´es de todas las capas de la red. Para facilitar el proceso de entrenamiento usaremos el algoritmo del backpropagation. El backpropagation es un algoritmo eficiente que pertenece al campo de la diferenciaci´on autom´atica y se usa para determinar las derivadas parciales de funciones estructuradas en forma de grafo como es el caso de las redes neuronales. Este algoritmo se basa en la regla de la cadena del c´alculo diferencial. En el contexto de las redes neuronales, la regla de la cadena calcula el gradiente de la funci´on de p´erdida en t´erminos de sumas de productos de gradientes locales de los diferentes caminos desde el nodo de la neurona correspondiente hasta el dato de salida. Sin embargo, las sumas tienen un n´umero exponencial de par´ametros, que podemos calcular eficientemente utilizando la Programaci´on Din´amica. El algoritmo del backpropagation como su propio nombre indica es Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 34 una aplicaci´on directa de la programaci´on din´amica, ya que optimiza el funcional, en este caso la funci´on de p´erdida, de atr´as hacia delante. Este contiene dos partes principales: fase hacia delante (forward) y hacia atr´as (backward). La fase hacia delante consiste en calcular los valores de salida y sus derivadas locales en los nodos, y la fase hacia atr´as se basa en acumular el producto de todos las derivadas locales alrededor de todos los caminos desde el nodo determinado hacia el dato (o datos) de salida. M´as concretamente, dichas fases son las que siguen. Fase hacia delante. En esta fase vamos a llevar a cabo el funcionamiento normal de la neurona, es decir, a trav´es de los datos de entrada y unos par´ametros de inicializaci´on (pesos y bias fijados previamente) vamos a recorrer la neurona desde la primera capa hacia la ´ultima a trav´es de las capas intermedias para obtener el dato de salida. Una vez obtenido el valor predecido por la red lo compararemos con el valor deseado utilizando una funci´on de p´erdida determinada previamente, originando as´ı el error. En nuestro caso, el valor predecido ser´a una aproximaci´on de un problema determinado, ypred, y el valor deseado la soluci´on exacta de ´el, y. Fase hacia detr´as. La finalidad principal de esta fase es minimizar la funci´on de p´erdida en funci´on de los par´ametros, ya que nuestro objetivo es conseguir los pesos y las bias que hagan m´ınimo el error cometido en nuestra aproximaci´on. Se llevar´a a cabo a trav´es del c´alculo de su gradiente respecto a los diferentes par´ametros. Para ello usaremos la regla de la cadena en direcci´on hacia atr´as, empezando desde el nodo de salida. Una vez que los gradientes son minimizados, se actualizan los par´ametros de la red neuronal finalizando as´ı el proceso de entrenamiento. Este algoritmo se llevar´a a cabo tantas veces como sean necesarias para conseguir un error suficientemente peque˜no. El algoritmo puede ser resumido por la combinaci´on de cuatro ecuaciones. Para ello, utilizaremos la notaci´on usada en [34]y[25]. Asumiremos que hay Lcapas en nuestra red y cada capa denotada por l, con l= 1, . . . , L. Consideraremos wl j,k el peso asociado a la j-esima neurona de la capa l-esima de la red neuronal con la k-esima neurona de la capa anterior y bl jes la bias asociada a la j−esima neurona de la l-esima capa. El dato de de salida de la neurona jen la capa lpuede expresarse como: al j=σ(zl j),(1.13) donde zl jest´a expresado usando los datos de salida de la capa anterior a(l−1) k: zl j=X k wjl kal−1 k+bl j.(1.14) En la ecuaci´on (1.13), σrepresenta la funci´on de activaci´on. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 35 Por definici´on, la funci´on de p´erdida Ldepende directamente de los datos de salida de la ´ultima capa aL jy como consecuencia de las ecuaciones (1.13) y (1.14), tambi´en de los datos de salida de todas las capas anteriores. Introduciremos una nueva variable δl jexpresada por: δl j=∂L ∂zl j .(1.15) An´alogamente, usando la regla de la cadena y (1.13) podemos definir δl j, como sigue: δl j=∂L ∂zl j =∂L ∂al j ∂al j ∂zl j =∂L ∂al j σ′(zl j),(1.16) donde σ′(zl j) es la derivada respecto a zl jde la funci´on de activaci´on. Ahora, reescribiremos ∂L ∂wl jk en funci´on de t´erminos conocidos, ∂L ∂wjl k =∂L ∂zl j ∂zl j ∂wjl k =δl j ∂ ∂wjl k X k wjl kal−1 k+bl j!=δl jal−1 k.(1.17) De manera similar, ∂L ∂bl j se puede escribir como ∂L ∂bl j =∂L ∂zl j ∂zl j ∂bl j =δl j ∂ ∂bl j X k wjl kal−1 k+bl j!=δl j.(1.18) Como δl jno solo depende de los valores de salida de la capa l, sino tambi´en de las pr´oximas capas, podemos expresarlo en funci´on de la capa l+ 1, obteniendo δl j=∂L ∂zl j =X k ∂L ∂zl+1 k ∂zl+1 k ∂zl j =X k δl+1 k ∂zl+1 k ∂zl j .(1.19) Teniendo en cuenta (1.13) y (1.14), obtenemos zl+1 k=X j wkl+1 jal j+bl+1 k=X j wkl+1 jσ(zl j) + bl+1 k.(1.20) Entonces, la derivada con respecto de zl jpuede ser calculada como ∂zl+1 k ∂zl j =wkl+1 jσ′(zl j).(1.21) Sustituyendo (1.21) en (1.19), tenemos δl j=X k wkl+1 jδl+1 kσ′(zl j).(1.22) Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 36 Por ´ultimo calculamos, la expresi´on δL jpara todas la neuronas de la ´ultima capa, L, usando (1.16) y obtenemos δL j=∂L ∂zL j =∂L ∂aL j ∂aL j ∂zL j =∂L ∂aL j σ′(zL j).(1.23) Una vez calculada la relaci´on (1.23), basta plantear el algoritmo como un problema de Programaci´on Din´amica e ir calculando δl jpara todas las capas hacia atr´as. En conclusi´on, el algoritmo del backpropagation se puede resumir en los siguientes pasos: 1. Activaci´on en la primera capa. Calculamos a1 jen todas las neuronas de la primera capa. 2. Fase hacia delante. Empezando por la primera capa, recorremos toda la red neuronal calculado zlyalpara l= 2, . . . , L. 3. Calculamos el error en la ´ultima capa. Usando la ecuaci´on (1.16), calculamos (1.23) para todas las neuronas de la ´ultima capa. 4. Aplicamos la programaci´on din´amica. Usando la ecuaci´on (1.22), calculamos δl jpara todas las capas aplicando la programaci´on din´amica. 5. Calculamos el gradiente. Usando las ecuaciones (1.18) y (1.17), calculamos ∂L ∂wjl k y∂L ∂bl j . Las cuatro ecuaciones a las que nos referiamos al inicio de la secci´on son (1.17), (1.18), (1.22) y (1.23). Posteriormente, aplicaremos un optimizador de los vistos en la Secci´on 1.3.5 para actualizar los par´ametros de la red. 1.5.1. Ejemplo simple de implementaci´on del algoritmo. Con el fin de entender mejor el algoritmo explicado y la notaci´on utilizada en el apartado anterior, daremos un ejemplo de implementaci´on de ´el. Para ello consideraremos los datos del Ejemplo 1.3. Ejemplo 1.15 En el Ejemplo 1.3 realizamos primera fase del algoritmo, fase hacia delante con los datos siguientes: x="x1 x2#="2 2#. Los datos asociados a la primera capa, (A1, b1): A1="w(1) 1,1w(1) 1,2 w(1) 2,1w(1) 2,2#="1 1 2 0 #, b1="b(1) 1 b(1) 2#="1 2#. y a la segunda capa, (A2, b2): A2=hw(2) 1,1w(2) 1,2i=h1 1 i, b2=hb(2) 1i=h3i. Cap´ıtulo 1. Conceptos previos sobre Redes Neuronales 37 Obtuvimos como resultado, los siguientes valores z(1) ="z(1) 1 z(1) 2#="5 6#,a(1) ="a(1) 1 a(1) 2#="25 36 #,z(2) =h64 i,a(2) =h4096 i. La funci´on de activaci´on aplicada es f(x) = x2, por tanto su derivada es f′(x)=2x. Tomaremos como funci´on de p´erdida L=1 2(y−aL)2, donde yes el valor que deseamos aproximar. Por simplicaci´on, en los c´alculos siguientes, tomaremos y= 4095. Conviene darse cuenta de lo siguiente ∂L ∂aL j =∂ ∂aL j1 2(y−aL)2=−(y−aL). Ahora, apliquemos los pasos vistos en la descripci´on del algoritmo: 1. Calculamos δ2 1, aplicando la ecuaci´on (1.23), δ2 1=∂L ∂z2f′(z2) = −(y−f(z2))f′(z2) = −(y−f(64))f′(64) = −(4095 −4096)128 = 128. 2. Calculamos δ1 1yδ1 2usando (1.22), δ1 1=w12 1δ2 1f′(z1 1)=1·δ2 1f′(5) = 1 ·128 ·10 = 1280. δ1 2=w12 2δ2 1f′(z1 2)=1·δ2 1f′(6) = 1 ·128 ·12 = 1536. 3. Calculamos las derivadas parciales, usando (1.17) y (1.18) ∂L ∂w11 1 =δ1 1x1= 1280 ·2 = 2560, ∂L ∂w11 2 =δ1 1x2= 1280 ·2 = 2560, ∂L ∂w21 1 =δ1 2x1= 1536 ·2 = 3072, ∂L ∂w21 2 =δ1 2x2= 1536 ·2 = 3072, ∂L ∂w12 1 =δ2 1a1 1= 128 ·25 = 3200, ∂L ∂w12 2 =δ2 1a1 2= 128 ·36 = 4608, ∂L ∂b1 1 =δ1 1= 1260, Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 44 Secci´on 1.3.2). Tras los resultados obtenidos, podemos concluir que la elecci´on adecuada de la red es muy importante. Tanto una elecci´on muy simple, como muy compleja puede llevar a no obtener una aproximaci´on adecuada e incluso llevarnos a fen´omenos como underfitting ooverfitting, mencionados en la Secci´on 1.3.2. Adem´as, esto har´a que se incremente el tiempo de ejecuci´on de la realizaci´on. 2.1.2. N´umero de iteraciones En este apartado estudiaremos como var´ıa la aproximaci´on de la ecuaci´on diferencial (2.3) en funci´on del n´umero de iteraciones que realicemos durante el entrenamiento, es decir, el n´umero de veces que pasan los datos de entrenamiento a trav´es de la red. Plantearemos, como en la secci´on anterior, los dos casos extremos: un valor muy peque˜no y muy alto de iteraciones. Las siguientes realizaciones planteadas del problema poseen los mismos hiperpar´ametros referidos en la Tabla 2.1 a excepci´on de las iteraciones. (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.4: Resultados con 1000 iteraciones y los hiperpar´ametros referidos en la Tabla 2.1. En primer lugar, planteamos el caso de un valor muy peque˜no de iteraciones (iter = 1000, una reducci´on considerable frente a la realizaci´on de la aproximaci´on adecuada). En la Figura 2.4 a la izquierda podemos ver como no obtenemos una buena aproximaci´on de la EDO y a la derecha, como el error disminuye de forma lineal a lo largo de las iteraciones, pero sigue siendo del orden de 103. Sin embargo, merece la pena notar la diferencia en el tiempo de ejecuci´on que se ha reducido hasta llegar a 0.96s, que no llega a ser 1s. Consideramos ahora la realizaci´on con un n´umero muy alto de iteraciones (iter = 1000000, 100 veces m´as que en la realizaci´on de la aproximaci´on adecuada). En la Figura 2.5a podemos observar que pr´acticamente obtenemos la misma predicci´on que en la realizaci´on adecuada descrita anteriormente en la Secci´on 2.1. A la parte derecha de la Figura 2.5, podemos observar como el error se estabiliza a partir de un cierto n´umero de iteraciones. Una diferencia importante Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 45 (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.5: Resultados con 1000000 Iteraciones y los hiperpar´ametros referidos en la Tabla 2.1. es el tiempo de ejecuci´on t= 1022.83s, que ha incrementado enormemente debido al elevado n´umero de iteraciones. Luego, la elecci´on de un elevado n´umero de iteraciones hace que obtengamos una predicci´on que bien podr´ıa ser lograda con un n´umero menor de ellas y un coste de tiempo mucho menor al obtenido. 2.1.3. Funci´on de activaci´on. En la Secci´on 1.3.1, vimos la importancia de la elecci´on de una funci´on de activaci´on apropiada a la hora de aproximar funciones. Debido al problema diferencial planteado (2.3), es fundamental la elecci´on de una funci´on de activaci´on no lineal, ya que como podemos ver en la Figura 2.1 su soluci´on exacta es exponencial, es decir, es una funci´on no lineal. En la realizaci´on de aproximaci´on adecuada realizada anteriormente, hemos elegido una funci´on de activaci´on no lineal, tanh. Estudiaremos ahora los resultados tras la aplicaci´on de una funci´on activaci´on lineal, como puede ser ReLU. Al ser una funci´on de activaci´on lineal, podemos observar en la Figura 2.6a como la predicci´on realizada por la red neuronal presenta trazos lineales y no consigue adaptarse uniformemente a la soluci´on exacta. A pesar de ello, obtenemos una buena aproximaci´on de la soluci´on exacta. Basta fijarse en la Figura 2.6b y ver como el error decrece de forma exponencial llegando a ser insignificante. Con respecto al tiempo, obtenemos t= 10.28spr´acticamente el mismo tiempo de ejecuci´on que la realizaci´on adecuada descrita al comienzo de la Secci´on 2.1. Sin embargo, si consideramos otra funci´on de activaci´on no lineal como es Sigmoide, podemos observar en la Figura 2.7 que obtenemos mejores resultados, incluso que en la realizaci´on considerada adecuada y el tiempo de ejecuci´on se ha reducido (t= 9.72s). En la parte derecha de la Figura 2.7 apreciamos como la predicci´on de la red neuronal es pr´acticamente la misma Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 46 (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.6: Resultados con la funci´on de activaci´on ReLU y los hiperpar´ametros referidos en la Tabla 2.1. (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.7: Resultados con la funci´on de activaci´on sigmoide y los hiperpar´ametros referidos en la Tabla 2.1. que la soluci´on exacta de (2.3) y en la parte izquierda, que el error decrece a lo largo de las iteraciones hasta a llegar a ser del orden de 10−2. Tras los resultados obtenidos en esta secci´on, podemos concluir que la elecci´on de la funci´on de activaci´on es fundamental para la obtenci´on de una buena aproximaci´on de nuestro problema. 2.1.4. N´umero de datos de entrenamiento Debido a que en este Cap´ıtulo estudiamos el ajuste de una EDO a trav´es de redes neuronales proporcionando ´unicamente unos datos conocidos de la soluci´on exacta, es de inter´es estudiar como var´ıa la predicci´on obtenida en funci´on de n´umero de datos proporcionados, es decir, del valor de Ndata en este caso. Como en los apartados anteriores, estudiaremos los dos casos extremos: un valor muy peque˜no (Ndata = 1) y otro muy grande (Ndata = 101). Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 47 Figura 2.8: Datos de entrenamiento tomados para el ajuste (Ndata = 1) Las realizaciones planteadas a continuaci´on poseen los mismos hiperpar´ametros mostrados en la Tabla 2.1, a excepci´on del n´umero de datos de entrenamiento, Ndata. Si consideramos Ndata = 1, en este caso el dato de entrenamiento ser´a ´unicamente el dato inicial del problema, como podemos observar en la Figura 2.8. Los resultados obtenidos de esta implementaci´on est´an representados en la Figura 2.9. En su parte izquierda, podemos observar que la predicci´on de la red neuronal se mantiene constante en y= 0, era de preveer, ya que el ´unico dato que le proporcionamos a la red es el dato inicial (y0= 0). A pesar de ello, podemos observar en la parte derecha de la Figura 2.9 como el error decrece r´apidamente, llegando a ser del orden de 10−13. (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.9: Resultados con Ndata = 1 y los hiperpar´ametros referidos en la Tabla 2.1. Si consideramos ahora el otro caso extremo con Ndata = 101, podemos notar en la Figura 2.10 que, pr´acticamente tenemos todos los puntos de la soluci´on exacta. Como era de esperar, conseguimos una predicci´on de la EDO obtenida por la red casi id´entica (lo podemos ver en la Figura 2.11a), siendo este mejor resultado que el obtenido al comienzo de Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 48 Figura 2.10: Datos de entrenamiento tomados para el ajuste (Ndata = 101) la Secci´on 2.1. En cuanto al error, podemos ver en la Figura 2.11b que decrece de forma lineal. llegando a ser del orden de 10−1. En ambas realizaciones, el tiempo de ejecuci´on ha disminuido frente a la realizaci´on considerada adecuada al comienzo de la Secci´on 2.1, llegando a ser t= 9.77s. Podemos concluir que en el problema de ajuste planteado, mientras mayor sea el n´umero de datos de entrenamiento tomados (Ndata, en este caso el n´umero de datos de la soluci´on exacta), mejor ser´a la predicci´on obtenida por la red neuronal. (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.11: Resultados con Ndata = 101 y los hiperpar´ametros referidos en la Tabla 2.1. 2.1.5. Optimizadores En este apartado estudiaremos como var´ıan los resultados tras la aplicaci´on de diferentes optimizadores: SGD yAdagrad (considerados en la Secci´on 1.3.5) y compararemos los resultados obtenidos con la aproximaci´on adecuada llevada a cabo con el optimizador Adam. Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 49 (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.12: Resultados con el optimizador SGD y los hiperpar´ametros referidos en la Tabla 2.1. En la Figura 2.12, podemos observar que obtenemos la misma predicci´on que la realizaci´on adecuada pero, el error cometido es algo mayor que con el optimizador Adam (siendo el de la realizaci´on adecuada de orden de 10−2y el obtenido en esta de 100). (a) Comparaci´on de la predicci´on de la red y la soluci´on exacta de (2.3). (b) Evoluci´on del error a lo largo del entrenamiento. Figura 2.13: Resultados con el optimizador Adagrad y los hiperpar´ametros referidos en la Tabla 2.1. En la aplicaci´on del optimizador Adagrad (Figura 2.13) obtenemos unos resultados bastantes peores a los obtenidos en la Tabla 2.1. A la izquierda en la Figura 2.13a, observamos que la predicci´on obtenida no se parece para nada a la soluci´on exacta. Esto aparece reflejado en la gr´afica de la derecha, obteni´endose un error del orden de 103. Por tanto, los mejores resultados los obtenemos con diferencia con el optimizador Adam, referidos para ´el en la Tabla 2.1. Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 50 2.2. Caso de un SDO: modelo Presa-Depredador Lo visto en el apartado anterior puede extenderse f´acilmente al caso de un sistema de EDO’s. La diferencia entre un sistema y una ´unica EDO estar´a en el n´umero de variables dependientes que tengamos, y en el contexto de las redes neuronales, en el n´umero de nodos de salida. Adem´as, la funci´on de p´erdida tendr´a ahora la diferencia entre los valores conocidos y la aproximaci´on de cada una de las EDO’s que conformen el sistema. Consider´aremos en este caso el modelo Lotka-Volterra (ver en [14], [9]), tambi´en conocido como el modelo Presa-Depredador. Este modelo consiste en un par de ecuaciones diferenciales de primer orden no lineales que describen la din´amica de dos poblaciones en un medio. Denotaremos P=P(t) al n´umero de depredadores en el instante t, y N=N(t) al n´umero de presas en el instante t, y planteamos el siguiente modelo:              P′=−λP(t) + αN(t)P(t), t ∈[t0, tf] N′=γN(t)−βN(t)P(t), P(t0) = P0, N(t0) = N0, (2.4) donde t0, tf∈Rson datos conocidos que representan el tiempo inicial y final, P0, N0>0 son los valores iniciales de las presas y los depredadores, respectivamente y, λ, α, γ, β son par´ametros positivos que representan las interacciones entre las dos especies. M´as concretamente, λdescribe el crecimiento de los depredadores, γel crecimiento de las presas, αcaracteriza el ´exito en la caza, que afecta al depredador y βrepresenta el ´exito en la caza del depredador, que afecta a la presa. Podemos observar que en ausencia de presas, la poblaci´on de depredadores decrece exponencialmente mientras que, en ausencia de depredadores, la poblaci´on de presas crece de manera exponencial. Como la soluci´on exacta para este sistema no lineal no la podemos determinar de forma expl´ıcita, la determinaremos con la funci´on odeint de Python. Denotaremos por φP(t) y φN(t) a la soluci´on exacta de (2.4) en el instante trepresentando los depredadores y las presas, respectivamente. De acuerdo al caso anterior, determinaremos la funci´on de p´erdida de forma similar a (2.2). Denotemos por yP,pred(ti), yN,pred(ti) a la salida de la red, que ser´a la predicci´on del SDO, es decir, las predicciones para las poblaciones de las presas y los depredadores para el valor ti. De esta manera, podemos definir la funci´on de p´erdida como sigue: L(θ) = Ldata(θ) = 1 Ndata Ndata X i=1 ∥yP,pred(ti)−φP(ti)∥2+1 Ndata Ndata X i=1 ∥yN,pred(ti)−φN(ti)∥2,(2.5) Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 51 donde Ndata es el n´umero de datos de entrenamiento proporcionados a la red y ∥·∥hace referencia a la norma eucl´ıdea. Es decir, igual que en el caso anterior, tomamos como funci´on de p´erdida el error cuadr´atico medio entre la soluci´on aproximada por la red y la soluci´on “exacta” en los puntos de entrenamiento. Los datos de entrenamiento proporcionados, {ti}Ndata i=1 , son puntos equidistantes del intervalo [t0, tf]. De esta forma, consideramos el valor de t0en la funci´on de p´erdida, tomando φP(t0) = P0 yφN(t0) = N0en (2.5). Para fijar ideas, consideraremos el siguiente caso particular :              P′=−P(t)+2N(t)P(t), t ∈[0,10], N′= 5N(t)−3N(t)P(t), P(0) = 2, N(0) = 0.5, (2.6) es decir, λ= 1, α= 2, γ= 5, β= 3, t0= 0, tf= 10, P0= 2 y N0= 0.5. Al igual que en el apartado anterior, la inicializaci´on de los par´ametros de la red se ha realizado con Xavier y las gr´aficas de evoluci´on del error a lo largo del entrenamiento ser´an vistas en la escala logar´ıtmica. En la Tabla 2.2 y en la Figura 2.14 quedan recogidos los hiperpar´ametros y resultados obtenidos para la realizaci´on que se ha considerado adecuada. Ndata Estructura Iteraciones Learning rate 21 [1,15,15,2] 28000 0.01 Optimizador Funci´on de activaci´on Error Tiempo Adam tanh 1.7210−535.606 s Tabla 2.2: Hiperpar´ametros utilizados para el ajuste adecuado del problema (2.6). Al estar ante un problema m´as complejo que el problema (2.3) anterior, hemos necesitado una estructura de red m´as compleja [1,15,15,2], dotada de 4 capas, capas donde sus capas intermedias poseen 15 neuronas cada una y por ello, el tiempo se ha incrementado hasta a llegar a 35.606s. Observando la Figura 2.14b podemos apreciar que con apenas Ndata = 15, podemos aproximar perfectamente soluci´on generada por odeint para el SDO (2.6) y en la Figura 2.14c como el error se minimiza adecuadamente. Al igual que en el apartado anterior, podr´ıamos hacer una discusi´on sobre como afectan la variaci´on de par´ametros en la aproximaci´on del problema (2.6), pero con el objetivo de no alargar en exceso la memoria, no lo presentamos aqu´ı. Cap´ıtulo 2. Redes Neuronales para el ajuste de EDO’s 52 (a) Datos de entrenamiento tomados para el ajuste. (b) Comparaci´on de la predicci´on de la red con los hiperpar´ametros referidos en la Tabla 2.2 y la soluci´on exacta del problema (2.6) (en azul y verde P(t) y en naranja y rojo N(t)). (c) Representaci´on de la evoluci´on del error y del learning rate durante el entrenamiento de la red con los hiperpar´ametros referidos en la Tabla 2.2 para el problema (2.6). Figura 2.14: Ajuste del problema (2.6). Cap´ıtulo 3 T´ecnicas de aproximaci´on de EDO’s: PINNs En la Secci´on 1.4 vimos la aplicaci´on de las Physics-Informed Neural Networks (PINNs) a la resoluci´on de Ecuaciones Diferenciales Ordinarias (EDOs). Se trata de una t´ecnica eficiente que se basaba en aproximar los problemas a trav´es del entrenamiento de las redes neuronales de forma que se minimizase la funci´on de p´erdida. A diferencia con el cap´ıtulo anterior, la funci´on de p´erdida incluir´a tanto las condiciones iniciales del problema, como la propia ecuaci´on diferencial. Recordemos que la funci´on de p´erdida ven´ıa dada por la siguiente expresi´on L(θ) = wFLF(θ) + wBLB(θ) + wdLdata(θ),(3.1) donde LF(θ), representa el error respecto de la EDO o sistema de EDO’s, LB(θ) describe el error de las condiciones de contorno (si existen) y Ldata(θ), la validaci´on de los valores conocidos de la soluci´on, que en este caso s´olo ser´an las condiciones iniciales del problema. Las constantes wF, wBywdson los pesos asociados a cada t´ermino, respectivamente. En nuestro caso, wB= 0 ya que las EDO’s de primer orden no poseen condiciones de contorno. A lo largo de este cap´ıtulo estudiaremos con detalle la aplicaci´on de esta t´ecnica a diversos problemas, como ser´a el caso de una EDO y un sistema de EDO’s. 3.1. Caso de una EDO En primer lugar, tomaremos un ejemplo sencillo para resolverlo mediante la t´ecnica de PINNs como es el caso del problema de Cauchy (2.1) visto en el cap´ıtulo anterior. Definimos el siguiente problema como sigue: (y′(t) = λy(t), t ∈[t0, tf], λ ∈R, y(t0) = y0,(3.2) 53 Cap´ıtulo 3. T´ecnicas de aproximaci´on de EDO’s: PINNs 60 En la Tabla 3.3 y en la Figura 3.4 aparecen recogidos los hiperpar´ametros utilizados y los resultados obtenidos para la aproximaci´on que se ha considerado adecuada. De nuevo, hemos utilizado la inicializaci´on Xavier para los par´ametros de la red y la gr´afica de la funci´on de p´erdida ser´a representada en escala logar´ıtmica para facilitar su interpretaci´on. Debido a la complejidad del sistema planteado, un sistema diferencial no lineal formado por tres ecuaciones, se ha necesitado unos hiperpar´ametros m´as complejos que en los casos estudiados anteriormente. La diferencia m´as notable es la estructura de la red necesitada, ([1,32,32,32,3]), con 5 capas y en cada una de sus capas intermedias 32 neuronas. Tambi´en, hemos necesitado 1000 puntos de colocaci´on para el entrenamiento. Debido a ello, el tiempo de ejecuci´on ha sido de 257.64s. En la Figura 3.4a, podemos observar como la predicci´on obtenida por la red, representada mediante l´ıneas discontinuas, es casi id´entica a la obtenida mediante la funci´on odeint de “Python”, representada mediante l´ıneas continuas. Por otra parte, en la Figura 3.4b vemos como el error decrece r´apidamente a lo largo del entrenamiento y en la Figura 3.4c, como el learning rate se mantiene constante. (a) Comparaci´on de la predicci´on de la red con la soluci´on exacta. (b) Representaci´on de la evoluci´on del error en escala logar´ıtmica. (c) Representaci´on de la evoluci´on del learning rate durante el entramiento de la red. Figura 3.4: Resoluci´on mediante PINNs del problema (3.10) con los hiperpar´ametros referidos en la Tabla 3.3 Cap´ıtulo 4 Resoluci´on de Ecuaciones en Derivadas Parciales (EDP’s). Tal y como hemos explicado en el cap´ıtulo anterior para la resoluci´on de EDO’s, podemos extender f´acilmente las PINNs para la resoluci´on de EDP’s. Para ello, deberemos de modificar la funci´on de p´erdida incluyendo tanto la ecuaci´on en derivadas parciales como las condiciones de contorno, que en el caso de las EDO’s no ten´ıamos, y las condiciones iniciales, s´ı se trata de um problema de evoluci´on. A lo largo de este cap´ıtulo, estudiaremos su implementaci´on el caso de la ecuaci´on del calor. 4.1. Ecuaci´on del calor La ecuaci´on del calor es una ecuaci´on en derivadas parciales que describe como se propaga el calor en un medio determinado a lo largo del tiempo. Debido a que el calor tiende a disiparse a medida que pasa el tiempo, esta ecuaci´on se aplica para modelar gran cantidad de fen´omenos “dispersos”, es decir, fen´omenos donde las part´ıculas o variables estudiadas tienden a dispersarse en el medio, como la evoluci´on de la concentraci´on de un producto sometido a una reacci´on qu´ımica o incluso, la evoluci´on de una poblaci´on de c´elulas tumorales (ver [15], [29]). Introduciremos la siguiente notaci´on: sea Ω un abierto conexo acotado de RN(con N= 1,2 ´o 3) cuya frontera ser´a representada por ∂Ω. En este caso, interpretaremos a Ω como el medio donde se encuentran las part´ıculas o variables estudiadas. Tambi´en, consideraremos el intervalo de tiempo (0, T), donde T > 0 es una constante positiva. De esta forma, denotaremos a Q= Ω×(0, T) y Σ = ∂Ω×(0, T ). Sea u:Q7→ Runa funci´on que representa la difusi´on de las part´ıculas en Ω durante el intervalo de tiempo (0, T ), podemos decir que usatisface el siguiente problema de Cauchy- Dirichlet: 63 Cap´ıtulo 4. Resoluci´on de Ecuaciones en Derivadas Parciales (EDP’s). 64                ut−k∆u=f(x, t),(x, t)∈Q, u(x, t) = uΓ(t),(x, t)∈Σ, u(x, 0) = u0(x), x ∈Ω, (4.1) donde k > 0 es una constante dada, f(x, t) es una funci´on conocida y uΓ, u0son dados que representan, respectivamente, las condiciones de contorno e iniciales para u. Se puede demostrar que bajo ciertas condiciones, el problema (4.1) posee una ´unica soluci´on (ver [7]). Para este caso particular, podemos interpretar a f(x, t) como una fuente externa de calor, k∆ucomo la difusi´on de la temperatura, siendo kla cantidad de flujo emitido, utla variaci´on de uen el tiempo, uΓ(t) como la temperatura en el medio sobre Σ y u0(x) como la temperatura de Ω en el instante inicial. Por tanto, escrita de la forma ut=k∆u+f(x, t) podemos decir que la variaci´on de la temperatura a lo largo del tiempo est´a determinada por la difusi´on de ella y su fuente externa asociada en todo (x, t)∈Q. Por simplicidad, para su implementaci´on pr´actica consideraremos el caso de una ecuaci´on del calor unidimensional (N= 1). Por ejemplo, estudiaremos el siguiente problema                ut−uxx = exp(−t)(sin (πx)−π2sin (πx)),(x, t)∈(−1,1) ×(0,1), u(−1, t) = 0 = u(1, t), t ∈(0,1), u(x, 0) = sin (πx), x ∈(−1,1), (4.2) es decir, en este caso tenemos que Ω = (−1,1), T= 1, f(x, t) = exp(−t)(sin (πx)−π2sin (πx)), uΓ= 0 y u0(x) = sin (πx). Es f´acil ver que su soluci´on exacta es u(x, t) = exp(−t) sin (πx),(x, t)∈(−1,1) ×(0,1), pues satisface tanto la ecuaci´on, como las condiciones iniciales y de contorno de (4.2). La funci´on de p´erdida asociada a este problema se puede definir como sigue: L(θ) = 1 Nc Nc X i=1 ∥ut,pred(xc i, tc i)−uxx,pred(xc i, tc i)−exp(−tc i)(sin (πxc i)−π2sin (πxc i))∥2 +1 Nb Nb X i=1 ∥upred(−1, tb i)∥2+∥upred(1, tb i)∥2+∥upred(xb i)−sin (πxb i)∥2!, donde upred es la soluci´on predecida por la red neuronal, ut,pred yuxx,pred representan, respectivamente, la derivada de primer orden respecto de ty de segundo orden respecto de xde la predicci´on de la red neuronal, que ser´a realizada mediante el algoritmo del backpropagation. Cap´ıtulo 4. Resoluci´on de Ecuaciones en Derivadas Parciales (EDP’s). 65 Al igual que en los casos estudiados anteriormente, Ncrepresenta los puntos de colocaci´on donde ser´a evaluada la ecuaci´on en derivadas parciales durante el entrenamiento de la red y Nb los puntos para evaluar las condiciones de contorno e iniciales del problema. Debido a la complejidad del problema, hemos elegido los puntos de colocaci´on para el entrenamiento de la red a trav´es del muestreo del hipercubo latino (Latin Hypercube sampling). El muestreo del hipercubo latino es una t´ecnica estad´ıstica que se basa en generar muestras aleatorias de forma controlada y eficiente, ya que aleatoriamente evita que se repitan puntos. Es una generalizaci´on n-dimensional del “cuadrado latino” en 2 dimensiones, una matriz de n×n elementos en la que cada casilla est´a ocupada por uno de los nelementos de tal modo que cada uno de ellos aparece exactamente una vez en cada columna y en cada fila (ver [42], [24]). Tambi´en exiten otros tipo de elecci´on de los puntos de colocaci´on como la forma aleatoria o los puntos de Sobol (ver [21]). Figura 4.1: Elecci´on de los puntos de colocaci´on. En la Figura 4.1, podemos ver como act´uan las distintas t´ecnicas de muestreo. Claramente, se observa como la t´ecnica m´as eficiente es el muestreo del hipercubo latino (a la derecha), ya que los puntos ocupan pr´acticamente todo el dominio y no hay espacios huecos, como en el muestreo de forma aleatoria (Random). NcNbEstructura Iteraciones Learning rate 10000 100 [2,22,22,1] 26000 0.001, x 0.99/20 iteraciones Optimizador wbFunci´on de activaci´on Error Tiempo Adam 1 tanh 4.210−41275.12 s Tabla 4.1: Hiperpar´ametros utilizados para el ajuste adecuado del problema (4.2). En la Tabla 4.1 y en las Figuras 4.2 y4.3 aparecen recogidos los hiperpar´ametros utilizados y la representaci´on de los resultados obtenidos para la aproximaci´on que se ha considerado adecuada. Al igual que en los casos anteriores, la inicializaci´on de los par´ametros de la red ha sido realizada con Xavier y las gr´aficas de la funci´on de p´erdida y el learning rate ser´an representadas en escala logar´ıtmica con el fin de facilitar su interpretaci´on. A pesar de ser (4.2) una EDP unidimensional, es un problema m´as complejo de resolver que los estudiados anteriormente. Esto aparece reflejado en los resultados obtenidos, principalmente Cap´ıtulo 4. Resoluci´on de Ecuaciones en Derivadas Parciales (EDP’s). 66 en el tiempo de ejecuci´on, t= 1275.12s, pues ha incrementado notablemente respecto a los casos anteriores. Tambi´en, aparece reflejado en el n´umero de puntos de colocaci´on para evaluar la EDP, Nc= 10000, y para evaluar las condiciones iniciales y de contorno, Nb= 100. En este caso, la estructura de la red, ([2,22,22,1]), est´a formada por cuatro capas con dos neuronas en la capa de entrada, pues tenemos dos variables independientes (x, t) en la EDP, 22 neuronas en las capas intermedias y una en la capa de salida. En la Figura 4.3, podemos observar como la soluci´on exacta (a la izquierda) y la predecida por la red (a la derecha) son pr´acticamente id´enticas. (a) Evoluci´on del error. (b) Evoluci´on del learning rate. Figura 4.2: Evoluci´on de error durante el entreamiento de la red para la resoluci´on mediante PINNs del problema (4.2) con los hiperpar´ametros referidos en la Tabla 4.1. A diferencia de los casos anteriores, como podemos ver en la Tabla 4.1 no hemos considerado el learning rate constante, sino variable de forma que decrezca a lo largo del entrenamiento. M´as concretamente, hemos considerado un par´ametro de aprendizaje variable que cada 20 iteraciones lo multipliquemos por 0.99. De forma gr´afica, podemos observar en la Figura 4.2b como decrece de forma lineal conforme avanzan las iteraciones. Tal y como explicamos en la Secci´on 1.3.5, esto nos va a permitir que la red se adapte m´as r´apidamente al problema a lo largo del entrenamiento. (a) Soluci´on exacta. (b) Soluci´on predecida por la red. (c) Soluci´on exacta. (d) Soluci´on predecida por la red. Figura 4.3: Resoluci´on mediante PINNs del problema (4.2) con los hiperpar´ametros referidos en la Tabla 4.1. Cap´ıtulo 5 Implementaci´on con Python En este ´ultimo cap´ıtulo, nos centraremos en la implementaci´on pr´actica de las redes neuronales para la resoluci´on de ecuaciones diferenciales. Hasta ahora, hemos estudiado los fundamentos de las redes neuronales y el Teorema de Aproximaci´on Universal que nos permite utilizarlas para su resoluci´on. Para su implementaci´on, utilizaremos uno de los lenguajes de programaci´on m´as potentes y famosos hoy en d´ıa, “Python” (v´ease [26]). Python es un lenguaje com´unmente usado en el campo del aprendizaje profundo (deep learning), debido a su sintaxis simple y legible, flexibilidad y versatilidad que permite a los desarroladores experimentar f´acilmente con diferentes arquitecturas de modelos, t´ecnicas de entrenamiento y conjuntos de datos. Tambi´en es conocido por ser un lenguaje de prop´osito general, pues destaca por su capacidad de ser usado en una gran variedad de aplicaciones y contextos, como la inteligencia artificial o, incluso el desarrollo de juegos. Adem´as, posee librer´ıas para el aprendizaje profundo como son TensorFlow, Keras y PyTorch (v´ease [6], [27], [1],[18], [36]). A lo largo de este cap´ıtulo, no s´olo presentaremos los conceptos fundamentales de PyTorch, sino que tambi´en proporcionaremos ejemplos pr´acticos y c´odigos sencillos para ilustrar su implementaci´on. Veremos c´omo definir modelos de redes neuronales, c´omo entrenarlos con datos de ecuaciones diferenciales y c´omo utilizarlos para realizar predicciones y resolver problemas espec´ıficos. 5.1. PyTorch PyTorch es una librer´ıa de Deep Learning de c´odigo abierto que se basa en el lenguaje de programaci´on de Python y la librer´ıa Torch. Torch es una librer´ıa de c´odigo abierto que se usa para crear redes neuronales profundas (deep neural networks). PyTorch fue desarollada por Facebook AI Research en 2017 y desde ah´ı no ha dejado de crecer. Al mismo tiempo que se cre´o aparecerieron otras librer´ıas similares como TensorFlow, 69 Cap´ıtulo 5. Implementaci´on con Python 76 4output_size = 1 5learnig_rate = 0.01 6model = PINN ( input_size , hidden_size , output_size ) #[1 ,5 ,1] 7optimizer = optim . Adam ( model . parameters () , lr = learnig_rate ) # optimizador escogido 8 9loss_fn = PINNLoss () # Definicion de la funcion de perdida 10 11 # Puntos de colocacion 12 t_p = torch . linspace (0 ,1 ,10) . view ( -1 ,1) . requires_grad_ ( True ) #10 puntos equidistanciados del dominio 13 14 num_epoch = 10000 # numero de iteraciones 15 16 # Registrar el tiempo de inicio , ya que vamos a llamar a la funcion 17 # para comenzar el entramiento de la red 18 start_time = time .time () 19 20 # Entrenar el modelo utilizando la funcion de entrenamiento 21 loss_history , learning_rates = train ( model , optimizer , loss_fn , t_data , t_p , num_epoch ) 22 23 24 # Registrar el tiempo de finalizacion , puesto que hemos finalizado el entrenamiento 25 end_time = time . time () 26 27 # Calcular el tiempo transcurrido 28 elapsed_time = end_time - start_time 29 print (f' ' ) 30 print (f'Tiempo de ejecucion : { elapsed_time } segundos ') En el siguiente c´odigo, obtenemos los resultados del entrenamiento de la red neuronal (los imprime la funci´on train vista anteriormente). En cada l´ınea nos indica la iteraci´on en la que nos encontramos y el error producido por la red en ella. Al final del c´odigo (l´ınea 7) nos imprime el tiempo de ejecuci´on de la red. 1Epoch [2000/10000] , Total Loss : 0.0018908255733549595 2Epoch [4000/10000] , Total Loss : 0.0010343898320570588 3Epoch [6000/10000] , Total Loss : 0.0007274135132320225 4Epoch [8000/10000] , Total Loss : 0.0004271422512829304 5Epoch [10000/10000] , Total Loss : 0.000105695151432883 6 7Tiempo de ejecucion : 15.440767049789429 segundos Tras este proceso, ya tenemos la red neuronal entrenada. Por tanto, es el momento de obtener la predicci´on de la red para nuestro problema. Para obtener una representaci´on gr´afica de los resultados, seguimos los siguientes pasos: Cap´ıtulo 5. Implementaci´on con Python 77 1. Evaluamos en el modelo: 1# Evaluamos el modelo en todos los datos del dominio , al igual que la solucion exacta 2with torch . no_grad () : 3ypred = model (t) En este caso, tson los puntos del dominio. 2. Representamos las gr´aficas, (estos son los c´odigos correspondientes a la Figura 3.2) 1# Comparacion de la solucion exacta y la prediccion de la red 2plt . figure () 3plt . plot (t , ypred . detach () , c=" black " , label =" Sol. predicta por NN ") 4plt. plot(t,y ,c="blue ", linestyle =" --", label =" Sol . exacta ") 5plt . legend () 6plt . title (f 'Prediccion de la red con n. iter : { num_epoch } ') 7plt. show () 8 9# Graficas de evolucion del error y el learning rate 10 plt . figure () 11 # Representacion de la evolucion de la funcion de perdida a lo largo del entrenamiento 12 plt . plot ( range (1 , len ( loss_history ) + 1) , loss_history , label = 'Training Loss ', color ='blue') 13 plt . xlabel ( 'Epoch ') 14 plt . ylabel ( 'Loss') 15 plt . title ( 'Error a lo largo del entrenamiento ') 16 plt . yscale ( 'log ')# graficas vistas en escala logaritmica 17 plt. show () 18 plt . figure () 19 20 # Representacion de la evolucion del learning rate a lo largo del entrenamiento 21 plt . plot ( range (1 , len ( loss_history ) + 1) , learning_rates , label = 'Learning Rate ', color ='red ') 22 plt . xlabel ( 'Epoch ') 23 plt . ylabel ( 'Learning Rate ') 24 plt . title ( 'Evolucion del Learning Rate ') 25 plt. show () Observaci´on 5.1 Para la resoluci´on de este problema hemos considerado un learning rate constante. Si quisi´eramos utilizar uno que var´ıe a lo largo del entrenamiento, bastar´ıa a˜nadir a la funci´on train mencionada anteriormente el siguiente c´odigo: 1scheduler . step () # paso para actualizar el learning rate 2 e indicar en el c´odigo donde a˜nadimos los hiperpar´ametros de la red, como va a variar. Por ejemplo, Cap´ıtulo 5. Implementaci´on con Python 78 1scheduler = optim . lr_scheduler . StepLR ( optimizer , step_size =20 , gamma =0.99) 2 en este caso cada 20 iteraciones multiplicaremos ηpor 0.99. 5.3.2. Caso de un sistema de EDOs En este apartado indicaremos las modificaciones que habr´ıa que realizar al c´odigo visto en la secci´on anterior para la resoluci´on de un SDO. En particular, consideremos el modelo SIR estudiado en la Secci´on 3.2. Recordemos que este modelo ven´ıa dado por el siguiente sistema:                          S′=−0.1IS, I′= 0.1IS −0.9I, t ∈[0,3] R′= 0.9I, S(0) = 100, I(0) = 50, R(0) = 0. Tal y como hemos mencionado anteriormente, lo ´unico que debemos modificar es la funci´on de p´erdida, pues debemos incluir tanto la definici´on del sistema de ecuaciones diferenciales, como las condiciones iniciales. No obstante, tambi´en cambiar´a nuestra estructura de la red, pues al tener tres ecuaciones tendr´a que tener tres neuronas en la capa de salida (es decir, output-size=3). A diferencia del caso anterior, ahora tenemos un problema con tres variables dependientes y por tanto, deberemos calcular la derivada con respecto a tde cada una de ellas. En el siguiente c´odigo, indicamos cu´al ser´ıa la funci´on encargada de calcular la derivada en este caso, 1# Funcion para calcular la derivada 2def compute_derivative ( model , x): 3g = x. clone () 4g.requires_grad_(True) # Activar el flag requires_grad para que se realice el seguimiento de gradientes 5y_pred = model (g)# Obtener las predicciones del modelo para las entradas 6s= y_pred [: ,0]. view (1000 ,1) # 1000 son los puntos de colocacion 7i= y_pred [: ,1]. view (1000 ,1) 8r= y_pred [: ,2]. view (1000 ,1) 9s_t = torch . autograd . grad (s , g , grad_outputs = torch . ones_like (s) , retain_graph = True, create_graph=True)[0] 10 i_t = torch . autograd . grad (i , g , grad_outputs = torch . ones_like (i) , retain_graph = True, create_graph=True)[0] 11 r_t = torch . autograd . grad (r , g , grad_outputs = torch . ones_like (r) , retain_graph = True, create_graph=True)[0] 12 return s,i,r,s_t , i_t , r_t Cap´ıtulo 5. Implementaci´on con Python 79 Al tener un modelo con 3 neuronas en la capa de salida, cada columna se corresponder´a con una variable dependiente. Por ello, ses primera columna (l´ınea 6), ies segunda columna (l´ınea 7) yres tercera columna (l´ınea 8). Tal y como vimos en la Secci´on 3.2, la funci´on de p´erdida est´a definida como la suma de los errores de cada una de las ecuaciones y condiciones iniciales que conforman el sistema. De esta manera, definimos as´ı el error como aparece entre las l´ıneas 7 y 12. 1class PINNLoss ( nn . Module ): 2def __init__ ( self ): 3super ( PINNLoss , self ). __init__ () 4self . mse_loss = nn. MSELoss () 5 6def forward ( self , model , s_t , i_t , r_t , s, i, r, t0 , y0): 7loss_pde = self . mse_loss (s_t , torch . Tensor (- beta *s*i)) 8loss_pde += self . mse_loss ( i_t , torch . Tensor ( beta *s*i - gamma *i)) 9loss_pde += self . mse_loss ( r_t , torch . Tensor ( gamma *i)) 10 y_0 = initial_condition (model , t0 ) 11 loss_initial_condition = self . mse_loss (y_0 , y0) 12 total_loss = loss_pde + loss_initial_condition 13 return total_loss Podemos observar como en las l´ıneas 7-9, realizamos el error entre la predicci´on de la red y la definici´on de la ecuaci´on para cada una de ellas. Los valores de β, γ son proporcionados como datos al inicio del c´odigo. El c´odigo completo que se ha necesitado para la resoluci´on de (3.10) puede encontrarse en el repositorio [2] con el nombre de SIR-PINNs.ipynb. 5.3.3. Caso de una EDP Al igual que en el apartado anterior, indicaremos las modificaciones a realizar en el c´odigo de la Secci´on 5.3.1. Consideraremos el caso de la ecuaci´on del calor estudiada en la Secci´on 4.1. Recordemos que ven´ıa dada por el siguiente problema de Cauchy-Dirichlet:                ut−uxx = exp(−t)(sin (πx)−π2sin (πx)),(x, t)∈(−1,1) ×(0,1), u(−1, t) = 0 = u(1, t), t ∈(0,1), u(x, 0) = sin (πx)x∈(−1,1). Debido a que en su resoluci´on hemos elegido los puntos de colocaci´on a traves del muestreo por el hipercubo latino, debemos importar, adem´as de todas las librer´ıas consideradas anteriormente, la librer´ıa donde se encuentra. Cap´ıtulo 5. Implementaci´on con Python 80 1# Importar librerias 2from pyDOE import lhs # Latin Hypercube Sampling Debido a que estamos resolviendo una EDP, construiremos una malla para definir el dominio Q. 1x= torch . linspace (x0 ,xf , Nx). view ( -1 ,1) # intervalo de la x ( -1 ,1) 2t= torch . linspace (t0 ,tf , Nt). view ( -1 ,1) # intervalo de tiempo (0, T) 3# Crear el mallado 4X,T= torch . meshgrid (x. squeeze (1) ,t. squeeze (1) ) # construye la malla Al igual que en el cap´ıtulo anterior, debemos modificar tanto la funci´on que se encarga de calcular las derivadas, pues en este caso necesitamos derivadas segundas y respecto a las dos variables independientes, como la funci´on de error, ya que tenemos otra ecuaci´on y en este caso, tambi´en condiciones de contorno. 1# Funcion para calcular la derivada 2def compute_derivative ( model , x): 3g = x. clone () 4g.requires_grad_(True) 5f = model (g) 6f_x_t = torch . autograd . grad (f ,g , torch . ones ([ g. shape [0] ,1]) . to( device ) , retain_graph=True, create_graph=True)[0] # primera derivada 7f_xx_tt = torch . autograd . grad ( f_x_t ,g , torch . ones (g . shape ). to ( device ) , retain_graph=True, create_graph=True)[0] # segunda derivada 8f_t = f_x_t [: ,[1]] 9f_xx = f_xx_tt [: ,[0]] 10 return f_t , f_xx , g 11 12 # Funcion para calcular las condiciones iniciales y de contorno 13 def initial_condition ( model , x_BC ): 14 return model ( x_BC ) En este caso, calculamos las derivadas respecto a las dos variables de forma conjunta. Es decir, f-x-t ser´a un tensor donde en la primera columna almacenar´a la derivada de fcon respecto de x,f-x, y en la segunda columna, f-t. De igual manera, ocurre con f--xx--tt. Por simplicidad, en este caso trabajaremos conjuntamente con las condiciones de contorno e iniciales, luego la funci´on initial-condition calcular´a tanto la predicci´on de la red para ambas condiciones del problema. El siguiente c´odigo define la funci´on de p´erdida: 1class PINNLoss ( nn . Module ): 2def __init__ ( self ): 3super ( PINNLoss , self ). __init__ () 4self . mse_loss = nn. MSELoss () 5 6def forward ( self ,model , x_BC ,u_BC ,f_t ,f_xx ,g): Cap´ıtulo 5. Implementaci´on con Python 81 7h = torch . exp (-g [: ,1]) *( torch . sin ( np. pi*g [: ,0]) - np. pi **2* torch . sin (np .pi *g [: ,0]) ) # fuerza externa 8fh =f_t - f_xx + h # definicion de la edp 9loss_pde = self . mse_loss (fh , f_hat ) 10 u_BN = initial_condition (model , x_BC ) # incluye tantos las condiciones de contorno como las inciales 11 loss_condition = self . mse_loss (u_BN , u_BC ) # condiciones de contorno e iniciales 12 total_loss = loss_pde + loss_condition 13 return total_loss A diferencia del c´odigo proporcionado en la Secci´on 5.3.1, en la definici´on de esta funci´on hemos incluido la definici´on de la EDP, como podemos observar en la l´ınea 8. En este caso, f-hat ser´a un tensor de la misma dimensi´on que fh donde todos sus elementos son 0. Al igual que en los casos anteriores, hemos considerado wF=wB=wd= 1 El siguiente c´odigo definiremos de forma adecuada las condiciones iniciales y de contorno del problema para poder entrenar la red neuronal, 1# Condiciones iniciales 2# X puntos del dominio 3# U valores que toma u 4# u(x ,0) = sin(pi*x) 5X_left = X [: ,0] 6T_0 = T[: ,0] 7left_X =torch . hstack ((X[: ,0][: , None],T [: ,0][: , None ])) # apila de forma horizontal 8left_U = left_Y = torch . sin (np .pi * left_X [: ,0]) . unsqueeze (1) 9 10 # Condiciones de contorno 11 # u(-1,t)=0=u(1,t) 12 # x =-1, primera fila , t( todas ) primera fila 13 bottom_X = torch . hstack (( X[0 ,:][: , None ],T[0 ,:][: , None ])) 14 bottom_U = torch . zeros ( bottom_X . shape [0] ,1) # vector de ceros 15 # x=1, ultima fila 16 top_X = torch . hstack ((X [ -1 ,:][: , None ],T[ -1 ,:][: , None ]) ) # vector de ceros 17 top_U = torch . zeros ( top_X . shape [0] ,1) 18 19 # Lo agrupamos conjuntamente 20 # [ iniciales , contorno ] 21 X_train = torch . vstack ([ left_X , bottom_X , top_X ]) # apila de forma vertical 22 U_train = torch . vstack ([ left_Y , bottom_U , top_U ]) Al haber definido anteriormente X,T como una malla, de tal forma que se definir´an los valores de Ω por filas y los instantes de tiempo, (0, T ), por columnas. La funci´on torch.hstack apila los tensores de forma horizontal, es decir, uno debajo de otro, y la funci´on torch.vstack de forma vertical (ver [27]). Esta ´ultima la utilizaremos para trabajar con las condiciones de contorno e iniciales de forma conjunta, como podemos observar Cap´ıtulo 5. Implementaci´on con Python 82 en las l´ıneas 21 y 22. Tal como explicamos en la Secci´on 4.1, elegiremos los puntos de colocaci´on a trav´es del hipercubo latino. Esto se implementa con el siguiente c´odigo, 1# Puntos de colocacion ( para evaluar la EDP), los elegimos a traves del hipercubo latino . 2# lb = primer valor del dominio ( -1 ,0) 3# ub = ultimo valor del dominio (1 ,1) 4X_train_Nf = lb +( ub - lb)* lhs (2 ,Nf ) # 2, pues tenemos dos variables (x,t) De esta forma elegimos, Nfpuntos de forma aleatoria en el dominio con el fin de evaluar la EDP durante el entrenamiento. Nfha sido definido con anterioridad en el c´odigo. El c´odigo completo para este problema puede encontrarse en el repositorio [2] con el nombre de Heat1D-PINNs.ipynb. Bibliograf´ıa [1] Aston Zhang, Zachary C. Lipton, Mu Li, Alexander J. Smola. Dive into Deep Learning. Junio, 2021. https://arxiv.org/abs/2106.11342. [2] Carlota Gordillo. Numerical aproximation of differential equations with neural networks. Repositorio GitHub, 2024. https://github.com/carlotagordillo2/ Numerical-approximation-of-differential-equations-with-neural-networks.git [3] Charu Aggarwal, Neural Networks and Deep Learning. Springer, 2018. [4] Daniel Bourke. Learn PyTorch for Deep Learning: Zero to Mastery book. Julio 2022. [5] Diederik P.Kingma, Jimmy Lei Ba. Adam: A method for Stochastic Optimization. International Conference on Learning Representations, 2015. [6] Eli Stevens, Luca Antiga, Thomas Viehmann. Deep Learning with PyTorch. Manning, 2020. [7] Evans, L.C. Partial differential equations. Second edition. Graduate Studies in Mathematics, 19. American Mathematical Society, Providence, RI, 2010. [8] Facundo Bre, Juan M.Gimenez. Prediction of wind pressure coefficients on building surfaces using Artificial Neural Networks. Artincle in Energy and Buildings, November 2017. https: //www.researchgate.net/publication/321259051_Prediction_of_wind_pressure_ coefficients_on_building_surfaces_using_Artificial_Neural_Networks. [9] Fred Brauer, Carlos Castillo - Chavez. Mathematical Models in Population Biology and Epidemiology. Springer, 2012. [10] G.Cybenko, Approximation by Superposition of a Sigmoidal Function. Springer, 1989. [11] Ha¨ım Br´ezis, Analyse Fonctionnelle. Th´eorie et applications. Masson, 1987. [12] Hubert Baty, L´eo Baty. Solving differentia equations using physics informed deep learning: a hand-on tutorial with benchmark tests. CERMICS, Ecole des Ponts, 2023. [13] Huili Yu. Weight Inizialization for Deep Neural Network. April, 2023. https://medium.com/@freshtechyy/ weight-initialization-for-deep-neural-network-e0302b6f5bf3. 83 Bibliograf´ıa 84 [14] J.D. Murray. Mathematical Biology. I: An Introduction. Springer, 2002. [15] Jackson, T., Komarova, N., Swanson, K.R., Mathematical oncology: using mathematics to enable cancer discoveries. Amer. Math. Monthly 121 (2014), no. 9, 840–856. [16] Jorge Nocedal, Stephen Wright. Numerical Optimization. Springer, 2006. [17] Juan Diego Toscano. Repositorio “GitHub”. 2024. https://github.com/jdtoscano94/ Learning-Scientific_Machine_Learning_Residual_Based_Attention_PINNs_ DeepONets. [18] Keras documentation. KERAS.https://keras.io/. [19] Maia Martcheva. An Introduction to Mathematical Epidemiology. Springer, 2015. [20] Mart´ın Hern´andez. Repositorio “GitHub”. 2023. https://github.com/DCN-FAU-AvH/ pinns_heat. [21] M. Mckay, R. Beckman and W. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics, vol. 21, no. 2, pp. 239-245, May. 1979 [22] Neha Yadav, Anupam Yadav y Manoj Kumar, An Introduction to Neural Network Methods for Differential Equations. Springer, 2015. [23] NumPy documentation. NUMPY.https://numpy.org/. [24] Olsson, A. et al. (2003). On Latin hypercube sampling for structural reliability analysis. Structural Safety: Volume 25, Issue 1, January, Pages 47-68. Retrieved January 5, 2018. [25] Pankaj Mehta, Ching-Hao Wang, Alexandre G. R. Day, and Clint Richardson. A high-bias, low-variance introduction to Machine Learning for physicists. Boston University, University of California, Unlearn.AI, San Francisco and City University of New York. May, 2019. [26] Python Software Foundation. PYTHON.https://www.python.org/. [27] PyTorch documentation. PYTORCH.https://pytorch.org/. [28] Randall J. LeVeque. Finite Difference Methods for Ordinary and Partial Differential Equations. Siam, 2007. [29] Salsa S., Diffusion. In “Partial Differential Equations in Action”, UNITEXT, vol 86. Springer, Cham, 2015. [30] Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi y Francesco Piccialli. Scientific Machine Learning through Physics-Informed Neural Networks: Where we are and What’s next. University of Naples Federico II, SISSA and University of Colorado, Junio 2022. TFG Carlota Gordillo ´ Alvarez. Grado en Matem´aticas Cap´ıtulo 5. Bibliograf´ıa 85 [31] Savyakhosla. Mini-Batch Gradient Descent with Python. Agosto, 2022. https://www. geeksforgeeks.org/ml-mini-batch-gradient-descent-with-python/. [32] SciPy documentation. SCIPY.https://scipy.org/. [33] Sebastian Ruder. An overview of gradient descent optimization algoritms. Insight Centre for Data Analytics, NUI Galway, 2017. [34] Stefan Kollmannsberger, Davide D’Angella, Moritz Jokeit y Leon Herrmann, Deep Learning in Computational Mechanics. Springer, 2021. [35] Sushant Patrikar. Batch, Mini Batch Stochastic Gradient Descent. Octubre, 2019. https://towardsdatascience.com/ batch-mini-batch-stochastic-gradient-descent-7a62ecba642a. [36] TensorFlow documentation. TENSORFLOW.https://www.tensorflow.org/. [37] Thomas J.R.Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications,INC. New York, 2000. [38] Vaibhav Haswani. Learning Rate Decay and methods in Deep Learning. Septiembre, 2020. https://medium.com/analytics-vidhya/ learning-rate-decay-and-methods-in-deep-learning-2cee564f910b#:~: text=Learning%20rate%20decay%20is%20a,help%20both%20optimization%20and% 20generalization.. [39] V.I.Gorikhovskii, T.O.Evdokimova y V.A.Polentamshy. Neural networks in solving differential equations. Journal of Physics: Conference Series,2022. [40] Wes McKinney. Python for Data Analysis O’Reilly Media, Inc, 2017. [41] William E.Boyce, Richard C.Diprima, Douglas B.Meade. Elementary Differential Equations and Boundary Problems. Wiley, 2017. [42] Xin, L. Numerical Methods for Engineering Design and Optimization: Latin Hypercube Sampling (LHS). 2014