Solving Partial Differential Equations using Artificial Neural Networks

Ph.D. Dissertation Doktorego tesia Solving Partial Differential Equations using Artificial Neural Networks Deribatu partzialetako ekuazioen ebazpena neurona-sare artifizialen bidez Carlos Uriarte 2023 (cc)2024 CARLOS URIARTE BARANDA (cc by-nc-sa 4.0) Ph.D. Dissertation Solving Partial Differential Equations using Artificial Neural Networks Carlos Uriarte Supervisors: David Pardo, Elisabete Alberdi 2023 Doktorego tesia Deribatu partzialetako ekuazioen ebazpena neurona-sare artifizialen bidez Carlos Uriarte Zuzendariak: David Pardo, Elisabete Alberdi 2023 This dissertation is written in English and Basque. In the part prior to the main body, the mixture of languages obeys the following rule: for contents with a single paragraph, two paragraphs will appear—the first one in English and the second one in Basque. For contents with more than one paragraph, they will be distinguished by subsections entitled in either English or Basque. In the body of the dissertation, odd chapters are developed in Basque while even chapters are in English. To internationalize the scope of the presented work, the appendices include the English translations of the odd chapters written in Basque. Tesi hau ingelesez eta euskaraz idatzita dago. Tesiaren gorputzaren aurreko zatian, hizkuntzen nahasketak honako arau hau beteko du: bi paragrafo agertuko dira (lehenengoa ingelesez eta bigarrena euskaraz) testu-paragrafo bakarra duten edukietan. Paragrafo bat baino gehiagoko edukiak aldiz, ingelesezko edo euskarazko izenburua duten azpiatalen bidez bereiziko dira. Tesiaren gorputzean, kapitulu bakoitiak euskaraz garatuko dira, eta kapitulu bikoitiak, ingelesez. Aurkeztutako lanaren irismena nazioartera zabaltzeko xedez, euskarazko kapitulu bakoitien ingelesezko itzulpenak eranskinetan jasoko dira. i Contents/Aurkibidea Acknowledgments/Esker onak ix Abstract/Laburpena xiii Acronyms/Akronimoak xv Chapters/Kapituluak 1 1 Sarrera 3 1.1 Motibazioa ............................... 3 1.2 Zenbakizko metodo tradizionalak .................. 5 1.3 Neurona-sareak hurbiltzaile unibertsal gisa ............. 6 1.4 Berrikuspen literarioa ......................... 8 1.4.1 NNak erabiliz PDEak ebazteko lehenengo lanak ...... 8 1.4.2 NNen bidez PDEak ebazten TensorFlow-tik aurrera . . . . 10 1.5 Ekarpen nagusiak ........................... 12 1.5.1 Deep elementu finituen metodoa ............... 12 1.5.2 Deep Ritz metodo bikoitza .................. 13 1.5.3 Memorian oinarritutako Monte Carlo integrazioa ..... 15 1.6 Egitura ................................. 15 2 Feed-Forward Neural Networks 17 2.1 Architecture .............................. 17 2.2 Continuum, parameterized, and discretized setups ......... 19 2.3 Gradient-based training ....................... 22 2.3.1 Continuum-level Gradient-Descent method ......... 23 2.3.2 Parameterized-level Gradient-Descent method ....... 24 2.3.3 Stochastic Gradient-Descent method ............ 25 2.4 A case of study ............................ 26 Appendix 2.A Interpretation of the set of realizations as a family of vector spaces ............................. 32 iii Acknowledgments/Esker onak enjoyed a lot. I also thank her for introducing me to Prof. Demkowicz, from whom I learned a great deal of knowledge I have applied in my research works. Similarly, I thank Ignacio, Paulina, Keko, and Patrick for their incredible kindness and hospitality during my research stay at the Institute of Mathematics of the Pontifical Catholic University of Valpara´ıso, Chile. I have many good anecdotes engraved in my memory in which they were taking part, particularly at moments of “carreteo”. It is difficult to sufficiently acknowledge all the lived experiences with so many colleagues during these four years. To avoid making this section too long, I briefly thank Jamie, Jon Ander, Julen, Lena, Oscar, Felipe, Mahdi, Ali, Jes´us, Manuela, and Tom´as. Although the greeting is brief, it has a lot of meaning. I thank them all a lot. I would also like to thank all the BCAM staff for assisting me during the bureaucratic difficulties encountered, which have not been few. In this regard, I want to thank Miguel especially. He has been infinitely patient and lovely with me on many occasions. A thousand thanks! Besides, during these four years, I have met people who, although not directly related to my research environment, have become friends of mine due to their relationships with workmates and significantly contributed to the easygoingness of this adventure. I particularly want to thank Claudio and Mario for our many good times together. Finally, I want to thank my family and friends for their support. In particular, I want to offer my eternal gratitude to my father, Enrique; my brother, Luis; and my partner, Iratxe, who unconditionally accompanied me and firsthand experienced my moments of tension and joy. Thank you for understanding the conflict of a researcher lifestyle—working at ungodly hours, leaving aside the usual healthy social habits, simply because I felt I could not procrastinate my moments of research inspiration. Heartfeltly, thanks. x Esker Onak Doktoretza egiten hasi nintzenetik, arrastoa utzi didaten hainbat pertsona ezagutu ditut bidean, eta guztiei adierazi nahiko nieke nire esker ona. Lehenik eta behin, nire esker onik handiena adierazi nahi diot Davidi. David pertsona berezia da, jarrera asertiboa eta axolagabea dauka, eta bizitza ezohiko eran ikusten du. Ezaugarri horiek, haren pazientziarekin, haren trebetasun interpertsonal harrigarriekin, eta bere-berea duen intuizio handiarekin batera, pertsona benetan zoragarria egiten dute berarekin lan egiteko orduan, harengandik aholkuak jasotzeko orduan, eta une dibertigarrietan eta ez hain dibertigarrietan gozatzea posible izaten da. Eskerrak eman nahi dizkiot, hain zuzen ere, haren etengabeko irakaspenengatik, gure elkarrizketa luzeengatik, eta egitea hainbeste kostatu zitzaizkidan lanak desegiteko duen berezko diziplinagatik (egia esan, arrazoia izaten zuen gehienetan). Beti egongo naiz berari eskertuta doktorego-bidaia liluragarri eta errepikaezin honetan aurkitu ditudan arrakasta eta behaztopengatik. Aldi berean, eskerrak eman nahi dizkiot nire beste zuzendariari, Elisabeteri, bere ardurapean tesia garatzeko aukera emateagatik. Batez ere, eskerrak eman nahi dizkiot tesi honen antolaketan eta idazketan egin duen jarraipenagatik. Bera gabe, askoz lan konplexuagoa izango zatekeen. Nire zuzendariei ezezik, laneko giro bikaina lankideei ere zor diet. Zehazki, Javier funtsezko lankidea izan zen nire doktoretzaren lehenengo fasean; beti egoten zen prest pandemian zehar era birtualean laguntzeko. Lan aspergarri eta nekagarrietan duen diziplinak liluratu egiten nau, eta ez ditut ahaztu nahi gogor lan egindako egunen amaieran “azken garagardoa” hartzeko proposamenak (gutxitan azkena izaten zena). Asko ikasi dut harengandik eta une ahaztezinak partekatu ditugu. Anari ere eskerrak eman nahi dizkiot, nire karrera aurredoktoralean hasieratik bukaeraraino ondoan egon den lankidea izan baita, eta une gozagarriak partekatu baititut hitzordu profesional eta ez hain profesional guztietan. Une ez oso profesionalen artean, Donejakue bideko abenturak eta Aste Nagusietako distiraz betetako denboraldiak dira aipagarrienetakoak. Asko eskertzen diot Juditi hainbeste une berezi partekatzegatik nire egonaldian zehar Oden Konputazioko eta Ingeniaritzako Institutuan, Texasko Unibertsitatean Austinen, Estatu Batuetan. Asko gustatu zitzaidan abentura izan zen. Halaber, eskerrak eman nahi dizkiot Demkowicz irakaslea aurkezteagatik, nire ikerlanetako askotan aplikatu dudan ezagutzaren garatzailea hura baita. Era berean, eskerrak eman nahi dizkiet Ignacio, Paulina, Keko eta Patrick ikerlariei, Matematikako Institutuan Valpara´ısoko Pontifizia Katolika Unibertsitatean nire egonaldian izan zuten adeitasun eta abegikortasunagatik. Anekdota on asko grabatuta dauzkat oroimenean, batez ere “carreteo” garaietakoak. Zaila da lau urte hauetan hainbeste lankiderekin bizitako esperientzia guztiak xi Acknowledgments/Esker onak behar bezala aitortzea. Atal hau gehiegi ez luzatzearren, labur bada ere, eskerrak eman nahi dizkiet Jamie, Jon Ander, Julen, Lena, Oscar, Felipe, Mahdi, Ali, Jes´us, Manuela eta Tom´as lankideei. Nahiz eta agurra laburra izan, esanahi handikoa da. Era berean, eskerrak eman nahi dizkiet BCAMeko administrazioko langile guztiei, hain gutxi izan ez diren zailtasun burokratikoen aurrean laguntza emateagatik. Bereziki, Migueli. Oso maitagarria eta pazientzia handikoa izan da nirekin askotan. Mila esker! Gainera, lau urte hauetan, nire ikerketako ingurunearekin lotura zuzenik izan ez arren, egindako adiskideei ere eskerrak eman beharrean nago. Bereziki, Claudiori eta Mariori elkarrekin izan ditugun momentu itzelengatik. Azkenik, eskerrak eman nahi dizkiet nire etxekoei eta adiskideei emandako laguntzagatik. Bereziki, esker onak nire aita Enriqueri, Luis anaiari eta nire bikotekide Iratxeri, baldintzarik gabe lagundu didatelako eta nire tentsioeta poztasun-uneak bertatik bertara bizi izan dituztelako. Hain zuzen, eskerrik asko ikertzaile baten bizimoduak dakarren gatazka ulertzeagatik; hau da, ezorduetan lan eginez, ohiko gizarte-ohitura osasungarriak alde batera utzita, ikerketainspirazioko uneak atzeratzerik sentitzen ez nuelako. Bihotz-bihotzez, eskerrik asko. xii Abstract/Laburpena Abstract Partial differential equations have a wide range of applications in modeling multiple physical, biological, or social phenomena. Therefore, we need to approximate the solutions of these equations in computationally feasible terms. Nowadays, among the most popular numerical methods for solving partial differential equations in engineering, we encounter the finite difference and finite element methods. An alternative numerical method that has recently gained popularity for numerically solving partial differential equations is the use of artificial neural networks. Artificial neural networks, or neural networks for short, are mathematical structures with universal approximation properties. In addition, thanks to the extraordinary computational development of the last decade, neural networks have become accessible and powerful numerical methods for engineers and researchers. For example, imaging and language processing are applications of neural networks today that show sublime performance inconceivable years ago. This dissertation contributes to the numerical solution of partial differential equations using neural networks with the following two-fold objective: investigate the behavior of neural networks as approximators of solutions of partial differential equations and propose neural-network-based methods for frameworks that are hardly addressable via traditional numerical methods. As novel neural-network-based proposals, we first present a method inspired by the finite element method when applying mesh refinements to solve parametric problems. Secondly, we propose a general residual minimization scheme based on a generalized version of the Ritz method. Finally, we develop a memory-based strategy to overcome a usual numerical integration limitation when using neural networks to solve partial differential equations. xiii Abstract/Laburpena Laburpena Deribatu partzialetako ekuazioak aplikazio ugari dituzte fenomeno fisiko, biologiko edo sozial anitzen modelizazioan. Horregatik, funtsezkoa da ekuazio horien soluzioen hurbilpenak konputazionalki egingarriak diren terminoetan adieraztea. Gaur egun, ingeniaritzan, deribatu partzialetako ekuazioak ebazteko zenbakizko metodo ezagunenen artean daude diferentzia finituen eta elementu finituen metodoak. Duela gutxi, deribatu partzialetako ekuazioak ebazteko ospea hartu duen zenbakizko metodo bat neurona-sare artifizialen erabilpena da. Neurona-sare artifizialak edo, laburtzearren, neurona-sareak hurbilketa-propietate unibertsalak dituzten egitura matematikoak dira. Gainera, azken hamarkadako garapen konputazional apartari esker, ingeniarientzako eta ikertzaileentzako zenbakizko metodo eskuragarri eta indartsu bihurtu dira. Adibidez, irudien eta hizkuntzen prozesamendua neurona-sare artifizialen gaur egungo aplikazioak dira, eta duela urte batzuk pentsaezina zen errendimendu bikaina erakusten dute. Tesi honetan neurona-sare artifizialen bidezko deribatu partzialetako ekuazioen zenbakizko ebazpena aztertuko dugu, honako helburu bikoitzarekin: alde batetik, neurona-sareen portaera ikertzea deribatu partzialetako ekuazioen soluzioen hurbiltzaile gisa, eta bestetik, neurona-sareetan oinarritutako metodoak proposatzea, ohiko zenbakizko metodoen bidez nekez ekin dakiekeen lan-esparruetarako. Neurona-sareetan oinarritutako proposamen berritzaile gisa, lehenik eta behin, elementu finituen metodoaren funtzionamenduan oinarritutako metodo bat aurkeztuko dugu, problema parametrikoak ebazteko diskretizazioan fintzeak aplikatzen direnean. Bigarrenik, hondarren minimizazio eskema orokor bat proposatuko dugu, Ritz-en metodoaren bertsio hedatu batean oinarritua. Azkenik, memorian oinarritutako zenbakizko integrazio teknika bat erakutsiko dugu. Teknika horren helburua da deribatu partzialetako ekuazioak ebazteko neurona-sareak erabiltzen direnean agertzen den ohiko muga gainditzea. xiv Acronyms/Akronimoak Throughout this dissertation, English acronyms will be used exclusively, regardless of whether they appear in English or Basque portions of text. In English, the first letter of each constituent word is capitalized. In Basque, only proper or foreign constituent words have the first letter capitalized (or at the beginning of a sentence). In the following list, acronyms appear in singular form in English, while in bare form and between parentheses in Basque.∗ Tesi honetan zehar, ingelesezko akronimoak bakarrik erabiliko dira, horiek agertzen diren testu zatia ingelesez ala euskaraz idatzita dagoen kontuan hartu gabe. Ingelesez, lehen letra larriaren irizpideari eutsiko diogu hitz eratzaile guztietan. Euskaraz, hitz eratzaileen lehen letrak larriz egongo dira izen berezien kasuan bakarrik (edo esaldiaren hasiera bada). Hurrengo zerrendan, akronimoak singularrean agertzen dira ingelesez, eta euskaraz, berriz, oinarriko forman eta parentesi artean.∗ 1D one-dimensional (dimentsio bakar) 2D two-dimensional (bi dimentsio) 3D three-dimensional (hiru dimentsio) AD Automatic Differentiation (Deribaketa automatiko) ANN Artificial Neural Network (Neurona-sare artifizial) BVP Boundary-Value Problem (Mutur-balioen problema) D2RM Deep Double Ritz Method (Deep Ritz metodo bikoitz) DeepFEM Deep Finite Element Method (Deep elementu finituen metodo) DNN Deep Neural Network (Neurona-sare sakon) DRM Deep Ritz Method (Deep Ritz metodo) FDM Finite Difference Method (Diferentzia finituen metodo) ∗With some exceptions/Salbuespenak salbuespen. xv FEM Finite Element Method (Elementu finituen metodo) FFNN Feed-Forward Neural Network (Aurrerantz elikatutako neurona-sare) GDRM Generalized Deep Ritz Method (Deep Ritz metodo orokortua) GPU Graphics Processing Unit (Prozesamendu grafikoko unitate) MC Monte Carlo NN Neural Network (Neurona-sare) ODE Ordinary Differential Equation (Ekuazio diferentzial arrunt) PDE Partial Differential Equation (Deribatu partzialetako ekuazio) PINN Physics-Informed Neural Network (Fisikatik informatutako neurona-sare) SGA Stochastic Gradient Ascent (Gradiente igoera estokastiko) SGD Stochastic Gradient Descent (Gradiente jaitsiera estokastiko) SGDM SGD with Momentum (SGD momentum-arekin) TF TensorFlow TF2 TensorFlow 2 UAT Universal Approximation Theorem (Hurbiltasun unibertsalaren teorema) WANs Weak Adversarial Networks (Sare aurkari ahul [pluralean]) xvi List of Figures/Irudien zerrenda 1.1 Aurreranzko eta alderantzizko problemen diagrama eremu elektromagnetiko batean. .......................... 4 2.1 Graph of a fully-connected FFNN with a three-dimensional input, two-dimensional output, depth three, and hidden layers of width five. .................................. 18 2.2 Some usual activation functions appearance: (a) sigmoid, (b) hyperbolic tangent, (c) rectified linear unit, and (d) softplus. . . . . 19 2.3 Sketch of the parameterization of a NN and its objective-function minimization formulation. ...................... 20 2.4 Gradient-descent performance sketch on a convex scenario. . . . . 24 2.5 Graph of F ◦ΦNN for (2.23). Top panel: 3D view. Bottom panel: Top view. ............................... 28 2.6 Graph of the partial derivatives of F ◦ ΦNN for (2.23). Dashed curves indicate vanishing regions. .................. 29 2.7 Different executions of the gradient-descent method for (2.23). We maintain constant learning rates in all cases and carry out 200 iterations. Top panel: 3D view. Bottom panel: Top view. ..... 31 3.1 FEMaren diskretizazio-sare fintzeen eta DeepFEMaren geruzen arteko lotura. .............................. 36 3.2 Denda itxuradun eta zatika linealak diren oinarriko funtzioen, ψjen, euskarria. ............................. 38 3.3 1Dko diskretizazio-sare baten fintze uniformeen eta nodoen zenbaketaren progresioa. .......................... 39 3.4 Eskema ez-parametrikoaren DeepFEM arkitektura dinamikoaren irudikapena koefiziente konstantedunetarako. Gezi grisek aldagai entrenagarriak adierazten dituzte, eta gezi beltzek, berriz, (entrenaezinak diren) zabaltze-operazioak. ................ 41 3.5 s. urratsean bloke entrenagarri baten arkitekturaren irudikapena. Gezi grisek aldagai entrenagarriak adierazten dituzte. ....... 42 xvii List of Figures/Irudien zerrenda 3.6 Eskema parametrikoaren DeepFEM arkitektura dinamikoaren irudikapena koefiziente konstantedun problemarako. Gezi grisek aldagai entrenagarriak adierazten dituzte, eta gezi beltzek, berriz, (entrenaezinak diren) zabaltze-operazioak. ............. 42 3.7 Eskema parametrikoaren DeepFEM arkitektura dinamikoaren irudikapena zatika konstanteak diren koefizienteak dituztenetarako. Gezi grisek aldagai entrenagarriak adierazten dituzte, eta gezi beltzek, berriz, (entrenaezinak diren) zabaltze-operazioak. ...... 43 3.8 (3.25). problemaren lehenengo lau urratsetako DeepFEMeko iragarpenak. u∗soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. .............. 50 3.9 Hasieratik amaierarainoko entrenamendua DeepFEMarentzako (3.25). probleman. Galera-funtzioa errorearen energia-norma da. ..... 52 3.10 Geruzaz geruzako entrenamendua DeepFEMean galera-funtzioan energia-norma erabiltzen denean (3.25). probleman. ........ 53 3.11 DeepFEMaren entrenamenduak (3.25). problemarentzat galera-funtzioan L2-norma erabiltzean. ......................... 54 3.12 Hasieratik amaierarainoko DeepFEMaren entrenamenduak (3.25). problemarentzat tamaina ezberdineko blokedun Jacobi aurrebaldintzatzaileak erabilita. ......................... 55 3.13 Energia-normaren entrenamendu luzea vs. energiaeta L2-normaren arteko trukea 3.12b. irudiaren bosgarren (s= 5) urratsaren erdian. 56 3.14 Poisson (3.27). probleman lortutako iragarpenak lau urratsetan zehar. u∗soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. ................. 58 3.15 Helmholtz (3.27). probleman lortutako iragarpenak lau urratsetan zehar. u∗soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. ................. 59 3.16 Hasieratik amaierarainoko entrenamenduak DeepFEMean hiru urratsetan zehar (3.27). probleman, 32 tamainako Jacobi blokedun aurrebaldintzatzaileak erabilita. ................... 60 3.17 (3.30). probleman lortutako iragarpenak lau urratsetan zehar. u∗ soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. ...................... 62 3.18 Hasieratik amaierarainoko entrenamenduak DeepFEMean lau urratsetan zehar (3.30). probleman. 48, 48, 96 eta 96 tamainako Jacobi blokedun aurrebaldintzatzaileak erabili dira 1., 2., 3. eta 4. urratsetan, hurrenez hurren. ....................... 63 xviii List of Figures/Irudien zerrenda 3.19 Lau urratsetan lortutako iragarpenak proba datu-basean (3.32). probleman, 0 < α < 200 denean. Galera-funtzioan energia-norma erabili dugu. uFEM(α) elementu finituetako soluzioa da eta uNN(α) DeepFEMaren iragarpena αparametroarentzako. ......... 65 3.20 Galera-funtzioaren eta errorearen energia-normaren eboluzioak lau urratsetan (3.32). probleman, 0 ≤α≤200 denean. ........ 67 3.21 Galera-funtzioaren eta errorearen energia-normaren eboluzioak lau urratsetan (3.32). probleman, −50 < α < −30 denean. ...... 68 3.22 Lau urratsetan lortutako iragarpenak proba datu-basean (3.32). probleman, −50 < α < −30 denean. Galera-funtzioan errorearen H1-norma erabili dugu. uFEM(α) elementu finituetako soluzioa da eta uNN(α) DeepFEMaren iragarpena αparametroarentzako. . . . 69 4.1 Implementation sketch of the general flux for the proposed methods at each training iteration. ...................... 83 4.2 Main model architecture for WANs. It consists of two independent NNs, uNN and vNN, combined via the loss function Lmax min ...... 84 4.3 Main model architecture for the DRM. It consists of a NN, uNN, composed with the trial-to-test operator Tand the loss function L.85 4.4 Main model architecture for the D2RM. It consists of two NNs, uNN and τNN, equipped with the loss functions LτNN and Lopt uNN .. . 85 4.5 Trial network predictions and errors in WANs, the DRM, and the D2RM at the end of the training in model problem (4.24) with exact solution u∗=x(x−1). ..................... 87 4.6 Loss evolution during the WANs training for model problem (4.24) with exact solution u∗=x(x−1). .................. 87 4.7 Loss evolution during the DRM training in model problem (4.24) with exact solution u∗=x(x−1). .................. 88 4.8 Loss evolution of the D2RM training for model problem (4.24) with exact solution u∗=x(x−1). ..................... 88 4.9 Trial network predictions and errors for WANs, the DRM, and the D2RM in model problem (4.24) with exact solution u∗ α=xα(x−1) for α∈ {2,5,10}........................... 91 4.10 Trial network predictions and errors for the DRM and the D2RM in model problem (4.24) with exact solution u∗ α=xα(x−1) for α∈ {0.6,0.7,0.8}. The last row is an augmented version of the third one in a reduced neighborhood of zero. ............ 93 4.11 Trial error functions in model problem (4.24) with with solution (4.29) at different stages of the training progress for the DRM and the D2RM. ............................... 95 xix 1. Sarrera∗ 1.1. Motibazioa Deribatu partzialetako ekuazioek (PDE-ek)1balio handia dute gizartearentzat fenomeno biologiko, fisiko edo sozial anitzak modelatzeko duten gaitasun zabalagatik [73,207,70,192]. Ekuazio horietan deribatuak erabiltzeak aukera ematen digu erlazio konplexuak eta haien denboraren eta/edo espazioaren gaineko aldaketaren tasak deskribatzeko. Adibidez, honako fenomeno hauek PDE-en bidez modelizatutako ereduak onartzen dituzte: bero-transferentzia [102,88], eremu elektromagnetikoak [106,221], fluidoen dinamika [23,14,223], populazioen eboluzioa [96,160,38], eta finantza [26,101] edo osasun [193,141,165] arloko iragarpenak. Hala ere, egoera fisiko edo sozial bat deskribatzeko ekuazio bat ezartzea lehen urratsa baino ez da. Arazo mota ezberdinak sortzen dira PDEak erabiltzean, eta horietara hurbiltzeko modua eskura dauden konputazio-baliabideen araberakoa da. Aurreranzko problemetan, PDE bat eta hasierako edo mugaldeko baldintza batzuk betetzen dituen funtzioa zehaztea da interesatzen zaiguna. Horrela, aurreranzko problemen soluzioek modelatutako fenomenoaren buruzko ezagutza eman diezagukete, eta hori erabil daiteke iragarpenak egiteko, normalean soluziotik hartutako post-prozesamenduetako neurketen bidez. Alderantzizko problemetan, soluzioaren neurriak ditugu, baina ekuazioan parametro espezifikoek duten balioa aurkitu nahi izaten da. Zeharkako arazo bat da, emandako neurketekin bat datorren ebazpen bat ematen duten parametroak zehaztu behar baitira [210]. 1.1. adibidean eremu elektromagnetikoetan dauden aurreranzko eta alderantzizko problemen eredu posible bat erakusten da [72]. ∗See Appendix A for the English version/Ikus A. eranskina ingelesezko bertsiorako. 1Ingelesetik, Partial Differential Equation (PDE). 3 1. Sarrera 1.1. adibidea (Eremu elektromagnetikoak).Maxwellen ekuazioek lau lege fisiko ezagunen arabera modelatzen dute eremu elektromagnetikoa [72,106]:            ∇×E=−jωµH−M,Faradayren Legea, ∇×H= (σ+jωε)E+J,Amper`eren Legea, ∇·(εE) = ρf,Gaussen Elektrizitatearen Legea, ∇·(µH)=0,Gaussen Magnetismoaren Legea. (1.1) Hemen, Eeta Helektrizitatearen eta magnetismoaren eremuak dira frekuentzien domeinuan, hurrenez hurren; Jeta Miturri terminoak dira; σingurunearen eroankortasun elektrikoa da; εpermitibitate elektrikoa da; µiragazkortasun magnetikoa da; ρfkarga elektrikoaren dentsitatea da; junitate irudikaria da; eta ωfrekuentzia angularra da. Orduan, eremu elektromagnetikoen aurreranzko eta alderantzizko eredu posible bat {σ,µ,ε}parametroen eta {Z}neurketen mende osa daiteke, Zinpedantzia tentsorea izanik eta honako era honetan definituta egonik: E=ZH. Xehetasun gehiagoren bila, ikus, adibidez, [9]. 1.1. irudian eredu honen laburpen grafikoa erakutsi da. Parametroak {σ,µ,ε} Soluzioa {E,H} Neurketak {Z} Ebatzi PDEa Post-prozesaketa Aurreranzko problema Alderantzizko problema 1.1. irudia: Aurreranzko eta alderantzizko problemen diagrama eremu elektromagnetiko batean. Aurreranzko problemak ebazteko, orokorrean, alderantzizko problemak ebazteko baino baliabide gutxiago behar dira. Behin parametroak aukeratuta, aurreranzko problema bat PDE simulazio bakar bat burutzean datza. Aitzitik, alderantzizko problema batek, aurreranzko-problemen urratsez urratseko ebazpena eskatzen du, urrats bakoitzean parametroen aukeraketa doituz aurreko aurreranzko simulazioan lortu den emaitzaren arabera [210,103]. Bi arazo motak lotuz, problema parametrikoak ditugu, non helburua soluzio edo neurketen portaera parametroen funtzio gisa aztertzea baita. Soluzioa parametro 4 1.2. Zenbakizko metodo tradizionalak baten edo batzuen araberakoa izan daiteke, eta parametro horien aldaketek sistemaren portaerari nola eragiten dioten azter dezakegu. Alternatiboki, ikuspegi parametrikoa aurreranzko problema gisa har daiteke, non ekuazioaren parametroek soluzioaren aldagaien rola hartzen baitute. Tipikoki, aurreranzko problema parametrikoa ez da lineala sartutako aldagai berriekiko; eta horrek problema ebazteko gaitasuna izugarri mailakatzen eta eragozten du. Tesi honetan, PDE lineal parametrikoak eta ez parametrikoak ebazteko zereginari soilik helduko diogu, non aplikazio eta lan-esparru askotan kritikoa baita. Gainera, PDEak esango diegu bai dimentsio bakarreko ekuazio diferentzialei — oro har, ekuazio diferentzial arruntak (ODEak)2deituak— bai eta dimentsio handiagokoei ere. 1.2. Zenbakizko metodo tradizionalak PDEak ebazteko hainbat metodo daude, oro har, bi kategoriatan sailka daitezkeenak: metodo analitikoak eta zenbakizkoak. Metodo analitikoek PDEaren soluzio zehatz bat aurkitzeko teknika matematikoak erabiltzea eskatzen dute, hala nola, aldagaien bereizketa, transformatu integralak edo analisi konplexua [205,87,213]. Hala ere, metodo horiek, askotan, PDE eta geometria nahiko sinpleetara mugatzen dira. Gainera, PDEa analitikoki ebatz daitekeenean ere, gerta daiteke emaitzaren adierazpena konplikatua izatea (adibidez, potentzia infinitu sorta baten moduan adieraztea). Horrelakoetan, zaila izan daiteke hori interpretatzea eta kasu praktikoetan aplikatzea. Bestalde, zenbakizko metodoak PDEaren soluzioa hurbiltzen dute algoritmo konputazionalen bidez. Oro har, metodo hauek erraz interpreta edo berreskura daitezkeen hurbilketak proposatzen dituzte (adibidez, funtzio sinpleen konbinazio lineal finituak erabilita). Jarraian, labur aztertuko ditugu PDEak ebazteko erabiltzen diren hiru zenbakizko metodo ezagunen alderdi nagusiak. Diferentzia finituen metodoak (FDMak)3PDEaren domeinu espaziala eta/edo denbora tartea azpidomeinuen kopuru finitu batean banatzen du [203,178,130]. Orduan, puntu diskretu horietako deribatuen ebaluazioak diferentzia finituak eta hurbileko puntuen balioak dituzten ekuazio aljebraikoak ebatziz hurbiltzen dira. FDMak PDE bat ekuazio linealen sistema batean bihurtzen du, aljebra matrizialeko tekniken bidez ebatz daitekeena. Kontzeptualki sinplea izan arren, desafiatzailea gerta daiteke arazo konplexuetarako edo domeinu irregularretarako diseinatzen denean. Zehazki, dimentsioen maldizio deritzona jasaten du [29,30], non inplikatutako matrizearen tamaina esponentzialki hazten baita problemaren dimentsioarekin. 2Ingelesetik, Ordinary Differential Equation (ODE). 3Ingelesetik, Finite Difference Method (FDM). 5 1. Sarrera Era berean, elementu finituen metodoak (FEMak)4[100,39,138,184] diskretizazio-sareetan oinarritutako eskema bat proposatzen du, non soluzio hurbildua aurrez ezarritako funtzio batzuen konbinazio lineal finitu baten forman adierazita baitago (normalean, zatikako polinomioak) euskarri lokalekin. Euskarria berariaz diseinatzen da domeinuko eskualde espezifikoetan, elementu deiturikoetan, emaitza bezala lortzen den matrizea sakabanatua izateko. Bereziki, FEMa egokia da PDEak formulazio bariazionaletan daudenean, baina FDMaren antzeko zailtasunak ditu. Metodo espektralak [81,36,46] zenbakizko metodo mota bat dira, eta euskarri globala duten oinarriko funtzio ortogonalen konbinazio linealen bidez hurbiltzen dituzte funtzioak, hala nola, Chebyshev polinomioak edo Fourier serieak. Metodo horiek oso zehaztasun handia lor dezakete eta bereziki erabilgarriak dira portaera oszilakor handiko soluzio leunetarako, baina konputazio-ikuspegitik bideraezinak izan daitezke geometria konplexuetarako. Funtsean, aipatutako zenbakizko metodoek oinarri finitu-dimentsional bat aurrez ezartzea eta dagokion konbinazio linealaren koefizienteen bidez hurbilpena parametrizatzea dute lan-eremu. Diskretizazio-sare bat erabiltzen denean (FDMetan eta FEMetan), horrek behar bezain fina izan behar du soluzio zehatzaren portaera konplexuak behar bezala antzemateko, baina gehiegizko finketaren eragozpen konputazionalean erori gabe. Metodo espektraletan, erronkarik handiena da oinarri ortogonaleko funtzioak domeinu arbitrarioetarako hautatzea. 1.3. Neurona-sareak hurbiltzaile unibertsal gisa Neurona-sare artifizialek (ANNek)5, edo, laburtzeko, neurona-sareek (NNek)6, funtzioak hurbiltzeko alternatiba bat eskaintzen dute, aurretik azaldutako ohiko konbinazio linealen ikuspegiarekiko ezberdina dena. NNen formarik oinarrizkoena aurrerantz elikatutako neurona-sareak (FFNNak)7 dira [190,105,12,77,195], garuneko NN biologikoen egituran eta funtzionamenduan inspiratzen direnak, eta elkarren artean konektatutako neuronen bidez osaturik daudenak. Neuronak geruzetan antolatzen dira, eta informazioa norabide batetik bestera doa, geruza batetik hurrengora. Neurona bakoitzak aurreko geruzako neuronen informazioa jasotzen du eta batura haztatu bat kalkulatzen du. Batura hori aurreko geruzako neurona guztietatik abiatuta kalkulatzen denean, FFNNa guztiz konektatuta dagoela esango dugu. Batuketa hori aktibazio ez-linealeko funtzio batetik pasarazten da, gero emaitza hurrengo geruzako neu4Ingelesetik, Finite Element Method (FEM). 5Ingelesetik, Artificial Neural Network (ANN). 6Ingelesetik, Neural Network (NN). 7Ingelesetik, Feed-Forward Neural Network (FFNN). 6 1.3. Neurona-sareak hurbiltzaile unibertsal gisa ronetara transmititzeko. Guztiz konektatutako FFNNen deskribapen formala 2. kapitulura atzeratuko dugu. Guztiz konektatutako FFNNak hurbiltzaile unibertsalak dira. Era arinean esanda, emandako ia edozein funtzio nahi den adina hurbiltzeko gaitasuna dute. Hurbiltasun unibertsalaren teoremaren (UATren)8lehen bertsioetako bat zabalera arbitrarioko kasua da. Bertan frogatzen da ezkutuko geruza bakarreko FFNNak funtzio jarraituen hurbiltzaile unibertsalak direla, behar bezain zabalak baldin badira. Emaitza hori paraleloki eta independenteki frogatu zen 1989an, “sigmoid” aktibazio funtzioetarako [54] eta aktibazio-funtzio ez-konstante, mugatu eta monotonoki gorakorretarako [98]. Azken lan horren egile nagusiak bi urte geroago frogatu zuen ez dela aktibazio-funtzioaren hautaketa espezifikoa FFNNei hurbiltzaile unibertsal izateko ahalmena ematen diena, baizik eta FFNNaren arkitektura bera [97]. Horrekin lotuta, geroago erakutsi zen baliokideak direla hurbilketa unibertsalaren ezaugarria eta aktibazio ez-polinomikoaren funtzioa izatea [129,175]. UATrako zabalera arbitrarioaren kasuaren bertsio “duala” zabalera mugatua finkatzean eta sakontasun arbitrarioa ezartzean datza. Zenbait egilek formulazio hori aztertu dute 2000ko hamarkadaren hasieratik [83,118], ReLU aktibaziofuntzioetan arreta berezia jarriz [230,140,90,89]. Era berean, zabalera eta sakonera mugatuko emaitzak aurkitzen ditugu. Adibidez, [146]-k baieztatzen du FFNN baterako existitzen dela bi sakonerako eta zabalera mugatuko hurbiltasun unibertsalaren ezaugarriak dituen sigmoidal moduko aktibazio funtzio bat. Gainera, [85]-ek frogatu zuen geruza ezkutu bakarreko eta zabalera mugatuko FFNNek aldagai-bakarreko funtzioetarako hurbiltasun unibertsalak izaten jarraitzen dutela, oro har, aldagai-anitzeko funtzioetarako ez dena betetzen. Azken hiru hamarkadetan, UATren hainbat hedapen garatu dira non ez-jarraituak diren aktibazio funtzioak [129], domeinu ez-trinkoak [118], eta arkitektura eta topologia alternatiboak landu baitituzte [118,136]. Zehazki, [22]-k UATren bertsio kuantitatibo bat ezartzen du FFNNen dimentsionaltasunaren maldizioa gainditzeko gaitasuna aztertuz, zatika linealak diren funtzioen espazioekin alderatuta. Gai horri buruzko eztabaida zabalagoak aztertzeko, [176,19,24,226] lanak eta horietan jasotako erreferentziak aipatzen ditugu. Azken urteotan, hainbat ekimen eta ikerketek ebidentzia edo froga enpirikoen bidez erakutsi dute arkitektura mota alternatiboek edo haien arteko konbinazio egokiek FFNNen hurbiltzeko eta kalkulatzeko gaitasuna gaindi dezaketela. Besteak beste, arkitektura konboluzionalak, errepikariak eta hondarrarenak (ingelesezko, convolutional,recurrent eta residual terminoetatik eta CNN, RNN eta ResNet laburduren bidez adierazitako NNak, hurrenez hurren) dira nabarmenenak [52,59,132,206,77,93,222,4,67,220]. 8Ingelesetik, Universal Approximation Theorem (UAT). 7 1. Sarrera Tesian zehar, guztiz konektatutako FFNN arkitekturetara mugatuko gara. Beraz, NN eta FFNN terminoak bereizi gabe irakur daitezke hemendik aurrera, kontrakorik esaten ez bada. Gainera, nabarmentzekoa da erabilitako arkitekturaren hurbilpen gaitasunen azterketa tesi honen irismenetik kanpo geratzen dela. Bestela esanda, (era inozoan) onartuko dugu esperimentazioan zehar hautatutako arkitekturek “behar adina gaitasun dutela hurbiltzeko”. 1.4. Berrikuspen literarioa Azken urteotan, PDEak NNen bidez ebazteko lan ugari sortu dira. 1.4.1. NNak erabiliz PDEak ebazteko lehenengo lanak 1.3. ataleko terminoei jarraituz NNen bidez PDEak ebazteko lehenengo lan dokumentatua 1994koa9da [63]. Aurreko lan batzuek PDEak ebazteko NNak erabili arren (ikus, adibidez, [127,216,152]), [63] proposamena izan zen PDE baten soluzio baten irudikapena ustiatu zuen lehenengoa, hurbilketa unibertsalaren propietatean oinarritutako NN baten bidez. Lan honek FFNN arkitektura proposatzen du, aktibazio sigmoidaleko funtzioekin, mutur-baldintzak dituen PDE baten soluzioa hurbiltzeko —normalean, mutur-balioen problema (BVPa)10 deitzen dena—. Espezifikoki, BVPa honela deskribatuta dago: ®Au =f, Ω-n, Du =g, ∂Ω-n,(1.2) non Ω domeinu ireki eta mugatua baita, ∂Ω domeinuaren muga da, femandako iturri-funtzio bat da, eta Aeta Doperadore (diferentzial/mugadun) jakin batzuk dira. Egileek honako hau adierazi zuten hitzez hitz (itzulita)11: “(1.2). ekuazioaren soluzio sendo bat Uespazioko funtzioen elementu bat da, PDEa eta mutur-baldintzak betetzen dituena. [...] Hala ere, Uespazio horretan hurbiltzaile ‘unibertsal’ bat baldin badago (hau da, behar adina parametro dituen eta espazioko edozein elementurekiko behar adina hurbil dagoen funtzio bat badago), orduan, hurbiltzaile hori bilatzen ari garen soluziorako aukera ona izan daiteke. Bi funtzioren arteko hurbiltasuna U-ko dentsitatearen kontzeptuaren bidez azaltzen da. Definizio matematiko zehatz baterako, ikus Hornik et al. [98]”. 9Argitaletxeak 1992ko abuztuan jaso zuen, 1993ko abuztuan berrikusi zen eta 1994ko martxoan argitaratu zen lehen bertsioa. 10Ingelesetik, Boundary-Value Problem (BVP). 11Jatorrizko sinboloak eta erreferentziak egokitu egin dira gure aurkezpenerako. [...] ikurrak esan nahi du jatorrizko testuaren zati bat ez dela jaso. 8 1.4. Berrikuspen literarioa Orduan, PDEaren ebazpen-zeregina hurrengo funtzio objektiboaren minimizazio bezala ezar daiteke: F(uNN) = ZΩ∥AuNN −f∥2+Z∂Ω∥DuNN −g∥2,(1.3) non uNN FFNN baita. Jatorriz, ∥·∥ zehaztugabea zen, baina, ziurrenik, 2-norma diskretutzat hartuko zen F-ren minimizazioa L2-minimizazio-eskeman uler dadin ([53]-n iradokitzen den moduan). Aurkezpen horren ondoren, egileek esan zuten12: “Orain, Au eta Du modu itxian kalkula daitezke xterminoaren mende. Ω-n eta ∂Ω-n puntuak hautatzen ditugu eta (1.3) formula aplikatzea proposatzen dugu. Beraz, funtzio objektiboa minimiza dezakegu [NNtik ikas daitezkeen parametroetarako balioak aurkitzeko (hortik dator (1.2)-ren soluzio hurbildua)]”. Egileek zehaztu zuten funtzio objektiboaren minimizazioa “ia-Newton metodo” baten arabera eta “diferentzia finituko gradienteak”erabiliz egin zutela (ikus, adibidez, [61] eta [156] erlazionatutako erreferentzia gisa). Harrezkero, hainbat lan sortu ziren NNak erabiliz PDE ebazpenaren esparruan. Jarraian, 1990eko hamarkadaren amaieratik 2010eko hasierara arteko garrantzitsuenetako batzuk errepasatuko ditugu. 1994an, [154,153,71] lanek FFNN bat erabiltzea proposatu zuten, geruza ezkutu bakarreko arkitekturarekin eta “muga gogorreko”aktibazio funtzioak erabiliz (jatorriz horrela deitu zieten egileek), PDE linealak eta ez-linealak ebazteko. Aktibazio-funtzio hori hautatzeak aukera ematen du NN guztiak zatika-linealak diren funtzio familia gisa ulertzeko, eta, horrela, PDEaren soluzioa termino horietan hurbiltzeko. 1998an, [124] eskema bat proposatu zuen mutur-baldintzak modu “gogor” batean inposatzeko, FFNN (muturrik gabeko) arkitektura filtro batekin konposatzean ([63]-n ez bezala, non mutur-baldintzak “leunki” zehazten baitziren helburu-funtzioaren bitartez). 1999an, [164]-k hainbat emaitza eta diseinu dokumentatu zituen NN bidezko inplementazio ezberdinak erabilita, aldagai anitzeko funtzioen eta horien deribatu partzialen zenbakizko hurbilketarako. Horien artean, PDE aplikazioak daude, fluxu biskoelastikoak aztertzeko inguraketa-elementuen metodo baten parte gisa. 2000. urtean, NNetan oinarritutako eredu bat aurkeztu zen lehen ordenako zenbait PDE ebazteko, kontrol-sistema ez-linealetan interes berezia zuena, eta [125]-ek FFNN motako bi arkitektura ezberdinez osatutako eredua aurkeztu zuen, geometria irregularretan mutur-baldintzak dituzten PDEak ebaztera iristeko. 2001etik 2010eko hamarkadaren erdialdera, arestian aipatutakoen jarraipentzat har daitezkeen hainbat lan daude. Zehazki, [124]-k eta [94]-k aurreka12Lehen bezala, jatorrizko sinboloak eta erreferentziak egokituta daude. Kortxete arteko testuzatiak, [testua], jatorrizko testua aldatu edo ordezkatu dugula adierazten du ulergarriagoa izateko asmoarekin. 9 1. Sarrera ri esanguratsuak ezarri zituzten arkitektura-diseinuetarako eta entrenamendukonfigurazioetarako, hurrenez hurren, honako lan hauetarako: [1]-ek PDEak ebazteko NNak eta algoritmo ebolutiboak eta horien mutur-baldintzak konbinatzen dituen metodo bat proposatu zuen; [202]-k NNak erabili zituen Kuramato-Sivashinsky eta Navier-Stokes ekuazio ez-linealen dinamika aztertzeko; [147]-k metodo hibrido bat aurkeztu zuen ordena handiko PDEak ebazteko; [199]-k Schr¨odingerren ekuazioa landu zuen FFNN erabilita, energia-funtzioaren eskema hobetu batekin eta gainbegiratu gabeko entrenamenduaren bidez; [28]-k minimizazio-tekniken eta kokapen-metodoen konbinazio bat proposatu zuen NNen bidez soluzioetarako hurbilpen itxiak ezartzeko; [215]-ek FFNNak erabili zituen bilakaera gramatikalarekin eta optimizazio lokalarekin batera, ODE eta PDE sistemak ebazteko (ikus [214] bilakaera gramatikalean oinarritutako NNan entrenamenduan izandako aurrekari baterako); eta [151]-k NNan oinarritutako metodo bat aurkeztu zuen, domeinu irregularrak dituzten PDEak ebazteko, mutur-baldintza mistoekin. Gainera, PDE-en bereizmen-esparruen bi luzapen mota nagusi ditugu NNetan oinarrituak: oinarri-erradiala duten NNak13 [142,143,108,144,109,113,145, 76,13,75,47,121]; eta FEMean inspiratutako NNak [31,60,234,183,148,15, 110,120] —batzuetan elementu finituen neurona-sareak deituak—. Aipatutako lan asko berrikusteko, [227] kontsulta daiteke. 1.4.2. NNen bidez PDEak ebazten TensorFlow-tik aurrera TensorFlow (TF)[3,2] kode irekiaren lehen bertsioaren aurretik 2015ean, diferentziaketa automatiko (ADa)14 eraginkorrik ez egoteak erronkak planteatzen zituen NNetan oinarritutako ikerketan. Gaur egun, ADa funtsezkoa da NNen gradienteak kalkulatzeko entrenamenduan zehar [82,27,150]. Hori ez zegoenean, ikertzaileek eta profesionalek eskuz diseinatu behar izaten zuten gradienteen kalkulu-eskema, eta horrek, normalean, denbora asko eramaten zuen eta akatsak egiteko joera areagotu egiten zen, batez ere ereduak konplexuak ziren heinean. Ondorioz, berrikuntza eta esperimentazioa oztopatu egiten ziren. Ideia berriak eta eredu mota desberdinak azkar probatzeko ezintasunak aurrerapena moteltzen zuen, ikerlari berriei eremu horretan sartzea zailduz. 2015ean TF [3,2] eta 2017an PyTorch [171,172] sartu zirenean, ADa eskuragarriago bihurtu zen ikertzaileentzat, eskuzko konputazioaren diseinuetatik askatuz, eta garapen eta berrikuntza azkarragoak bultzatuz. NNetara bideratutako plataforma horien agerpenari esker, zeinak ADetan oinarrituta baitzeuden eta erabilerrazak baitziren, 2010eko hamarkadaren amaieratik aurrera PDE-en ebazpenaz arduratutako komunitateetan NNetan oinarritutako bideak ikertzeko 13Ikusi [170] oinarri-erradial NNen UAT aurrekari baterako. 14Ingelsetik, Automatic Differentiation (AD). 10 1.4. Berrikuspen literarioa interesa piztu zen. Zehazki, fisikatik informatutako neurona-sareak (PINNak)15 PDE-en ebazpenaz arduratutako komunitatean NNak erabiltzearen aurrekari garrantzitsuenetakotzat har daitezke, bereziki azken bost urteetan. PINNen lehen bertsioa bi zatitan argitaratu zen 2017an [180,181]. Han, egileek datuetan oinarritutako alternatibak proposatu zituzten PDEak ebazteko ([180]) edo aurkitzeko ([181]). “Ebaztea” eta “aurkitzea” bata bestearengandik bereizteko, pentsa dezagun koefiziente (parametro) jakin batzuen mende dagoen operadore diferentzial bat duen PDE batean. Alde batetik, koefiziente guztien balioak ezagutzen baditugu, PDEan inplikatutako lege diferentzialerako sarbide osoa eskuragarri izango dugu. Beraz, lehen aipatutako “lehen laneko” [63]16 terminoetara murritz gaitezke PDEaren soluzioa hurbiltzeko (gogoratu 1.4.1. ataleko lehenengo zatia). Bestalde, koefiziente batzuk ezezagunak badira, baina PDEaren soluzioaren lagin batzuk izan baditzakegu espazio-denborako hainbat puntutan, ikaskuntza-eskema gainbegiratu bat ezar dezakegu, non NN bat (PDEaren soluzioa irudikatzen duena) eta parametro entrenagarri gehigarri batzuk (PDEaren koefiziente ezezagunak adierazten dituztenak) PDEa eta etiketatutako datuak asetzeko murrizketarekin entrenatzen baitira. PINN terminologiarekin jarraituz, formulazio bariazionalak jorratzen dituzten lanak aurkitzen ditugu, PINN (hp-)bariazionalak ((hp-)VPINNak) deituak [115,116,188]; kontserbazio-legekoak, PINN kontserbatzaileak (cPINNak) deituak [107]; zatikizko ordena duten PDE-entzako, PINN zatikidunak (fPINNak) deituak [168]; eta datu zaratatsuetarako proposamen bayesiarrak, PINN bayesiarrak (B-PINNak) deituak [229]. PINNek mundu errealeko fenomenoetan dituzten aplikazioen artean, honako hauek nabarmentzen ditugu: [149]-k abiadura handiko fluxuetarako, [157,99]-k potentzia-sistemetarako, [200]-ek ultrasoinuen bidezko kuantifikazio ez-suntsitzailerako, [43]-k beroa transferitzeko eta [42,111]- k fluidoen mekanikarako. [53,126,35]-etara jotzen dugu PINNetan lotutako berrikuspen bibliografiko zehatzagoetarako. Bestela, PDE-en NNetan oinarritutako bereizmenaren helburuari jarraituz, deep izena hartzen duten lanak ere proposatu dira. Batzuk aipatzekotan: [68,69]- kDeep Ritz metodoa (DRMa) aurkezten dute, PDE-en bariazional-formulazioan inplikatutako Ritz-en energia-funtzioa minimizatzen duena ([139,65,135,217] lanak aipatzen ditugu analisi eta hedapen gehigarrietarako); [201]-ek Deep Galerkin 15Ingelesez, Physics-Informed Neural Networks (PINNs). 16Bitxia bada ere, ez [180,181]-ek ezta antzeko egiletza-berrikuspenek ere, esaterako [182, 114], ez dute aipatzen [63]-ren lana. Berrikuspen orokorragoak kontsultatu behar ditugu (egiletzarekiko independenteak) erlazio hori aurkitzeko. Adibidez, [53]-k baieztatzen du (itzulita): “Ikaskuntza automatikoko algoritmo batean aurretiazko ezagutza txertatzearen kontzeptua ez da guztiz berria. Izan ere, Dissanayake eta Phan-Thien-ren [63] lana lehen PINNetakotzat har daiteke”. 11 2. Feed-Forward Neural Networks x=y0y1y2y3uNN Input Hidden Output x(1) x(2) x(3) y(1) 1 y(2) 1 y(3) 1 y(4) 1 y(5) 1 y(1) 2 y(2) 2 y(3) 2 y(4) 2 y(5) 2 y(1) 3 y(2) 3 y(3) 3 y(4) 3 y(5) 3 u(1) NN u(2) NN Figure 2.1: Graph of a fully-connected FFNN with a three-dimensional input, two-dimensional output, depth three, and hidden layers of width five. the parameterization domain, we redefine the FFNN as follows: uNN :X×Θ−→ Y, x , θ 7−→ uNN(x;θ),(2.3a) where θ:= {W1,b1,W2,b2,...,WK,bK,W}(2.3b) stands for the set of learnable parameters with domain Θ := Rn1×n0×Rn1×Rn2×n1×Rn2×···×RnK+1×nK×RnK+1 ×RnK+1 .(2.3c) The functions arising from this parameterization, also known as realizations, belong to a space that, in general, is not a vector space but a manifold [174]. Following the above notation, we denote it by UNN := {uNN(·;θ) : X−→ Y}θ∈Θ,(2.4) and by ΦNN : Θ −→ UNN the realization map that relates each configuration of parameters with the corresponding realization for the given NN architecture, i.e., ΦNN(θ) := uNN(·;θ) : X−→ Y, θ ∈Θ.(2.5) We will drop the boldface notation when the NN is specified as a scalar-valued function, i.e., we will write uNN instead of uNN whenever nK+1 = 1. The set of realizations produced by a given FFNN defined as in (2.3) could be interpreted as a family of finite-dimensional vector spaces. See Appendix 2.A for details. 18 2.2. Continuum, parameterized, and discretized setups −8−4 4 8 −1 1 x φ(x) (a) φ(x) = tanh(x) = ex−e−x ex+e−x −8−4 4 8 0.5 1 x φ(x) (b) φ(x) = σ(x) = ex ex+1 −8−4 4 8 2 4 6 8 x φ(x) (c) φ(x) = ReLU(x) = max{0, x} −8−4 4 8 2 4 6 8 x φ(x) (d) φ(x) = softplus(x) = ln(1 + ex) Figure 2.2: Some usual activation functions appearance: (a) sigmoid, (b) hyperbolic tangent, (c) rectified linear unit, and (d) softplus. 2.2. Continuum, parameterized, and discretized setups Although NNs are often venerated for their tremendous approximation capabilities, they also possess some noteworthy topological undesirable properties. For example, given a NN architecture, the set of realizations is generally non-closed and non-convex [174,40]. This typically translates into possibly being unable to establish an “optimal” realization (that is better than any other) or that multiple ones may exist. This is a substantial departure from the usual approach where the optimal approximation candidate is typically well-defined as the projection onto a pre-established finite-dimensional vector subspace. Aware of these ill-posed omens, we formally present the approximation task via NNs as follows. Let u∗:X−→ Ybe a function that belongs to a certain suitable space of functions U, e.g., an infinite-dimensional vector space, and let F:U−→ Rbe an objective function that characterizes u∗as its well-defined minimizer, i.e., u∗= arg min u∈UF(u).(2.6) 19 2. Feed-Forward Neural Networks Let uNN :X×Θ−→ Ybe a NN architecture according to definition (2.3), and assume that the corresponding set of realizations UNN is embedded in U. Then, we define an optimal approximation of u∗via uNN as u∗ NN ∈Usuch that F(u∗ NN) = inf uNN∈UNN F(uNN).(2.7) Although u∗ NN may not belong to UNN, by assuming that both ΦNN and Fare continuous on Θ and U, respectively, we obtain that both UNN and F(UNN) are (path-)connected. Then, there always exists a sequence of elements in UNN that converges to the accumulation point u∗ NN. In other words, u∗ NN can be approximated with arbitrary precision by means of UNN, as well as F(u∗ NN) by means of F(UNN). Figure 2.3 illustrates this setup. U UNN u∗ u∗ NN u∗ NN Θ ΦNN F F(UNN)F(u∗)F(u∗ NN) R Figure 2.3: Sketch of the parameterization of a NN and its objective-function minimization formulation. So far, we presented the parameterization in NNs as an embedding of a finitedimensional Euclidean space Θ into Uvia the realization mapping ΦNN. Some authors might refer to this process as a discretization due to the utilization of a finite number of parameters to represent the approximating functions. However, to avoid confusion, we will avoid the latter term to distinguish it from the one concerning the input domain Xthat we discuss below. Assume that the objective function possesses the following integral form: F(u) := ZX I(u)(x)dx, u ∈U,(2.8) 20 2.2. Continuum, parameterized, and discretized setups where I(u) stands for the integrand that consists of transformations of u. Then, the natural restriction of Fto UNN is by means of the parameterization ΦNN, i.e., FNN(ΦNN(θ)) = ZX I(ΦNN(θ))(x)dx, θ ∈Θ.(2.9) To address the challenges associated with the exact integration of NN, we consider L: Θ ×XN−→ R(2.10a) defined by L(θ;x1, x2, . . . , xN) = N X i=1 ωiI(ΦNN(θ))(xi),(2.10b) where xi∈Xand ωi>0 are conveniently selected integration points and weights, respectively, so as to resemble a suitable quadrature rule to approximate F◦ΦNN, i.e., F(ΦNN(θ)) ≈ L(θ;x1, x2, . . . , xN), θ ∈Θ, xi∈X, 1≤i≤N. (2.10c) Thus, we obtain a parameterized and discretized representative Lof the objective function F|UNN intended to be computable. Moreover, we may endow a stochastic nature to (2.10) by sampling the integration points {xi}N i=1 ⊂Xaccording to a random distribution every time we estimate the value of the integral. Example 2.1 shows a widely used stochastic integral approximation approach when using NNs. Example 2.1 (Monte Carlo integration).For Vol(X) := RX1dx < ∞and {xi}N i=1 ⊂Xfollowing a random uniform distribution, Monte Carlo (MC) integration consists of selecting ωi= Vol(X)/N in (2.10) for all 1 ≤i≤N, i.e., F(ΦNN(θ)) = ZX I(ΦNN(θ))(x)dx, θ ∈Θ,(2.11a) ≈Vol(X) N N X i=1 I(ΦNN(θ))(xi) = L(θ;{xi}N i=1), θ ∈Θ.(2.11b) Henceforth, we will refer to the parameterized, discretized, and possibly stochastic rendition of the objective function as the loss function. In this way, we distinguish the continuum-level functional that characterizes our minimization formulation, F, from its “computationally feasible” counterpart, L. In data science, it is common to directly encounter the minimization formulation in a discretized and stochastic form. Example 2.2 presents a so-typical supervised learning approach with NNs. 21 2. Feed-Forward Neural Networks Example 2.2 (Supervised learning with NNs).Let D={xi, yi}i≥1⊂X×Ybe an available labeled large database, and let uNN be a given NN model. Then, the approximation task is typically presented as the minimization of L(θ;{xi, yi}N i=1) = 1 N N X i=1 ∥uNN(xi;θ)−yi∥,(2.12) where ∥ · ∥ is a pre-established discrete norm and {xi, yi}N i=1 is a stochastically selected subsample from D, usually called batch, such that N << |D|<∞. This loss-function-level presentation could be viewed as a Monte Carlo approximation of the continuum-level objective function F(u) = 1 Vol(X)ZX∥uNN(x;θ)−u∗(x)∥dx, (2.13) where u∗stands for the continuum-level function from which we extracted the labeled database, i.e., yi=u∗(xi). Overfitting [92,51,77] is a common issue during training, which makes the NN learn to perform exceptionally well on the training data but struggles to generalize to new (unseen) data. In our integration context, this usually means that the NN focuses on satisfying the quadrature requirement only at the integration points (the training data) but possibly misbehaves elsewhere. In [186], the authors show situations where integrating the NN with few fixed points leads to greedy behavior during training, producing extremely undesirable solutions. To overcome this, MC integration is typically utilized due to its stochastic nature. Unfortunately, MC integration demands an immense sample size to avoid making large integration errors (recall Section A.5.3/1.5.3). 2.3. Gradient-based training To carry out the minimization, we adopt first-order gradient-based methods. Here, the term “first order” specifies that we only calculate gradients with respect to the learnable parameters (and not higher-order derivatives). The cost of computing the gradient of a scalar-valued function is known to be bounded by the cost of evaluating the function itself times a constant when using a usual backpropagation algorithm via Automatic Differentiation (AD) [82, 27,150]. Indeed, such a constant is known to be smaller than or equal to three under specific theoretical conditions (Baur-Strassen Theorem [25]). Hence, the computational cost of calculating the gradient of the objective/loss function is comparable to simply evaluating it. 22 2.3. Gradient-based training However, this fact cannot be extrapolated to vector-valued functions. There, the cost of computing its gradient behaves as the cost of computing the componentwise gradients, i.e., it is comparable to the cost of evaluating the vector function multiplied by the number of its (scalar-valued) output components. Then, second-order minimization methods for NNs have a cost proportional to evaluating the loss/objective function times the number of learnable parameters. Since NNs typically follow an over-parameterized design to exploit their approximation capacities, second or higher-order optimization methods become extremely expensive. In what follows, we present the gradient-descent method that serves as the foundation of the majority of currently employed optimizers when training NNs. We first introduce this method at a continuum level to later deduce its parameterized, discretized, and stochastic counterpart. 2.3.1. Continuum-level Gradient-Descent method Let Ube a Hilbert space and, in addition to the well-posedness, assume that F:U−→ Ris convex, i.e., for u, w ∈Uand 0 ≤λ≤1, it satisfies F(λu + (1 −λ)w)≤λF(u) + (1 −λ)F(w),(2.14) and sufficiently differentiable. Then, the gradient1operator ∇F is monotonically increasing (Kachurovskii’s Theorem [112]), and thus provides an ideal framework to carry out a continuum-level gradient-descent minimization [204] defined as follows: ut=ut−1−λ∇F(ut−1),(2.15) where λ > 0 is the learning rate and ut∈Uis the function at the tth iteration. Intuitively, the gradient-descent method begins selecting an initial candidate u(0) ∈Uand, iteratively, follows steps against the gradient, which by an appropriate (adaptive) control of the size of the learning rate, it is well known that we can approach u∗as much as desired [32,204]. Figure 2.4 illustrates a graphical performance of the gradient-descent method on a convex scenario, and Example 2.3 briefly analyzes a simple continuum-level case of study. Example 2.3. Let ∥·∥Udenote the norm induced from the inner product of the Hilbert space U. Consider the objective function defined by F(u) := ∥u∥U. Then, Fis convex on Uand Fr´echet differentiable on U\ {0}with gradient 1In Hilbert spaces, we refer to the operator defined via the Riesz representatives of the Fr´echet derivative at each point, i.e., at a Fr´echet differentiable point u∈U,∇F(u) denotes the unique element in Uthat satisfies DuF(w)=(w, ∇F(u))Ufor all w∈U, where DuFdenotes the usual Fr´echet derivative at u∈Uand (·,·)Ustands for the inner product in U. 23 2. Feed-Forward Neural Networks u∈U F(u)∈R Figure 2.4: Gradient-descent performance sketch on a convex scenario. ∇F(u) = u ∥u∥U. In consequence, given an initial candidate u0∈U\ {0}, the continuum-level gradient-descent iterative method described in (2.15) becomes ut=ut−1−λut−1 ∥ut−1∥U =Å1−λ ∥ut−1∥Uãut−1.(2.16) Selecting 0 < λ =λ(t)≤ ∥ut−1∥U, the iterative method converges to 0 ∈U, which is the global minimum of Fon U. Indeed, if at any iteration we select λ(t) = ∥ut−1∥U, then we would reach the minimum at that moment. 2.3.2. Parameterized-level Gradient-Descent method Because the continuum-level gradient-descent method is computationally intractable, to turn it into a practical method, we resort to the composition of the objective function with the available realization mapping, F ◦ΦNN : Θ −→ R,(2.17) and apply the corresponding gradient operator ∇(F ◦ΦNN). Because we are now in a finite-dimensional Euclidean domain, we identify the gradient with the vector of partial derivatives given by ∂(F ◦ΦNN) ∂θ := ï∂(F ◦ΦNN) ∂θ(i)òdim(Θ) i=1 ,(2.18) 24 2.3. Gradient-based training where θ(i)denotes the ith learnable parameter of θ. Consequently, the gradientdescent iterative method transletes into the domain of the set of learnable parameters as follows: θt=θt−1−λ∂(F ◦ΦNN) ∂θ (θt−1),(2.19) where θt∈Θ stands for the set of learnable parameters at the tth iteration. While (2.19) can be seen as a parameterized-level version of (2.15), it is essential to note that some of the ideal conditions observed at the continuum level in the gradient-descent method, such as the convexity, do not necessarily satisfy at the parameterized level when dealing with NNs. Indeed, (2.19) merely aims to converge to parameters in Θ whose partial derivatives tend to vanish, which may lead to undesirable results (e.g., when falling in local minima or saddle points). This supposes another departure from the usual vector-space-based approach where the convexity condition of Fis preserved in linear parameterizations (see Example 2.4). Example 2.4 (Convexity preservation under linear transformations).For simplicity, we first present the case of linear combinations. Let Fbe convex on U and consider ΦNN defined by ΦNN(θ) = dim(Θ) X i=1 θ(i)u(i), θ ={θ(i)}dim(Θ) i=1 ,(2.20) where {u(i)}dim(Θ) i=1 ⊂Uis a pre-established set of functions. Then, UNN is a vector space of dimension at most dim(Θ) and F ◦ΦNN is convex on Θ. Proof. The dimension of UNN is the number of linearly independent functions in {u(i)}dim(Θ) i=1 ⊂U. To prove the convexity of F ◦ΦNN, let 0 ≤λ≤1 and θ, ϑ ∈Θ. By linearity, ΦNN(λθ + (1 −t)ϑ) = λΦNN(θ) + (1 −λ)ΦNN(ϑ). We obtain the desired result by applying Fand the convexity definition. Following a similar reasoning, it is straightforward to check that the convexity condition is also preserved when replacing the linear combination with a general linear transformation. In particular, the restriction of F ◦ΦNN to the learnable parameters of the output layer in a FFNN according to Section 2.1 is convex whenever Fis convex. 2.3.3. Stochastic Gradient-Descent method To avoid the computational challenges associated with the integral form of F, we resort to the loss-function approximation introduced in (2.10). Then, the 25 2. Feed-Forward Neural Networks so-called Stochastic Gradient Descent (SGD) optimizer arises as θt=θt−1−η∂L ∂θ (θt−1;xt),(2.21) where xt∈XNdenotes the stochastically sampled set of integration points at the tth iteration. Note that (2.21) introduces discrete dependencies and uncertainty during training compared to (2.19). Multiple variants of stochastic gradient-based optimizers have been designed and exploited in recent years with the aim of improving many of the aforementioned poor conditions in which NNs are immersed. Here are some of the most popular today: SGD with momentum [179] and/or Nesterov acceleration [162], Adagrad [66], Adadelta [233], RMSprop [95], Adam/AdaMax [119], and Nadam [64]. We refer to [191] for an overview. 2.4. A case of study To illustrate many of the (undesirable) properties pointed out so far, we discuss the following easy-to-analyze case of study. At the continuum level, let X:= (−1,1), Y:= R, and u∗(x) := sgn(x) =      −1,if x < 0, 0,if x= 0, 1,if x > 0. (2.22) Consider U:= L2(X) = {u:X−→ R:RXu2<∞} and let F(u) := ∥u−u∗∥2 L2(X)=ZX (u(x)−u∗(x))2dx, (2.23) whose minimum in Uis zero and is attained only at u∗. At the parameterization level, let Θ := R2and let uNN :X×Θ−→ Rbe defined by uNN(x;θ) := ctanh(ax),(2.24) where θ={a, c}denotes the set of learnable parameters. Then, UNN = ΦNN(Θ) = {ΦNN(a, c) = ctanh(ax) : X−→ Y}a,c∈R,(2.25) 26 2.4. A case of study and routinary calculations yield (F ◦ΦNN)(a, c) = 2(1 + c2)−4clog(cosh(a)) + 2c2tanh(a) a,(2.26) ∂(F ◦ΦNN) ∂a (a, c) = 2c2 log(cosh(a)) + (c−2a) tanh(a)−ac(sech2(a) a2,(2.27) ∂(F ◦ΦNN) ∂c (a, c) = 4 (ac −log(cosh(a)) −ctanh(a)) a,(2.28) for a= 0, (F ◦ ΦNN)(0, c) = 2, ∂(F◦ΦNN) ∂c (0, c) = −2c, and ∂(F◦ΦNN) ∂a (0, c) = 0. Figure 2.5 shows the graph of F ◦ ΦNN, and Figure 2.6 shows the graphs of its partial derivatives. The following statements expose a series of properties of this case of study: 1. The architecture uNN posseses arbitrary precision approximation capacity, which means that u∗ NN =u∗. Proof. Let c= 1. Then, it is straightforward to check that (F ◦ΦNN)(a, 1) converges to zero as a→ ∞, which implies that ΦNN(1, a) converges to u∗ as a→ ∞. 2. The set of realizations UNN is non-closed. Proof. The optimal realization u∗ NN ∈L2(X)\C∞(X) is unattainable because ΦNN(a, c)∈C∞(X) for all a, c ∈R. 3. The realization map ΦNN is non-injective. In particular, the sequence of learnable parameters whose images through ΦNN converge to u∗ NN is nonunique, although the limit realization u∗ NN is unique. Proof. The odd symmetry of tanh(x) implies ΦNN(a, c) = ΦNN(−a, −c) for all (a, c)∈Θ, which shows the non-injectivity. Then, for c= 1, we find distinct sequences {(a, 1)}a≥1and {(−a, −1)}a≥1of learnable parameters with distinct limits but whose images through ΦNN converge to u∗ NN as a→ ∞. The uniqueness of u∗ NN follows from item 1. and the well-posedness of the objective-function minimization. 4. The set of realizations UNN is non-convex. Proof. We need to find 0 <λ<1 and a1, a2, c1, c2∈Rsuch that f(x) = λc1tanh(a1x) + (1 −λ)c2tanh(a2x) is unrealizable via ΦNN, i.e., it is not of the form ΦNN(a, c) = uNN(x;a, c) = ctanh(ax) for any a, c ∈R. Indeed, let λ= 1/2, c1=c2= 2, a1= 2, and a2=−1. Then, f(x) = tanh(2x)− tanh(x), which is neither constant nor strict monotonic on X. Because uNN(x;a, c) is either strict monotonic or constant on X,f(x) is unrealizable via ΦNN. 27 3. Deep elementu finituen metodoa∗ Laburpena. Deribatu partzialetako ekuazio parametriko linealak ebazteko ikasketa sakoneko arkitektura dinamiko bat aurkezten dugu, elementu finituen metodoan oinarrituta dagoena. Arkitekturako neuronen arteko loturek elementu finituetako grafoetakoak antzeratzen dituzte, horietan diskretizazio-sarea fintzen denean. Lan honetan, zenbait galera-funtzio aztertu ditugu, konbergentzia hobetzeko hainbat aurrebaldintzatzaile eta norma erabiliz. Sinpletasunagatik, metodoaren inplementazioa dimentsio bakarreko (1Dko) espazioko eremuan gauzatu dugu, nahiz eta 2D eta 3D problemetara ere heda daitekeen. Egin ditugun zenbakizko esperimentuek, orokorrean, agerian uzten dute metodoak balio duela problema simetriko positibo definituetarako edo problema zehaztugabeetarako, bai ekuazio parametrikoetan baita ez-parametrikoetan ere. Hala ere, kasu batzuetan, ganbiltasun ezak eragotzi egiten du zehaztasun handiko soluzioak lortzea. Argitaratutako bertsiorako, ikus [218]. 3.1. Sarrera Berrikuspen literarioko lan ugaritan lantzen da PDE (parametriko) baten soluzioa hurbilduko duen funtzio jarraitu bat (NNa) aurkitzeko ideia (gogoratu 1.4. atala). Oro har, diseinu horiek NNa eremuko edozein puntutan ebaluatzeko aukera ematen dute, hau da, diskretizazio-sarearen mendekotasunik ez duen egitura daukate. Hala ere, badituzte muga batzuk. Honako bi hauek nabarmenduko ditugu: (a) gehienetan, NN arkitekturek ez daukate azaltzeko gaitasunik [194,21], eta (b) galera-funtzioaren barruan zenbakizko integrazio arauak diseinatzea ez da izaten erraza, [115]-en esaten den bezala, eta xehetasun handiagoz [186]-n. Aipatutako bi zailtasunak gainditzeko, PDE parametrikoen ebazpenerako FEMean oinarritutako ikaskuntza sakoneko metodo bat proposatzen dugu, Deep elementu finituen metodoa (DeepFEM)2edo elementu finituen metodo sakona deitua. Proposatutako NNaren arkitekturak FEMa aplikatu ostean lortzen den ekuazio linealetako sistema parametrikoaren ebazle gisa jardungo du, diskretizazio- ∗See Appendix B for the English version/Ikus B. eranskina ingelesezko bertsiorako. 2Ingelesetik, Deep Finite Element Method (DeepFEM). 35 3. Deep elementu finituen metodoa sarea fintzean elementu finituen konektagarritasun-grafoa imitatuz. NNaren geruza bakoitzak ResNet diseinua du [93] eta diskretizazio-sare zakarragoetako soluzioak diskretizazio-sare finagoetara hedatzen ditu. 3.1. irudiak FEMeko diskretizazio-sareen fintzeen eta NNko geruzen arteko erlazioa erakusten du. Horrela, arkitekturak neurri bateko azalgarritasuna du, eta sarearen iragarpena diskretizaziosare fineko nodoetako balioez osatutako bektorea da. Gainera, ikuspegi diskretu horrek zenbakizko integrazio zehatza ahalbidetzen du entrenamenduan zehar, NNak iragartzen duen soluzioa zatika definitutako polinomioen espazioan dagoelako. FEMaren diskretizazio-sare fintzeak 4. geruza 3. geruza 2. geruza 1. geruza DeepFEMaren geruzak 3.1. irudia: FEMaren diskretizazio-sare fintzeen eta DeepFEMaren geruzen arteko lotura. DeepFEMak, lehendabizi, hasierako arkitektura bat ezartzen du. Arkitektura hori, entrenamenduaren ostean, gai izango da soluzio zakarrak iragartzeko. Gero, prozesu iteratibo eta dinamiko bat jarraituz, arkitekturan geruzak sartuz joango gara, aurretik entrenatutako aldagaiak3mantenduz eta berriak gehituz. Ondoren, eredu berriaren aldagaiak berriro entrenatuko ditugu; eta prozesu bera errepikatuko da nahi dugun zehaztasuneko soluzioa lortu arte. Horrela, proposatutako NNak aukera ematen digu diskretizazio-sarearen konbergentzi-azterketa egiteko. Gure inplementazioa dimentsio bakarreko (1Dko)4problemetara mugatzen da, zatika linealak diren hurbilpenen bidez diskretizazio-sare uniformeen gainean. Hala ere, lana zuzenean heda daiteke dimentsio handiagoko problemetara, diskretizazio-sare handiagoko hurbilketa polinomikoetara, edota egokitutako diskretizaziosareetara. Hala ere, inplementazio landuagoa behar da irizpide geometrikoak kontuan izateko eta diskretizazio-sarearen fintzearen eraginez nodoen zenbaketa 3Kapitulu honetan zehar, 2. kapituluan erabilitako parametro entrenagarrien terminologia aldagai entrenagarrietara aldatuko dugu, PDEaren koefizienteak eta NNaren pisuak eta alborapenak (ingelesetik, bias terms) bereizteko. 4Ingelesetik, one-dimensional (1D). 36 3.2. Eredu problema eta elementu finituen formulazioa egiteko. Lan honetan ez dugu gai horietan sakonduko. Hainbat eredu problematan lortu ditugun zenbakizko emaitzak erakutsiko ditugu, bai parametro konstanteak erabili direnean baita zatika kostanteak erabili direnean ere. Aurkeztutako teknologiaren ekarpen nagusia da NNak problema parametrikoak ebazteko duen gaitasuna agerian uztea. Adibide gisa, lehenengo, parametrikoak ez diren kasuak aurkeztuko ditugu, eta, ondoren, parametrikoak. Gauza bera egingo dugu zenbakizko emaitzak deskribatzean; metodoak ingurune ezparametrikoetan duen ulergarritasunak eta mugak kasu parametrikora heda baitaitezke. Kapitulu honen gainerako zatiak honako hauek dira. 3.2. atalean gure intereseko problema eta hari dagokion formulazio bariazionala azalduko dira. 3.3. atalak DeepFEM arkitektura deskribatuko du eta 3.4. atalak hautatutako galerafuntzioa definituko du. 3.5. atalak TF liburutegia [3,2] erabiltzean aurkitu ditugun inplementazio-mugak eta xehetasunak erakutsiko ditu. Amaitzeko, 3.6. atalean lortutako zenbakizko emaitzak eztabaidatuko ditugu. 3.2. Eredu problema eta elementu finituen formulazioa BVP parametriko jakin bat hartuko dugu, nahiz eta aurkeztuko dugun metodologia FEMa erabiliz ebatz daitezkeen beste problema batzuetarako ere aplikagarria izan. Izan bedi Ω eremu leun bat eta ezar dezagun honako BVP parametrikoa:      −∇·σ∇u+αu =f, Ω-n, u= 0,ΓD-n, −σ∂u ∂n =g, ΓN-n, (3.1) non σ > 0 eta α∈Rparametroak zatika konstanteak diren funtzioak baitira, f iturria, eta gNeumannen datua. ΓDeta ΓNDirichlet eta Neumann mugaldeak dira, hurrenez hurren. nkanporanzko bektore normala da ΓN-ren puntu bakoitzean, eta ∂u/∂n := ∇u·n. Kontuan izan difusio-erreakzio eredu orokor honek Poissonen (α= 0) eta Helmholtzen (α < 0; σ= 1) ekuazioak barne hartzen dituela. Goiko BVPren formulazio bariazional bat honako hau da: ßAurkitu u∗∈H1 0(Ω) izanik, honako hau beteko duena: (σ∇u∗,∇v)Ω+ (αu∗, v)Ω= (f, v)Ω−(g, v)ΓN,∀v∈H1 0(Ω), (3.2) non (u, v)Ω:= ZΩ u·v, (3.3) 37 3. Deep elementu finituen metodoa eta hurrengo espazioak kontuan hartuta: L2(Ω) = {u: Ω −→ R: (u, u)Ω<∞},(3.4a) H1(Ω) = {u∈L2(Ω) : (∇u, ∇u)Ω<∞},(3.4b) H1 0(Ω) = {u∈H1(Ω) : u|ΓD= 0}.(3.4c) Prozesu osoan zehar suposatuko dugu feta gnahiko erregularrak direla, hau da, f∈L2(Ω) eta g∈H1/2(∂Ω). Sobolev espazio mugatuei/zatikariei buruzko eztabaida zehatz baterako, ikus [74]. 1Dko FEM formulazioa erabilita, honako forma hau daukan soluzioa aurkitu nahi dugu: uFEM(x;σ, α) := J X j=0 uFEM,j(σ, α)ψj(x),(3.5) non uFEM,j(σ, α) koefiziente ezezagunak izanik eta ψj(x) zatika linealak diren denda itxuradun oinarriko funtzioak (ikus 3.2. irudia). x0 ψ0(x) x1x2x3x4. . . xJ−1xJ ψ1(x)ψ2(x)ψ3(x)ψJ(x) 3.2. irudia: Denda itxuradun eta zatika linealak diren oinarriko funtzioen, ψjen, euskarria. (3.5) eta (3.2) ekuazioak konbinatuz eta ψjoinarriko funtzioekin testatuz, honako ekuazio linealetako sistema honetara iritsiko gara: Au =f,(3.6) non A=A(σ, α) := [(σψ′ j, ψ′ i)Ω+ (αψj, ψi)Ω]J i,j=0,u:= uFEM(σ, α) := [uFEM,j]J j=0 koefiziente ezezagunen bektorea izanik, eta f:= [(f, ψj)Ω−(g, ψj)ΓN]J j=0 karga bektorea baita. 3.3. Arkitektura dinamikoa Lehenengo, problema ez-parametrikoetarako deskribatuko dugu proposaturiko arkitektura, eta, gero, kasu parametrikora hedatuko dugu. Azkenik, parametroak konstanteak edota zatika konstanteak direneko kasuak aztertuko ditugu. 38 3.3. Arkitektura dinamikoa 3.3. irudiak erakusten du fintze uniformeak egiterakoan nodoak zenbatzeko aukeratu dugun modua. Horren arabera, Ezabaltze-operadorea matrize sakabanatu baten bidez emana dator. Matrize hori 1 edo 1/2 zenbakiez beteta dago, nodo bakoitza diskretizazio-sare zakarragotik finagora hedatzen den moduaren arabera. Kontuan izan, antzeko zabaltze-operadoreak existitzen direla bi dimentsioko (2Dko) 5eta hiru dimentsioko (3Dko)6problemetan, edota ordena handiagoko elementuetan, baita H(div), H(curl) eta L2diskretizazioetan ere [57]. n0 n1 n0 n1 n2 n0 n1 n2 n3 n4 n0 n1 n2 n3 n4 n5 n6 n7 n8 diskretizazio-sare finagoak 3.3. irudia: 1Dko diskretizazio-sare baten fintze uniformeen eta nodoen zenbaketaren progresioa. 3.3.1. Eskema ez-parametrikoa koefiziente konstantedunetarako Izan bitez σand αkonstante errealak. Elementu bakar bateko diskretizaziosarerako, honako bi-geruzako sakonera duen arkitektura hartuko dugu: u(1) NN(σ, α) = W(1) φÅW(1) σ,α ïσ αò+b(1) σ,αã∈R2,(3.7) non φaktibazio funtzioa baita eta output bektorea bi dimentsiokoa izanik, (elementu bakarreko) diskretizazio-sare zakarrean FEMaren soluzioa hurbiltzea helburu duena (s= 1 urratsa), hau da, u(1) FEM(σ, α) = ñu(1) FEM,0(σ, α) u(1) FEM,1(σ, α)ô≈ñu(1) NN,0(σ, α) u(1) NN,1(σ, α)ô=u(1) NN(σ, α).(3.8) 5Ingelesetik, two-dimensional (2D). 6Ingelesetik, three-dimensional (3D). 39 3. Deep elementu finituen metodoa Entrenagarriak diren aldagaien multzoa honako hau da: θ(1) ={W(1) σ,α,b(1) σ,α,W(1)}.(3.9) Gero, diskretizazio-sarea finduko dugu bi elementuko diskretizazio-sarea lortu arte (s= 2 urratsa). Orduan, sarrera-mendekotasuna daukan ResNet [93] motako geruza bat sartuko dugu u(1) NN arkitekturan u(2) NN lortzeko: r(2) NN(σ, α) = W(2) φÅW(2) σ,α ïσ αò+b(2) σ,αã∈R3,(3.10a) u(2) NN(σ, α) = E2 1u(1) NN(σ, α) + r(2) NN(σ, α)∈R3,(3.10b) non E2 1zabaltze-matrizea baita, u(1) NN bektorea diskretizazio-sare finean (bi elementukoan) egon dadin. Orain, aldagai entrenagarrien θ(2) multzoa honako hau da: {W(2) σ,α,b(2) σ,α,W(2)}edo θ(1) ∪{W(2) σ,α,b(2) σ,α,W(2)}(3.11) entrenamendua geruzaz geruza edo hasieratik amaieraraino nola egiten dugun kontuan hartuz, hurrenez hurren. Kontuan izan aldagai berriak r(2) NN-ri bakarrik dagozkiola. Berrentrenamenduaren hasieran r(2) NN = 0 izan dadin aldagai berriak hasieratuko ditugu (adibidez, W(2),b(2) = 0 hasieratuz), eta aurreko urratsean ikasitako aldagaiak, θ(1), mantenduko ditugu s= 2 urratsaren hasieran. u(2) NN berriro entrenatuko dugu (geruzaz geruza edo hasieratik amaierarainoko entrenamenduekin) FEMaren soluzioa hurbil dezan (bi elementuko) diskretizazio-sare finean, hau da, u(2) FEM(σ, α) =    u(2) FEM,0(σ, α) u(2) FEM,1(σ, α) u(2) FEM,2(σ, α)   ≈   u(2) NN,0(σ, α) u(2) NN,1(σ, α) u(2) NN,2(σ, α)   =u(2) NN(σ, α).(3.12) Prozesu hori urratsez urrats errepikatuko dugu, NNaren sakontasuna handituz, diskretizazio-sarearen elementu kopurua soluzio analitikoa zehazki hurbiltzeko gai den arte. 3.4. irudian ikus daiteke arkitektura dinamiko hori. 3.3.2. Eskema parametrikoa koefiziente konstantedunetarako Eskema ez-parametrikoa dimentsio oso txikiko pisuak eta alborapenak ezartzean zetzan (adibidez, W(1),W(1) σ,α ∈R2×2eta W(2) ∈R3×2,W(2) σ,α ∈R2×2), koefizienteen lagin bakar batean entrenatu nahi baikenuen. Problema parametrikoetarako (koefizienteen laginen datu-base batean entrenatzeko), sakonera eta zabalera gehituko dizkiegu arkitektura dinamikoaren atal 40 3.3. Arkitektura dinamikoa u(1) NN σ α 0 1 u(2) NN 0 1 2 1 0.5 0.5 1 u(3) NN n0 n1 n2 n3 n4 0 1 2 3 4 1 0.5 0.5 1 0.5 0.5 1 r(s) NN, s = 2,3 3.4. irudia: Eskema ez-parametrikoaren DeepFEM arkitektura dinamikoaren irudikapena koefiziente konstantedunetarako. Gezi grisek aldagai entrenagarriak adierazten dituzte, eta gezi beltzek, berriz, (entrenaezinak diren) zabaltze-operazioak. entrenagarriei, honela: u(1) NN =FC(1) NN(σ, α),(3.13a) r(s) NN =FC(s) NN(σ, α), s ≥2,(3.13b) u(s) NN =E(s) s−1u(s−1) NN +r(s) NN, s ≥2.(3.13c) Hemen, FC(s) NN-k adierazten du aktibatu gabeko azken geruza batekin guztiz konektatutako FFNN bat (gogoratu 2.1. atala), eta horren irteera-dimentsioa zehaztuko da FEM diskretizazioaren dimentsioekin bat etortzeko eran. Sinpletasunagatik, bloke entrenagarri deituko diogu FC(s) NN ereduei. 3.5. irudiak bloke entrenagarriaren arkitektura erakusten du eta 3.6. irudiak horiek irudikatzen ditu DeepFEM eskema parametrikoaren barruan. 3.3.3. Eskema parametrikoa zatika konstanteak diren koefizienteak dituztenetarako Azkenik, zatika konstanteak diren koefizienteen portaera erantsiko diegu parametroei. Horretarako, pentsatuko dugu parametroek balio konstanteak dituztela hasierako diskretizazio-sareko elementu bakoitzean. Horrela, arkitekturaren aldaketek soilik sarrerako geruzan dute eragina: σeta αbalio errealen ordez, hasierako 41 3. Deep elementu finituen metodoa σ α . . . FC(s) NN 3.5. irudia: s. urratsean bloke entrenagarri baten arkitekturaren irudikapena. Gezi grisek aldagai entrenagarriak adierazten dituzte. u(1) NN σ α 0 1 2 3 4 0 1 2 0 1 FC(1) NN FC(3) NN FC(2) NN u(2) NN 1 0.5 0.5 1 u(3) NN n0 n1 n2 n3 n4 1 0.5 0.5 1 0.5 0.5 1 3.6. irudia: Eskema parametrikoaren DeepFEM arkitektura dinamikoaren irudikapena koefiziente konstantedun problemarako. Gezi grisek aldagai entrenagarriak adierazten dituzte, eta gezi beltzek, berriz, (entrenaezinak diren) zabaltze-operazioak. diskretizazio-sarearen elementu kopuruarekin bat datorren dimentsioko bektoreak erabiliko dira. 3.7. irudiak arkitekturaren grafo dinamiko eta parametrikoa erakusten du, bi elementuko hasierako diskretizazio-sarearentzat. 42 3.4. Galera-funtzioa eta entrenamendua u(1) NN σ1 σ2 α1 α2 FC(1) NN 0 1 2 FC(2) NN 0 1 2 3 4 u(2) NN 1 0.5 0.5 1 0.5 0.5 1 u(3) NN n0 n1 n2 n3 n4 n5 n6 n7 n8 0 1 2 3 4 5 6 7 8 FC(3) NN 3.7. irudia: Eskema parametrikoaren DeepFEM arkitektura dinamikoaren irudikapena zatika konstanteak diren koefizienteak dituztenetarako. Gezi grisek aldagai entrenagarriak adierazten dituzte, eta gezi beltzek, berriz, (entrenaezinak diren) zabaltze-operazioak. 3.4. Galera-funtzioa eta entrenamendua 3.3. atalean urratsez urrats deskribatutako metodologiaren barruan, “(s)” goiindizea ezabatuko dugu sinpletasunagatik, hautatutako aldagai entrenagarriak, galera-funtzioak eta normak urratsaren mende daudela ulertuta. NNak FEMa aplikatu ostean sortzen den ekuazio linealetako sistema parametrikoaren soluzioa hurbil dezan, hau da, uNN ≈uFEM,(3.14) honako galera-funtzioa aukeratuko dugu: L(θ, {σi, αi}N i=1) := 1 N N X i=1 ∥A(σi, αi)uNN(σi, αi;θ)−f∥,(3.15) non uNN(σi, αi;θ) DeepFEM ereduaren iragarpena baita i. laginerako, A(σi, αi) FEM matrizea —gogoratu (3.6)—, A(σi, αi)uNN(σi, αi;θ)−fhondar bektorea, eta ∥·∥ ezarritako norma bektorial bat (adibidez, 2-norma). 3.4.1. Gradientean oinarritutako entrenamendua Galera-funtzioa minimizatzeko, Adam optimizatzailea [119] aplikatuko dugu urrats bakoitzaren hasieran. Exekuzioa geldituko dugu iterazio kopurua gehienezko 43 3. Deep elementu finituen metodoa u∗uFEM uNN 0 0.25 0.5 0.75 1 0 0.5 1 u(x) (a) s= 1 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 0 0.5 1 (b) s= 1 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 0 0.5 1 u(x) (c) s= 2 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 0 0.5 1 (d) s= 2 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 0 0.5 1 u(x) (e) s= 3 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 0 0.5 1 (f) s= 3 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 0 0.5 1 x u(x) (g) s= 4 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 0 0.5 1 x (h) s= 4 urratsa (entrenatu ondoren). 3.8. irudia: (3.25). problemaren lehenengo lau urratsetako DeepFEMeko iragarpenak. u∗soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. 50 3.6. Zenbakizko esperimentuak nez, 2.000 eta 4.000 iterazio ezarri ditugu Adam-entzako eta Adalr-entzako fase bakoitzeko exekuzioetan, hurrenez hurren. Hasieratik amaierarainoko entrenamenduak egiten ditugu NNan. Konbergentzia-desbideraketa nabaritzen dugu urrats-kopurua (eta diskretizazio-sarearen tamaina) handitzean. Nolanahi ere, azken errorea 10−8baino txikiagoa da urrats guztietan zehar (ikus 3.9c. irudia). 3.6.1.1. Hasieratik amaierarainoko entrenamendua vs. geruzaz geruzakoa Eraikuntzarengatik, DeepFEMaren geruza bakoitzak diskretizazio-sarearen elementu bakoitzarekin elkartutako oinarriko funtzioekin lotutako koefizienteak iragartzen ditu. Diskretizazio-sareak zenbat eta finagoak izan, orduan eta lokalagoak dira oinarriko funtzioei lotutako euskarriak. Alde batetik, hasieratik amaierarainoko entrenamenduak koefizienteen hierarkia osoaren entrenagarritasuna ahalbidetzen du. Beste alde batetik, geruzaz geruzako entrenamenduak, bakarrik azken diskretizazio-sarearekin (finena dena) lotuta dauden koefizienteen entrenagarritasuna hornitzen du. 3.10a. irudiak galera-funtzioaren konbergentzia erakusten du, goiko baldintza berberetan, baina geruzaz geruzako entrenamendua egiten denean. Ikusten den moduan, konbergentzia desorekatu egiten da iterazio kopurua handitzen den heinean, hasieratik amaierarainoko entrenamenduarekin konparatuta (gogoratu 3.9b. irudia). Hori gertatzen da galera-funtzioak errorearen gradientea dakarrelako energia-norma erabiltzen denean. 3.10b. irudiak errorefuntzioa irudikatzen du, ia konstantea dena. Horrela, errorearen deribatua ia zero da. Orduan, erronka izaten da errorearen energia-norma minimizatzea, lokalak diren euskarri funtzioei lotutako koefizienteak doituz. Horrek, askotan, errorearen deribatua handitzea dakar. Konbergentzia azkartu egiten da euskarri globaleko eta lokaleko oinarriko funtzioak erabiltzen direnean, multidiskretizazio-sareetan oinarritutako metodoekin bezala (ingelesetik, multigrid methods) [37]. Gradienteei jaramonik egiten ez dien norma bat aukeratzen badugu, adibidez, L2-norma, euskarri lokaleko oinarriko funtzioei lotutako entrenamenduak emaitza bikainak ematen ditu (geruzaz geruzako entrenamendua egitean), 3.11a. irudian ikusten den moduan. Hala ere, ekuazio diferentzialak ebazten direnean ez da gomendagarria L2-norma erabilita optimizatzea. Nahiz eta hasieratik amaierarainoko entrenamendua aukerarik onena izan hondarren/erroreen gradientea darabilten galera-funtzioetan, galera-funtzioaren desbideraketa ikusten dugu iterazioek aurrera egin ahala (gogoratu 3.9b. irudia). Pentsatzen dugu hori normaren aukeraketarekiko independentea dela, baina aldagai entrenagarri askoren koexistentzia kontrajarriaren ondorioa dela. Gertaera hori argitzeko, L2-normaren goiko kasua aztertuko dugu, non geruzaz geruzako entrenamendua nahikoa baitzen konbergentzia bikaina lortzeko, eta hasieratik amaieraraionoko entrenamendua egingo dugu. 3.11b. irudiak galera-funtzioaren konbergentzia desbideraketa erakusten du. 51 3. Deep elementu finituen metodoa 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 10−13 10−6 101 s= 1 s= 2 s= 3 s= 4 Adam Adalr iterazioa galera galera = ∥AuNN −f∥A−1= errorearen energia-norma (a) Galera-funtzioaren eboluzioa lehenengo lau urratsetan. 0 5,000 10,000 15,000 20,000 25,000 30,000 10−13 10−6 101 iterazioa galera (b) Galera-funtzioaren eboluzioa hamaika urratsetan. 0 0.25 0.5 0.75 1 −1 −0.5 0·10−8 x uNN −uFEM (c) Errore-funtzioa entrenamenduaren amaieran. 3.9. irudia: Hasieratik amaierarainoko entrenamendua DeepFEMarentzako (3.25). probleman. Galera-funtzioa errorearen energia-norma da. 52 3.6. Zenbakizko esperimentuak 0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 10−13 10−6 101 iterazioa galera galera = ∥AuNN −f∥A−1= errorearen energia-norma (a) Galera-funtzioaren eboluzioa hamaika urratsetan. 0 0.25 0.5 0.75 1 −1.5 −1 −0.5 0·10−3 x uNN −uFEM (b) Errore-funtzioa entrenamenduaren amaieran. 3.10. irudia: Geruzaz geruzako entrenamendua DeepFEMean galera-funtzioan energia-norma erabiltzen denean (3.25). probleman. Jarraian, bakarrik hasieratik amaierarainoko entrenamenduak egingo ditugu. 3.6.1.2. Aurrebaldintzatzaileen ekintzak Orain, (3.19). ekuazioak adierazten duen galera-funtzioa aztertuko dugu hiru aurrebaldintzatzailerekin: (a) Pidentitate-matrizea izanik (3.12a. irudia); (b) Pblokekako Jacobi aurrebaldintzatzailea, bi tamainako blokeekin (3.12b. irudia); eta (c) Pblokekako Jacobi aurrebaldintzatzailea egokitutako tamainako blokeekin, haien tamaina diskretizazio-sarearen elementu kopuruaren erdia izanik (3.12c. irudia). Kasu guztietan, galera-funtzioaren eboluzioa eta errorearen energia-normarena erakusten ditugu. Zenbait ezberdintasun ikusten ditugu errorearen energia-normaren eta galerafuntzioaren artean; espero bezala, galera-funtzioa handitu egiten da diskretizaziosarea handitzen denean. Jacobiren blokea zenbat eta handiagoa izan, galerafuntzioaren eta errorearen normaren arteko ezberdintasuna txikiagoa da. Gaine53 3. Deep elementu finituen metodoa 0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 10−13 10−6 101 iterazioa galera galera = ∥A−1(AuNN −f)∥M= errorearen L2-norma (a) Galera-funtzioaren eboluzioa geruzaz geruzako entrenamenduak egiten direnean. 0 10,000 20,000 50,000 10−13 10−6 101 iterazioa galera (b) Galera-funtzioaren eboluzioa hasieratik amaierarainoko entrenamenduak egiten direnean. 3.11. irudia: DeepFEMaren entrenamenduak (3.25). problemarentzat galerafuntzioan L2-norma erabiltzean. ra, galera-funtzioa errorearen energia-normatik zenbat eta urrunago egon, galerafuntzioak arinago lortzen du konbergentzia (konparatu 3.12d eta 3.12e. irudietako erroreak). Horrek pentsarazten digu energia-normako erroreak induzitutako galera-funtzioa sinplifikatutako beste aldaera batzuk baino ganbilagoa dela aldagaiekiko, besteak beste hondar bektorearen 2-normak induzitutako galerafuntzioa baino ganbilagoa. 54 3.6. Zenbakizko esperimentuak galera = ∥AuNN −f∥Perrorearen energia-norma = ∥uNN −uFEM∥A 0 5,000 10,000 15,000 20,000 25,000 30,000 10−13 10−5 103 iterazioa galera/errorea (a) Aurrebaldintzatzaile barik, hau da, P=I. 0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 10−13 10−5 103 iterazioa galera/errorea (b) Bi tamainako blokeak erabilita. 0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 10−13 10−5 103 iterazioa galera/errorea (c) Diskretizazio-sarearen elementu kopuruaren erdiko tamainako blokeak erabilita. 0 0.25 0.5 0.75 1 0 0.5 1 1.5·10−2 x uNN −uFEM (d) (a). kasuaren errore-funtzioa. 0 0.25 0.5 0.75 1 −1.5 −1 −0.5 0·10−7 x (e) (c). kasuaren errore-funtzioa. 3.12. irudia: Hasieratik amaierarainoko DeepFEMaren entrenamenduak (3.25). problemarentzat tamaina ezberdineko blokedun Jacobi aurrebaldintzatzaileak erabilita. 55 3. Deep elementu finituen metodoa 3.6.1.3. Normaren trukaketa entrenamenduan zehar Aukeratzen dugun normaren arabera, galera-funtzioaren ganbiltasun mota ezberdinak agertuko dira DeepFEMaren aldagaiekiko. Sarri gertatzen da aldagaien domeinuaren zenbait zatitan galera-funtzio bat beste batzuk baino ganbilagoa izatea, eta horrek zuzeneko eragina du konbergentziaren optimizatzailean. Gelditzea saihesteko, eta ganbiltasuna hobetzeko, optimizazioan zehar norma aldatzea proposatzen dugu. Horretarako, esaterako, honako era honetan definituko dugu norma: L(σ, α;θ) = CE∥AuNN −f∥P+CL2∥P(AuNN −f)∥M,(3.26) non CE, CL2∈ {0,1}balio ezberdinak baitira, eta elkarrekin trukatuko dira konbergentzia gelditzea gertatzen baldin bada. Goian aipaturiko ideia argitzeko, 11b. irudiaren kasua aztertuko dugu. Galerafuntzio berbera mantenduko dugu lehenengo lau urratsetan. Gero, bosgarren urratsaren galera-funtzioa ezartzeko bi aldaera hartu ditugu kontuan: (a) aurreko urratseko galera-funtzio berbera (CE= 1 ezarriz (3.26). ekuazioan) baina 12.000 iteraziorekin; eta (b) (3.26). ekuazioan CE= 1 ezarriz 2.000 iteraziotan, gero CL2= 1 ezarri dugu beste 8.000 iteraziotan, eta azkenean CE= 1 ezarpenera bueltatu gara azkenengo 2.000 iterazioetan. Nahiz eta aldaera bietako iterazio kopuru osoa berdina izan, galera-funtzioaren balore txikiagoa lortzen dugu bosgarren urratsa amaitu ondoren norma aldatzen dugunean (10−10 inguru) aldaketarik egiten ez denean baino (10−8inguru) —ikus 3.13. irudian—. Galera-funtzioaren eboluzioaren malda altuagoa da L2-norma darabilgunean, eta horrek galera baxuagoko balio batetik abiatzea ahalbidetzen du energia-normara itzultzean, entrenamendu osoan zehar berbera mantentzen denean baino. 10,000 12,000 14,000 16,000 18,000 20,000 22,000 10−13 10−7 10−1 CE= 1 CL2= 1 CE= 1 iterazioa galera galera = ∥AuNN −f∥P galera = CE∥AuNN −f∥P+CL2∥P(AuNN −f)∥M 3.13. irudia: Energia-normaren entrenamendu luzea vs. energiaeta L2normaren arteko trukea 3.12b. irudiaren bosgarren (s= 5) urratsaren erdian. 56 3.6. Zenbakizko esperimentuak 3.6.2. Soluzio sinusoidalak Poisson and Helmholtz ekuazioetan Izan bitez ®−u′′ = 100π2sin(10πx), u(0) = 0, u′(1) = 10π, ®−u′′ −100π2u= 0, u(0) = 0, u′(1) = 10π. (3.27) Hemen, soluzio zehatz berdina duten bi problema ezberdin ditugu (lehenengoa simetrikoa eta positibo-definitua da, eta bigarrena, berriz, indefinitua): u∗(x) = sin(10πx).(3.28) Biak ebatziko ditugu DeepFEM erabilita, eta galera-funtzioa aurrebaldintzatutako hondarraren H1-norma bezala hautatuko dugu, hau da, L(θ;σ, α) = ∥P(AuNN −f)∥K+M.(3.29) 32 elementuko diskretizazio-sare batekin hasiko gara eta hiru fintze egingo ditugu. Paurrebaldintzatzailea hautatuko dugu 32 elementurekin eta 1 ≤s≤4 urratsetan. Optimizatzaileak eta NNaren arkitekturak 3.6.1. atalekoak dira. 3.14. eta 3.15. irudiek Poisson eta Helmholtz problemekin lotutako DeepFEM iragarpenak erakusten dituzte, hurrenez hurren. Poissonen, FEMaren soluzioek u∗-rekin bat egiten dute nodoetan, baina hori ez da gertatzen Helmholtz ekuazioaren kasuan. Aurrebaldintzatzaile gisa alderantzizko matrizea hartuz entrenatzen denean, DeepFEMak beherapen monotonikoa erakutsiz konbergitzen du. Hala ere, alderantzizkoa ez den aurrebanditzatzaile bat hartzen denean, konbergentzia gelditu egiten da pare bat magnitude-orden murriztu ondoren (ikus 3.16a. eta 3.16b. irudiak). Adalr optimizatzailearen entrenamendu-faseetan, galera-funtzioa batbatean murrizten da lehen iterazioetan, baina gero lautu egiten da, H1-normari dagokion errore murrizketa hutsala eginez. Bi esperimentuetan, 10−3inguruko erroreak lortzen ditugu (ikus 3.16c. eta 3.16d. irudiak). Emaitza horiek erakusten dute nolabaiteko zehaztasun-maila (baina ez handia) duten soluzio hurbilduak lor daitezkeela, ziur asko galerak dituen ganbiltasun ezagatik (gogoratu 2.3. atala). 57 3. Deep elementu finituen metodoa u∗uFEM uNN 0 0.25 0.5 0.75 1 −1 0 1 u(x) (a) s= 1 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −1 0 1 (b) s= 1 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −1 0 1 u(x) (c) s= 2 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −1 0 1 (d) s= 2 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −1 0 1 u(x) (e) s= 3 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −1 0 1 (f) s= 3 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −1 0 1 x u(x) (g) s= 4 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −1 0 1 x (h) s= 4 urratsa (entrenatu ondoren). 3.14. irudia: Poisson (3.27). probleman lortutako iragarpenak lau urratsetan zehar. u∗soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. 58 3.6. Zenbakizko esperimentuak u∗uFEM uNN 0 0.25 0.5 0.75 1 −2 0 2 u(x) (a) s= 1 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −2 0 2 (b) s= 1 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −2 0 2 u(x) (c) s= 2 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −2 0 2 (d) s= 2 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −1 0 1 u(x) (e) s= 3 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −1 0 1 (f) s= 3 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −1 0 1 x u(x) (g) s= 4 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −1 0 1 x (h) s= 4 urratsa (entrenatu ondoren). 3.15. irudia: Helmholtz (3.27). probleman lortutako iragarpenak lau urratsetan zehar. u∗soluzio zehatza da, uFEM elementu finituetako soluzioa, eta uNN DeepFEMaren iragarpena. 59 3. Deep elementu finituen metodoa α0 3 15 50 200 ∥A(α)uNN(α)−f∥P(α)0.0015 0.013 0.0046 0.012 0.082 ∥uNN(α)−uFEM(α)∥A(α)0.0160 0.014 0.0061 0.014 0.095 3.1. taula: Hondarraren eta errorearen energia-normaren balioak proba datubasearen lagin bakoitzerako entrenamenduaren amaieran, (3.32). probleman eta 0 < α < 200 denean. Zortzi tamainako blokeak erabili ditugu aurrebaldintzatzaileetan. 1., 2., 3.eta 4.urratsetan, hurrenez hurren), errorearen energia-norma eta galerafuntzioa antzera txikitzen dira (ikus 3.20b. irudia). 3.2. taulak galera-funtzioaren eta errorearen energia-normak erakusten ditu proba datuetan ebaluatuta, aurrebaldintzatutako entrenamenduaren amaieran. α0 3 15 50 200 ∥A(α)uNN(α)−f∥P(α)0.0022 0.0059 0.0043 0.0086 0.011 ∥uNN(α)−uFEM(α)∥A(α)0.0053 0.0075 0.0058 0.0101 0.013 3.2. taula: Hondarraren eta errorearen energia-normaren balioak proba datubasearen lagin bakoitzerako entrenamenduaren amaieran, (3.32). probleman eta 0 < α < 200 denean. Tamaina gorakorreko blokeak erabili ditugu aurrebaldintzatzaileetan. 3.6.4.2. Helmholtzen ekuazio parametrikoa: −50 < α < −30 Kasu honen soluzio analitikoa honako hau da: u∗(x) = Csin(√αx), C=2π √αcos(√α) izanik. Ehun lagineko datu-base bat aukeratuko dugu NNa entrenatzeko, eta honako probarako datu-base hau hautatuko dugu: {−50,−45,−40,−35,−30}.H1norma aukeratuko dugu galera-funtzioarentzako eta goian erabili dugun blokeka entrenagarria zen arkitektura bera. NNaren entrenamendua goiko bloke tamainen progresioa erabilita egiten badugu (ikus 3.21a. irudia), ohartuko gara oraingo honetan galerak ez duela errorearen H1-norma gutxitzen. Alderantzizko matrizeak aurrebaldintzatzaile gisa erabiltzen baditugu (ikus 3.21b. irudia), emaitzak onargarriak dira. 3.22. irudian urratsez urratseko iragarpenak aurkezten dira, entrenamendua 3.21b. irudiaren arabera egiten denean. Kontuan izan baliokideak direla alderantzizkoak erabiltzea eta honako galerafuntzio hau aplikatzea: L(θ, D) = 1 |D|X α∈D∥uNN(α)−uFEM(α)∥K+M,(3.33) 66 3.6. Zenbakizko esperimentuak galera = 1 |D|Pα∈D∥A(α)uNN(α)−f∥P(α) errorearen H1-norma = 1 |D|Pα∈D∥uNN(α)−uFEM(α)∥A(α) 0 2,000 4,000 6,000 8,000 10,000 10−2 10−1 100 iteration galera/errorea (a) Zortzi tamainako aurrebaldintzatzaileak erabilita. 0 2,000 4,000 6,000 8,000 10,000 10−2 10−1 100 iteration galera/errorea (b) 1., 2., 3., eta 4. urratsetan 8, 8, 16 eta 32 tamainako blokeak erabilita, hurrenez hurren. 3.20. irudia: Galera-funtzioaren eta errorearen energia-normaren eboluzioak lau urratsetan (3.32). probleman, 0 ≤α≤200 denean. 67 3. Deep elementu finituen metodoa galera = 1 |D|Pα∈D∥P(α){A(α)uNN(α)−f}∥K+M errorearen H1-norma = 1 |D|Pα∈D∥uNN(α)−uFEM(α)∥A(α) 0 2,000 4,000 6,000 8,000 10,000 10−2 10−1 100 101 iterazioa galera/errorea (a) 1., 2., 3., eta 4. urratsetan 8, 8, 16 eta 32 tamainako blokeak erabilita, hurrenez hurren. 0 2,000 4,000 6,000 8,000 10,000 10−2 10−1 100 101 iteration galera/errorea (b) Matrizeen alderantzizkoak erabilita aurrebaldintzatzaile gisa. 3.21. irudia: Galera-funtzioaren eta errorearen energia-normaren eboluzioak lau urratsetan (3.32). probleman,−50 < α < −30 denean. 68 3.6. Zenbakizko esperimentuak uFEM(−50) uFEM(−45) uFEM(−40) uFEM(−35) uFEM(−30) uNN(−50) uNN(−45) uNN(−40) uNN(−35) uNN(−30) 0 0.25 0.5 0.75 1 −2 −1 0 1 2 u(x) (a) s= 1 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 (b) s= 1 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 u(x) (c) s= 2 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 (d) s= 2 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 u(x) (e) s= 3 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 (f) s= 3 urratsa (entrenatu ondoren). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 x u(x) (g) s= 4 urratsa (entrenatu aurretik). 0 0.25 0.5 0.75 1 −2 −1 0 1 2 x (h) s= 4 urratsa (entrenatu ondoren). 3.22. irudia: Lau urratsetan lortutako iragarpenak proba datu-basean (3.32). probleman,−50 < α < −30 denean. Galera-funtzioan errorearen H1-norma erabili dugu. uFEM(α) elementu finituetako soluzioa da eta uNN(α) DeepFEMaren iragarpena αparametroarentzako. 69 3. Deep elementu finituen metodoa non uNN(α) ikurrak DeepFEMaren iragarpen bektoriala adierazten baitu eta uFEM(α) ikurrak FEMaren soluzio bektorea αparametroarentzako. Bektore horiek guztiak aurretiaz kalkulatuz, eta (3.33). ekuazioan dagoen galera-funtzioa erabiliz, entrenamendua alderantzizko matrizeak erabilita baino askoz merkeagoa izango litzateke. Sinpletasunagatik, kapitulu honetan zehar alderantzizkoak kalkulatu ditugu aurkeztutako arrazoiketarekin koherentzia mantentzeko, eta ez dugu exekuzio denborarik konparatu; ebazle eraginkor batek inoiz ez bailuke matrize baten alderantzizkoa esplizituki kalkulatuko [5,169,158]. 70 4. The Deep Double Ritz Method Summary. Residual minimization is a widely used technique for solving Partial Differential Equations in variational form. It minimizes the dual norm of the residual, which naturally yields a saddle-point (min–max) problem over the so-called trial and test spaces. In the context of Neural Networks, we can address this min–max approach by employing one Neural Network to seek the trial minimum while another Neural Network seeks the test maximizers. However, the resulting method is numerically unstable as we approach the trial solution. To overcome this, we reformulate the residual minimization as an equivalent minimization of a Ritz functional fed by optimal test functions computed from another Ritz functional minimization. We call the resulting scheme the Deep Double Ritz Method, which combines two Neural Networks for approximating trial and optimal test functions along a nested double Ritz minimization strategy. Numerical results on different diffusion and convection problems support the robustness of our method up to the approximation properties of the considered Neural Networks and the training capacity of the optimizers. Refer to [217] for the published version. 4.1. Introduction Within the variational framework introduced in Section A.5.2/1.5.2, residual minimization reads as a saddle-point (min-max) problem as follows: min u∈Umax v∈V\{0}⟨Bu −l, v⟩V′×V ∥v∥V ,(4.1) where Uand Vare the trial and test spaces, B:U−→ V′is the differential operator governing the considered BVP in variational form, V′is the dual space of V, and l∈V′is the right-hand side. In [232,20], the authors proposed addressing this optimization scheme employing Generative Adversarial Networks (GANs) [78,79] by approximating uand vwith two NNs. Unfortunately, this approach presents a severe numerical limitation: the Lipschitz continuity constant of the test maximizers with respect to the trial functions might become arbitrarily 71 4. The Deep Double Ritz Method large when approaching the exact solution. Moreover, the corresponding test maximizer is highly non-unique in the limit. In consequence, we end up with an inherent lack of numerical stability that is easily confirmed by numerical experiments with simple model problems. To overcome the above limitations, we reformulate residual minimization as a minimization of a Ritz functional fed by optimal test functions [55,58]. Since optimal test functions are generally unavailable, we compute them for each trial function using another Ritz method. Thus, the resulting scheme is a nested double-loop Ritz minimization method: the outer loop seeks the trial solution, while the inner loop seeks the optimal test function for each trial function. We call the resulting scheme the Double Ritz Method. In some occasions, the trial-to-test operator that maps each trial function with the corresponding optimal test function is available. For example, when the problem is symmetric and positive-definite, and we consider the norm induced by the bilinear form for the trial and test spaces, the trial-to-test operator is the identity; or when selecting the strong variational formulation, the trial-to-test operator is the one given by the PDE operator. In these cases, the Double Ritz Method reduces to a single-loop Ritz minimization (this will be further discussed in Section 4.2.4). Thus, the Double Ritz Method is a general method for solving PDEs in different variational forms, which in some particular cases simplifies into a single-loop Ritz minimization method. Thanks to NNs, we find a simple and advantageous computational framework to approximate and connect the trial and test functions between the innerand outer-loop minimizations in the Double Ritz Method, a task that is challenging to tackle with traditional numerical methods. Herein, we propose using one NN to represent the trial functions and another NN to represent the local actions of the trial-to-test operator. Thus, the composition of both NNs represents the (optimal) test functions, and we preserve the trial dependence of the test functions during the entire process. We call the resulting NN-based method the Deep Double Ritz Method (D2RM). While the D2RM replicates existing residual minimization methods in the context of NNs, we fall short of providing a detailed mathematical convergence analysis due to the manifold structure of NNs that departs from the traditional vector-space-based mathematical approach (recall Chapter 2). Related to this, we encounter the usual drawbacks of lack of convexity between the objective/loss function with respect to the trainable parameters, which prevents us from making a proper diagnosis of the optimizer during training. The remainder of this chapter is organized as follows. Section 4.2 formalizes the variational setting introduced in Section A.5.2 and derives the Double Ritz Method at the continuum level. Section 4.3 introduces the D2RM within the NN 72 4.2. From residual to Ritz minimizations framework. Section 4.4 provides implementation details of the addressed methods and Section 4.5 develops on numerical experimentation. 4.2. From residual to Ritz minimizations We introduce the residual minimization framework, followed by a saddle-point reformulation and an alternative Double Ritz scheme at the continuum level. Subsequently, we describe three particular cases for which the Double Ritz method simplifies into a single Ritz method. 4.2.1. Residual minimization Let ßFind u∗∈Usuch that b(u∗, v) = l(v),∀v∈V,(4.2a) where Uand Vare real Hilbert trial and test spaces, respectively, b:U×V−→ R is a bilinear form, and l:V−→ Ris a continuous linear functional. Equivalently, in operator form: ßFind u∗∈Usuch that Bu∗=l,(4.2b) where B:U−→ V′is the operator defined by ⟨Bu, v⟩V′×V:= b(u, v), V′denotes the topological dual of V, and l∈V′. The equivalent residual minimization formulation reads as u∗= arg min u∈U∥Bu −l∥V′,(4.3) where Bu −l∈V′is the residual for each trial function u∈U. To guarantee well-posedness of (4.2)–(4.3), we assume the hypotheses of the Babuˇska–Lax–Milgram Theorem [17], which according to our presentation translate into that Bis an isomorphism that is bounded from above and below, i.e., there exist some positive constants γ≤Msuch that it satisfies γ∥u∥U≤ ∥Bu∥V′≤M∥u∥U, u ∈U.(4.4) Then, the error in Uis equivalent to the residual in V′in the following sense: 1 M∥Bu −l∥V′≤ ∥u−u∗∥U≤1 γ∥Bu −l∥V′, u ∈U.(4.5) Below, we examine two alternatives to evaluate and minimize the residual in its dual norm. 73 4. The Deep Double Ritz Method 4.2.2. Saddle-point problem (min-max optimization) The norm in V′is defined in terms of the norm in Vas ∥t∥V′:= sup v∈V\{0} ⟨t, v⟩V ∥v∥V = sup v∈V\{0}≠t, v ∥v∥V∑V′×V , t ∈V′.(4.6) Hence, combining (4.3) and (4.6) yields u∗= arg min u∈Umax v∈V\{0}Fmax min (u, v) (4.7a) with Fmax min (u, v) := ≠Bu −l, v ∥v∥V∑V′×V =bÅu, v ∥v∥Vã−lÅv ∥v∥Vã.(4.7b) When Fmax min is fed by the exact solution u∗, the resulting operator becomes the null functional, i.e., Fmax min (u∗,·)=0∈V′, which implies a highly non-uniqueness of the test maximizer in the limit—indeed, any element in V\{0}is a test maximizer for u∗. Outside this singular case, the operator that maps each trial function u∈U to its unitary test maximizer is well-defined but not Lipschitz continuous when approaching the exact solution u∗, leading to an unstable numerical method. We formalize this in the following two items: •Let vmax :U\ {u∗} −→ Vbe the mapping that for each trial function returns the test maximizer of the actions of the residual Bu −l∈V′over the unitary sphere, i.e., vmax(u) := arg max ∥v∥V=1⟨Bu −l, v⟩V′×V.(4.8) Then, vmax(u)∈Vis unique for each u∈U\{u∗}. Proof. Let u∈U\{u∗}. By the Riesz Representation Theorem, there exists a unique ru∈Vsuch that (ru, v)V=⟨Bu −l, v⟩V′×V,∀v∈V,(4.9) and ∥ru∥V=∥Bu −l∥V′= max ∥v∥V=1⟨Bu −l, v⟩V′×V.(4.10) Let vmax =vmax(u) be a test maximizer of (4.10) in the unitary sphere of V, i.e., ⟨Bu −l, vmax⟩V′×V= max ∥v∥V=1⟨Bu −l, v⟩V′×V.(4.11) 74 4.2. From residual to Ritz minimizations By the Cauchy-Schwarz inequality: 0<∥ru∥V=∥Bu −l∥V′=⟨Bu −l, vmax⟩V′×V= (ru, vmax)V(4.12a) ≤ ∥ru∥V∥vmax∥V,(4.12b) where the equality holds if and only if vmax =λrufor some λ > 0. Then, vmax =ru ∥ru∥Vis unique. •vmax is not Lipschitz continuous around any (reduced) neighborhood of u∗, i.e., there does not exist a constant 0 < C < ∞such that ∥vmax(u1)−vmax(u2)∥V≤C∥u1−u2∥U,∀u1, u2∈U\{u∗}.(4.13) Proof. Assume by contradiction that there exists 0 < C < ∞such that (4.13) holds, and let u2= 2u∗−u1with u1=u∗. Then, 2 = ∥vmax(u1)−vmax(u2)∥V≤C∥u1−u2∥U= 2C∥u1−u∗∥U.(4.14) Letting u1→u∗, we obtain C→ ∞. 4.2.3. Double Ritz Method with optimal test functions The Riesz Representation Theorem allows us to work isometrically in the test space instead of in its dual. In particular, for the residual, we have ∥Bu −l∥V′=∥R−1 V(Bu −l)∥V, u ∈U,(4.15) where RV:V∋v7−→ (v, ·)V∈V′denotes the Riesz operator. This relation suggests considering the trial-to-test operator T:U−→ Vdefined by T:= R−1 VB [58] since it relates the error in Uwith the Riesz representative of the residual in V, T(u−u∗) = R−1 V(Bu −l)∈V, u ∈U.(4.16) The images of trial functions through Tare known as optimal test functions [55], and they allow us to rewrite (4.2) in terms of the following symmetric and positive-definite variational problem: ßFind u∗∈Usuch that (Tu∗, Tu)V=l(Tu),∀u∈U,(4.17) Proof. On the one hand, b(u, v) = ⟨Bu, v⟩V′×V= (Tu, v)Vfor all u∈Uand all v∈V. On the other hand, Tis an isomorphism because so is B. Hence, testing with all v∈Vis equivalent to testing with Tu ∈Vfor all u∈U. 75 4. The Deep Double Ritz Method Algorithm 5: Training of the Deep Double Ritz Method (D2RM) Initialize θu∈Θuand θτ∈Θτ; /* Outer-loop */ while not converged do Randomly sample {xj}N j=1 ⊂Ω; θu=θu−λu ∂LτNN ∂θu (θu, θτ;{xj}N j=1); /* Inner-loop */ while not converged do Randomly sample {xj}N j=1 ⊂Ω; θτ=θτ−λτ ∂Lopt uNN ∂θτ (θu, θτ;{xj}N j=1); return θu Algorithm 6: Training of the Generalized Deep Ritz Method (GDRM) Initialize θu∈Θu; while not converged do Randomly sample {xj}N j=1 ⊂Ω; θu=θu−λu ∂LT ∂θu (θu;{xj}N j=1); return θu and outer loops. We establish an accumulated maximum number of iterations for both nested loops, and we fix four inner-loop iterations for each outer-loop iteration4(as considered in [232,20] for WANs5) unless otherwise specified. 4We made this decision motivated by the continuity of the elements involved during the process: a slight modification in θutranslates into a slight variation in uNN (continuity of the realization mapping—recall Section 2.1) that produces a slight variation in T uNN (continuity of T). Therefore, we expect (have the hope) that τNN(uNN) can provide a good approximation of T uNN after a small number of iterations in θτ. 5The lack of Lipschitz continuity in the min-max approach suggests that this is a poor strategy for optimizing WANs. 82 4.4. Implementation 4.4. Implementation We use the TensorFlow 2 (TF2) library [3,2] within Python to implement our neural network architectures and loss functions, and to manage the random creation and flow of data. Specifically, we accommodate all of our implementations to Keras (tf.keras). 4.4.1. Samples generation, input batch flow, encapsulation of models, and optimization Our inputs to networks are samples over the domain Ω. After feeding our networks with the batch of inputs, we combine the batch of outputs in a single loss prediction. We encapsulate the networks and losses inside a main model whose input and output are the batch of samples and the loss prediction, respectively. From the loss prediction, we optimize the learnable parameters of the networks (i.e., we fit the learnable parameters to the data) with the Keras built-in Adam optimizer. Figure 4.1 illustrates the described process at each training iteration. random batch generator of samples from Ω x1 x2 x3 x4 x5 . . . xN Main model loss prediction optimizer readjust learnable parameters Figure 4.1: Implementation sketch of the general flux for the proposed methods at each training iteration. For the maximization involved in WANs, we reverse the sign of the gradients so that when feeding the (default) gradient-descent-based optimizer, it instead performs a gradient ascent. 83 4. The Deep Double Ritz Method 4.4.2. Design of models by methods We implement networks, losses, and operators as models and layers in TF2 from the redefinitions of the corresponding Keras base classes (tf.keras.Model and tf.keras.Layer). In WANs, we implement uNN and vNN as two independent models that are combined via a common non-trainable layer for the loss (see Figure 4.2). In the DRM, we implement uNN as a model that subsequently connects with two non-trainable layers for the trial-to-test operator and the loss function, respectively (see Figure 4.3). In the D2RM, we implement uNN and τNN as two sequential models whose output feeds into two separate losses (see Figure 4.4). The loss functions are implemented as latent outputs of the main model, which, together with a TF2-suitable boolean variable, activate and deactivate alternatively. Despite the significant compilation time that the twobranch model (for the D2RM) takes compared with one-branch models (for WANs and the DRM)6, its fitting execution in graph mode is as fast as that of one-branch models. Main model x1 x2 x3 x4 x5 . . . xN . . .. . . uNN(·;θu) : Ω −→ Ru1 u2 . . . uN uj:= uNN(xj;θu) Lmax min . . .. . . v1 v2 . . . vNvj:= vNN(xj;θv) vNN(·;θv):Ω−→ R non-trainable Lmax min θu, θv;{xj}N j=1 Figure 4.2: Main model architecture for WANs. It consists of two independent NNs, uNN and vNN, combined via the loss function Lmax min . 4.4.3. Graph-mode execution dynamics We employ callbacks to avoid interrupting the graph execution mode carried out by the Keras fitting instruction (.fit). Callbacks act during fitting and enable 6Around one minute for the D2RM vs. a couple of seconds for WANs and the DRM. 84 4.4. Implementation Main model x1 x2 x3 x4 x5 . . . xN . . .. . . uNN(·;θu) : Ω −→ R u1 u2 u3 u4 u5 . . . uN uj:= uNN(xj;θu) T non-trainable T1 T2 T3 T4 T5 . . . TN Tj:= T(uj) LLTθu;{xj}N j=1 non-trainable Figure 4.3: Main model architecture for the DRM. It consists of a NN, uNN, composed with the trial-to-test operator Tand the loss function L. Main model x1 x2 x3 x4 x5 . . . xN . . .. . . uNN(·;θu) : Ω −→ R u1 u2 u3 u4 u5 . . . uNuj:= uNN(xj;θu) . . .. . . τNN(·;θτ) : uNN(Ω; θu)−→ R τ1 τ2 τ3 τ4 τ5 . . . τNτj:= τNN(uj;θτ) LτNN Lopt uNN LτNN θu;{xj}N j=1 Lopt uNN θτ;{xj}N j=1 non-trainable if uNN trainable if τNN trainable Figure 4.4: Main model architecture for the D2RM. It consists of two NNs, uNN and τNN, equipped with the loss functions LτNN and Lopt uNN . accessing certain elements of (main) models and modifying them. We utilize callbacks for loss monitoring (for WANs, the DRM, and the D2RM), to activate or deactivate the trainability of the networks and switch between optimizers (for WANs and the D2RM), or to interchange losses (for the D2RM) during training when iterating either over the outeror the inner-loop. 85 4. The Deep Double Ritz Method 4.5. Numerical results We show numerical experiments to compare the methods introduced above. Section 4.5.1 considers a simple model problem and compares WANs, the DRM, and the D2RM. Section 4.5.2 makes a more profound comparison considering a parametric model problem. Section 4.5.3 and Section 4.5.4 consider sources that lead to singular problems in pure diffusion and convection equations, respectively. Finally, Section 4.5.5 considers a pure convection equation in 2D. 4.5.1. Initial comparison of WANs, the DRM, and the D2RM on a simple problem We select model problem (4.24) in weak form with source f=−2, so the exact solution is u∗=x(x−1). Hence, U=H1 0(0,1) = Vand Tis the identity operator. We solve this problem by employing WANs, the DRM, and the D2RM. We select a two-layer fully-connected NN with 20 neurons on each layer and tanh activation functions for the architectures of uNN,vNN, and τNN. We perform 200 iterations for uNN in the three methods: WANs, the DRM, and the D2RM. Since in WANs and the D2RM we established four iterations to approximate the test maximizers, we end up with a total of 1,000 training iterations (200 for the trial function, and 800 for the test functions). Moreover, we select batches of size 200 for the training and a uniform distribution for the sample generation. Figure 4.5 shows the uNN network predictions and errors of the three methods at the end of the training. We observe that WANs produce a larger error than the DRM and the D2RM. Figure 4.6 shows the loss evolution for WANs. Every five iterations, the loss decreases (minimization in uNN) and increases in the remaining iterations (maximization in vNN with fixed uNN). From iteration 500 onwards, the loss stops improving and oscillates above the optimal value. Figure 4.7 shows the evolution of the loss for the DRM. Here, we have a single minimization, so the observed noisy behavior of the loss towards the end of the training is attributed to the optimizer performance. Due to the lower complexity of the training, we achieve better convergence performance than with WANs. Figure 4.8 shows the loss evolutions for the D2RM. Here, we have a nested min-min optimization. At each iteration, we evaluate both LτNN and Lopt uNN , even if we are only optimizing with respect to one of them. We superimpose both losses, each one with its own scale (the left vertical axis corresponds to LτNN and the right vertical axis to Lopt uNN ). Both losses exhibit a decreasing staircase shape with downward-sloping steps. Jumps occur when the optimization is performed with respect to LτNN , which suggests that Lopt uNN takes longer to converge as it depends on the convergence of LτNN (innervs. outer-loop). 86 4.5. Numerical results u∗uNN WANs uNN DRM uNN D2RM 0 0.33 0.66 1 −0.2 −0.1 0 x uNN 0 0.33 0.66 1 −0.2 −0.1 0 x 0 0.33 0.66 1 −0.2 −0.1 0 x 0 0.33 0.66 1 −6 −4 −2 0 2·10−3 x uNN −u∗ 0 0.33 0.66 1 −10 −5 0 5·10−4 x 0 0.33 0.66 1 −10 0 10 ·10−4 x Figure 4.5: Trial network predictions and errors in WANs, the DRM, and the D2RM at the end of the training in model problem (4.24) with exact solution u∗=x(x−1). Lmax min (uNN, vNN)0 0 100 200 300 400 500 600 700 800 900 1,000 0 0.2 0.4 iteration loss 0 20 40 60 80 100 0.3 0.4 iteration loss 500 750 1,000 0 0.5 1·10−2 iteration Figure 4.6: Loss evolution during the WANs training for model problem (4.24) with exact solution u∗=x(x−1). 87 4. The Deep Double Ritz Method LT(uNN)FT(u∗) = −1/6 0 50 100 150 200 −0.15 −0.1 −0.05 iteration loss 160 180 200 −0.16668 −0.16666 −0.16664 −0.16662 iteration Figure 4.7: Loss evolution during the DRM training in model problem (4.24) with exact solution u∗=x(x−1). LτNN (uNN)Lopt uNN (τNN(uNN)) Fopt u∗(Tu∗) = FT(u∗) = −1/6 0100 200 300 400 500 600 700 800 900 1,000 −0.15 −0.1 −0.05 0 iteration loss 0 10 20 30 40 50 −6 −4 −2 0·10−2 minimizing LτNN minimizing Lopt uNN iteration loss −3 −2 −1 0 ·10−2 500 750 1,000 −0.1675 −0.167 −0.1665 −0.166 iteration Figure 4.8: Loss evolution of the D2RM training for model problem (4.24) with exact solution u∗=x(x−1). 88 4.5. Numerical results The relative errors of the trial network predictions7at the end of the training are 3.12%, 0.99%, and 1.31% in WANs, the DRM, and the D2RM, respectively. 4.5.2. Comparison of WANs, the DRM, and the D2RM on smooth problems Now, we focus on the evolution of the relative error for the previous weak formulation of (4.24), but selecting the source so that the solution varies according to a parameter: u∗ α=xα(x−1) ∈H1 0(0,1), α > 1/2.(4.28) 4.5.2.1. Without singularities: α≥1 We experiment individually for α∈ {2,5,10}with 5,000 training iterations for uNN in WANs, the DRM, and the D2RM. Note that this corresponds to a total of 25,000 iterations for WANs and the D2RM when taking into account the iterations dedicated to the test maximizers. Table 4.2 displays the relative errors along different stages of the training, and Figure 4.9 shows the trial network predictions and error functions at the end of training. WANs show poor results in all the cases, with a clear non-convergent tendency as the training progresses, possibly justified by the unstable behavior of the method at the continuous level. Thus, we discard the WANs for the remainder of the experiments and focus on the other two methods. In the DRM and the D2RM, we observe a decreasing behavior of the relative error during training. We highlight the nearly identical behavior of the relative norm errors of uNN and τNN(uNN) in the D2RM, which suggests that the D2RM behaves as the DRM, as desired. 7To approximate ∥uNN −u ∗∥U ∥u∗∥U×100, we perform a composite intermediate-point rule with 104 integration nodes for the numerator, and analytically calculate ∥u∗∥U=√3/3 for the denominator. 89 4. The Deep Double Ritz Method Training progress 4%20%40%60%100% Method α∥uNN−u∗∥U ∥u∗∥U×100 WANs 2 2.49% 3.43% 5.92% 7.40% 10.07% 5 58.69% 41.78% 63.03% 74.06% 40.33% 10 93.15% 93.95% 89.04% 68.85% 367.61% DRM 2 3.40% 2.47% 1.17% 0.30% 0.23% 5 50.50% 8.53% 2.83% 2.27% 1.59% 10 58.55% 9.51% 3.39% 2.83% 1.69% D2RM 2 3.27% 1.84% 0.31% 0.27% 0.54% 5 59.13% 12.93% 2.92% 2.52% 1.56% 10 82.31% 20.18% 6.24% 2.68% 2.60% Method α∥τNN(uNN)−T u∗∥V ∥Tu∗∥V×100 D2RM 2 3.36% 1.94% 0.48% 0.45% 0.58% 5 59.13% 12.93% 2.90% 2.11% 1.55% 10 83.51% 20.18% 6.23% 2.68% 2.61% Table 4.2: Relative errors of uNN (in WANs, the DRM, and the D2RM) and τNN(uNN) (in the D2RM) along different stages of the training progress in problem (4.24) with exact solution u∗ α=xα(x−1) and α∈ {2,5,10}. 90 4.5. Numerical results u∗ 2u∗ 5u∗ 10 uNN,2uNN,5uNN,10 0 0.33 0.66 1 −1 0 1·10−1 x uNN,α 0 0.33 0.66 1 0 −1 ·10−1 x 0 0.33 0.66 1 0 −1 ·10−1 x 0 0.33 0.66 1 −1 0 1·10−1 x uNN,α −u∗ α 0 0.33 0.66 1 −4 −2 0 2·10−4 x 0 0.33 0.66 1 −1 0 1·10−3 x 0 0.33 0.66 1 −1 0 1 x (uNN,α −u∗ α)′ (a) With WANs. 0 0.33 0.66 1 −2 −1 0 1·10−2 x (b) With the DRM. 0 0.33 0.66 1 −3 −2 −1 0 1·10−2 x (c) With the D2RM. Figure 4.9: Trial network predictions and errors for WANs, the DRM, and the D2RM in model problem (4.24) with exact solution u∗ α=xα(x−1) for α∈ {2,5,10}. 91 4. The Deep Double Ritz Method Training progress 4%20%40%60%100% ∥vNN−Tu∗∥V ∥Tu∗∥V×100 19.03% 4.58% 3.55% 3.20% 2.83% ∥A′vNN−u∗∥U ∥u∗∥U×100 21.95% 5.28% 4.10% 3.69% 3.27% Table 4.5: Relative errors of vNN and A′vNN along different stages of the training progress in problem (4.30) with exact solution (4.33). Training progress 4%20%40%60%100% ∥uNN−u∗∥U ∥u∗∥U42.05% 29.69% 18.42% 11.20% 8.93% ∥τNN(uNN)−T u∗∥V ∥Tu∗∥V37.27% 28.25% 19.02% 9.55% 6.13% Table 4.6: Relative errors of uNN and τNN(uNN) along different stages of the training progress in problem (4.30) with exact solution (4.33). u∗uNN uNN −u∗ Tu∗τNN(uNN)τNN(uNN)−Tu∗ 0 0.33 0.66 1 0 0.5 1 x uNN 0 0.33 0.66 1 −0.5 0 0.5 x uNN −u∗ 0 0.33 0.66 1 0 0.2 0.4 x τNN(uNN) 0 0.33 0.66 1 −1 0 ·10−3 x τNN(uNN)−Tu∗ 0 0.33 0.66 1 −0.5 0 0.5 x (τNN(uNN)−Tu∗)′ Figure 4.13: uNN and τNN(uNN) predictions, errors, and derivative of the errors for the D2RM in problem (4.30) with exact solution (4.33). 98 4.5. Numerical results 4.5.5. Pure convection in 2D Let    ∂u ∂x +∂u ∂y =kπ sin (kπ(x+y)) ,in Ω = (0,1) ×(0,1), u(x, 0) = u(0, y) = 0,0≤x, y ≤1, (4.35) and consider its strong variational formulation b(u, v) := ZΩÅ∂u ∂x +∂u ∂y ãv, l(v) := kπ ZΩ sin(kπ(x+y))v, (4.36) for u∈Uand v∈Vsuch that U=ßu∈L2(Ω) : ∂u ∂x +∂u ∂y ∈L2(Ω) and u(x, y) = 0 when xy = 0™,(4.37a) V=L2(Ω).(4.37b) Its exact solution is u∗= sin(kx) sin(ky) and the trial-to-test operator is the PDE operator. For k= 3/2, we perform 200,000 iterations in the D2RM with a training regime of nine iterations in the inner loop for each iteration in the outer loop. Exceptionally, we run the inner loop for 2,000 iterations before the first outer-loop iteration. For integration, we consider 50 nodes on each axis (i.e., 250 integration points on the entire domain due to the cartesian-product structure). We increase the NN architecture to three layers of 50-neuron width. Figure 4.14 shows the trial and optimal test predictions with corresponding error functions at the end of the training. The resulting relative errors are 3.62% and 1.70% for uNN and τNN(uNN), respectively. 99 4. The Deep Double Ritz Method 00.33 0.66 100.330.66 1 −1 0 1 xy uNN −1−0.5 0 0.5 00.33 0.66 100.330.66 1 −5 0 5 xy τNN(uNN) −4−2 0 2 4 00.33 0.66 100.330.66 1 −4 −2 0 ·10−2 xy uNN −u∗ −2.5·10−20 00.33 0.66 100.330.66 1 −0.2 0 0.2 xy τNN(uNN)−Tu∗ −0.25 0 0.15 00.33 0.66 100.330.66 1 0 0.2 xy |∇(uNN −u∗)| Figure 4.14: Trial and optimal test predictions, errors, and derivative of the trial error for the D2RM in model problem (4.30) with exact solution (4.33). 100 5. Memorian onarritutako Monte Carlo integrazioa∗ Laburpena. Monte Carlo integrazioa oso erabilia den koadraturaerregela bat da neurona-sareen bidez deribatu partzialetako ekuazioak ebazten direnean, gain-doikuntza gabeko hurbilpenak bermatzen dituelako eta dimentsio handiko eskalagarritasun gaitasuna duelako. Hala ere, metodo estokastiko horrek galerak eta gradiente zaratatsuak eragiten ditu entrenamenduan zehar, konbergentzia-diagnostiko egokia zailago eginez. Normalean, hori gainditzeko integrazio-puntu kopuru izugarria erabiltzen da, nahiz eta horrek entrenamenduaren errendimendua hondatu. Lan honetan, memoria oinarri duen Monte Carlo integrazio-metodo bat proposatzen da, koadratura-erregela zehatza sortzea ahalbidetzen duena eta entrenamendurako ez duena behar lagin handien prozesakuntzak eragiten duen kostu konputazional garestirik. Argitaratutako bertsiorako, ikus [218]. 5.1. Sarrera NNaren portaera txarra ekar dezake integrazio-puntuetatik urrun integrazio-puntu finkoak dituen koadratura-erregela deterministiko bat erabiltzeak, gaindoikuntzako fenomenoagatik [186]. Horrek, integrazio-errore esanguratsuak eragiten ditu eta, ondorioz, pobreak dira lortzen diren hurbilpenak. Hori gainditzeko, egokia da Monte Carlo (MC) integrazioa koadratura-erregelaren aukeraketa, laginketa estokastikoa eta sarerik gabeko egitura baititu [128,84,62,48]. Hala ere, Monte Carlo integrazio-errorea honako ordena honetakoa da: O(1/√N), non Nintegrazio-puntuen kopurua baita [163]. Orduan, praktikan, hamarnaka mila edo ehunka mila integrazio-puntu behar dira MC integrazioaren bidez erroremaila onargarri bat lortzeko (baita dimentsio bakarreko problemetarako ere), entrenamendu-abiaduraren kaltetan. Lan honetan, memorian oinarritutako koadratura-erregela bat proposatuko dugu, NNak inplikatzen dituen integral definituak hurbiltzen dituena eta entrenamenduan zehar aurreko iterazioetan lortutako informazioa erabiliko duena. ∗See Appendix C for the English version/Ikus C. eranskina ingelesezko bertsiorako. 101 5. Memorian onarritutako Monte Carlo integrazioa Integral horietatik espero den balioa nabarmen aldatzen ez den bitartean, teknika honek espero den integrazio-errorea murriztu egiten du eta koadratura-erregela hobeak eragiten ditu. Gainera, gradientea ere integral definitu baten bidez adieraz daitekeenez, ikuspegi hori gradientearen konputazioei ere aplikatuko diegu, optimizatzailearen hiperparametroak behar bezala doituz. Horrela, momentum metodo ezagunaren berrinterpretazio bat lortuko dugu [177]. Kapitulu honen gainerako partea honela dago antolatuta. 5.2. atalean hurbilketaeta optimizazio-esparruak berrikusiko dira, 5.3. atalean memorian oinarritutako integrazioeta optimizazio-estrategiak proposatuko dira, eta 5.4. atalean proposamen horiek momentum metodoarekin lotuko dira. 5.2. Hurbiltzeko, diskretizatzeko, parametrizatzeko eta optimizatzeko konfigurazioen berrikuspena 2. kapituluan landutako kontzeptu asko laburbilduko ditugu notazioa (bir)definitzeko, eta irakurlearentzat modu arinagoan eta errazagoan aurkeztuko dugu metodoa. Izan bedi u∗= arg min u∈UF(u),(5.1) ondo definituta dagoen minimizazio problema bat, non U-k funtzioen bilaketaespazioa adierazten baitu, Ω domeinua dutenak. F:U−→ Rminimizazio problema zuzentzen duen funtzio objektiboa da eta u∗haren soluzio zehatz bakarra da. Izan bedi uθ: Ω −→ Rneurona-sareko arkitektura bat, hurbilpenaren parametrizazioa adierazten duena, non θ∈Θ entrenagarriak diren parametroen multzoa baita (Θ-k parametroen domeinua adierazten du). Orduan, (5.1) problemaren hurbilpena egiteko uθ-ren bidez, Ubilaketa-espazioren ordez UΘ={uθ:θ∈Θ} espazio parametrizatura mugatuko gara. Minimizazioa gauzatzeko, lehen mailako gradiente-jaitsiera eskema bat erabiliko dugu. Hona hemen eskema hori: θt+1 =θt−λ∂F ∂θ (vθt),(5.2) non λ > 0 ikaskuntza-tasa baita eta θt-k entrenagarriak diren parametroen multzoa adierazten baitu t. iterazioan. Funtzio objetiboak, F-k, integral definituaren itxura badu, hau da, F(vθ) = ZΩ I(vθ)(x)dx, (5.3) 102 5.3. Memorian oinarritutako integrazioa eta optimizazioa orduan koadratura-erregela baten bidez hurbilduko dugu, Lgalera-funtzioa sortuz. MC koadratura-erregela hautatzen badugu, galera-funtzioak honako itxura hau du: F(vθ)≈ L(vθ) := Vol(Ω) N N X i=1 I(vθ)(xi),(5.4) non {xi}N i=1-k integrazio-puntuen multzoa adierazten baitu, Ω-ko zorizko banaketa uniforme batetik lagindutakoa. Era berean, gradienteei dagokienez, honako itxura hau dugu: ∂F ∂θ (vθ)≈g(vθ) := ∂L ∂θ (vθ) = Vol(Ω) N N X i=1 ∂I(vθ) ∂θ (xi).(5.5) Honela, (5.5)-ek (5.2)-ren bertsio diskretizatu bat adierazten du, gradiente jaitsiera estokastikoa (SGDa)2deitua [187], eta honela deskribatuta dagoena: θt+1 := θt−λg(vθt).(5.6) Hemendik aurrera, F(θt), ∂F ∂θ (θt), L(θt) eta g(θt) idatziko dugu, sinpletasunagatik, F(vθt)-ren, ∂F ∂θ (vθt)-ren, L(vθt)-ren eta g(vθt)-ren ordez, hurrenez hurren3. 5.3. Memorian oinarritutako integrazioa eta optimizazioa Sarea (5.6)-ren arabera entrenatzen badugu, portaera zaratatsua eta oszilakorra lortzen dugu galeran eta gradientean. Hori gertatzen da Monte Carlo integrazioak entrenamenduko iterazio bakoitzean txertatzen duen erroreagatik. 5.1. irudiak Monte Carlo integrazioaren portaera zaratatsua erakusten du (kurba urdina), Fzehazki kalkulatzea ahalbidetzen duen parametro entrenagarri bakarreko sare batean (kurba beltza). Atal hau hiru zatitan antolatuta dago. 5.3.1. atalean esperimentaziorako aukeratutako eredu problema aurkezten da, 5.3.2. atalean proposatutako memorian oinarritutako MC integrazio-araua deskribatzen da, eta 5.3.3. atalean memorian oinarritutako eskema hori gradienteetara hedatzen da entrenamendurako. 2Ingelesetik, Stochastic Gradient Descent (SGD). 3Ohartu errealizazio-maparen terminologia alde batera utzi dugula —gogoratu (2.5)— parametroen funtzio objektiboaren mendekotasuna adieraztean. Hau da, F(θ)-k (F ◦ ΦNN)(θ) ordezkatzen du, 2. kapituluan garatutako notazioaren arabera. 103 5. Memorian onarritutako Monte Carlo integrazioa 5.3.1. Eredu problema Izan bedi −u′′ = 4δ1/2, Ω = (0,1) domeinuan eta ∂Ω = {0,1}-en Dirichlet muturbalio homogeneoekin. Orduan, bere formulazio bariazional ahulak honako forma bilineal eta lineal hauek ditu: b(u, v) = Z1 0 u′v′, l(v) = v(1/2), u ∈U=H1 0(0,1) = V∋v, (5.7) non bere soluzio zehatza honako hau baita: u∗(x)=2x, 0 ≤x≤1/2 bada; eta u∗(x) = 2(1 −x), 1/2≤x≤1 bada. Gainera, Ritz motako minimizaziobirformulazio bat onartzen du, honela: u∗= arg min u∈U 1 2Z1 0 [u′(x)]2dx −u(1/2).(5.8) Orain, 2 tanh(θ(x−1/2))-k (u∗)′arbitrarioki hurbiltzeko gaitasuna duela ohartuz —gogoratu 2.4. atala—, analitikoki integratuz eta mutur-balioak ezarriz, hurrengo NN arkitektura (atipikoa) lortzen dugu uθ-k hurbilketa arbitrarioa izan dezan u∗-rekiko: xuθ(x) = 2 θnlog cosh(θ/2)−log cosh(θ(1/2−x))o.(5.9) Ondorioz, kalkulu errazen bidez honako hauek lortzen ditugu: F(θ)=2−4 θtanh(θ/2) −8 θlog cosh(θ/2),(5.10a) ∂F ∂θ (θ) = −4(θ−1) tanh(θ/2) −2θsech2(θ/2) + 8 log[cosh(θ/2)] θ2,(5.10b) θ= 0 denean. Konturatu konbergentzia/hurbiltasuna lortzen dela θ→ ∞ denean, hurrengo ondoriozkoekin: F(θ)→ −2 eta ∂F ∂θ (θ)→0. 5.3.2. Integrazioa Entrenamenduan integrazio-errorea murrizteko, (5.4)-ren ordez, honako errepikapenprozesu hau erabiliko dugu: F(θt)≈ Lt:= ®L(θ0), t = 0 bada, αtL(θt) + (1 −αt)Lt−1, t ≥1 bada,(5.11a) non L(θt) MCren integrazio zenbatespena baita t. iteraziorako, eta {αt}t≥0 aukeratutako koefizienteen segida bat izanik 0 < αt≤1 eta α0= 1 betetzen dituena. Modu hedatuan, honako itxura hau du: Lt= t X l=0 αl t−l Y s=1 (1 −αl+s)!L(θl).(5.11b) 104 5.3. Memorian oinarritutako integrazioa eta optimizazioa Horrek erakusten du LtF(θt)-ren zenbatespena oraingo eta aurreko MC integrazio zenbatespen guztien konbinazio lineala dela. Ohartu, t≥0 guztietarako αt= 1 bada, orduan memoriarik gabeko ohiko MC integrazio kasua berreskuratzen dugula. 5.1. irudiaren kurba gorriak erakusten du memorian oinarritutako Ltgalerafuntzioaren eboluzioa entrenamenduan zehar, SGD optimizazio arruntaren, (5.6)- ren, arabera egiten denean eta hurrengo koefizienteen hautaketarekin: αt= e−0.001t+ 0.001. Hasieran, errore handiekin integratzen dugu (αtia bat da, eta, beraz, ia ez dago memoriarik Lt-n). Hala ere, entrenamenduan aurrera egin ahala, memoria gehitzen diogu Lt-ri (αttxikituz), eta, ondorioz, integrazioerrorea murriztu egiten da. Lt-k hurbilketa zehatzagoak egiten ditu L(θt)-k baino, konbergentziaren jarraipen hobea egiteko aukera eskainiz, adibidez, entrenamenduan zehar gelditzeko irizpide egokiak ezartzekoa. 0 2,000 4,000 6,000 8,000 10,000 −2 −1.5 −1 t galera L(θt)LtF(θt)balio zehatz optimoa 5.1. irudia: 5.3.1. atalean aukeztutako NNaren entrenamendua, parametro entrenagarri bakarrekoari dagokiona, eta haren arkitekturak F zehazki kalkulatzeko aukera ematen duena. Entrenamendua (5.6)- ren arabera egiten da, eta, beraz, L(θt) da entrenamendurako galera. Lteta F(θt) bakarrik jarraipena egiteko kalkulatzen dira. 5.3.3. Optimizazioa Nahiz eta proposatutako eskemak F-ren hurbilketa hobetzen duen, erabiltzen den SGD eskema (5.11)-n SGD klasikoaren baliokidea da ikaskuntza-tasa αt-rekin biderkatzen denean. Hala ere, memorian oinarritutako integrazioaren eskema gradienteei eman diezaiekegu, g(θt) MC integrazioaren bidez lortzen delako —gogoratu (5.5)—. Beraz, 105 5. Memorian onarritutako Monte Carlo integrazioa gradienteentzako memoria-eskema honela geratuko litzateke: ∂F ∂θ (θt)≈gt:= ®g(θ0), t = 0 bada, γtg(θt) + (1 −γt)gt−1, t ≥1 bada,(5.12a) non {γt}t≥0aukeratutako koefizienteen segida baita 0 < γt≤1 eta γ0= 1 izanik. Orduan, memorian oinarritutako SGD optimizatzaile bat lortzen dugu, g(θt) terminoaren ordez gtterminoa erabiltzen duena —gogoratu (5.6)—; hau da, θt+1 := θt−λgt.(5.12b) Modu hedatuan, honako itxura hau du: gt= t X l=0 γl t−l Y s=1 (1 −γl+s)!g(θl).(5.12c) 5.2. irudiak gradientearen eboluzioa erakusten du 5.1. irudiko baldintzetan. Lehen bezala, gainbehera esponentziala erabili dugu, αt=γtaukeratuz t≥ 0 denean. gt-k gradiente zehatzen hurbilketa zehatzagoak ematen ditu g(θt)-k baino, portaera zaratatsua minimizatuz. 0 2,000 4,000 6,000 8,000 10,000 −0.6 −0.4 −0.2 0 t gradientea g(θt)gt∂F ∂θ (θt)balio zehatz optimoa 5.2. irudia: Gradientearen eboluzioa 5.1. irudiko entrenamenduaren baldintzetan. Optimizazioa (5.6)-ren arabera egin da g(θt) erabiliz. gteta ∂F ∂θ (θt) kalkulatu dira jarraipenerako. αteta γtkoefizienteen doikuntza egokia funtsezkoa da integrazioaren errendimendua maximizatzeko. Koefizienteak altuak (memoria baxua) izan behar dira inplikatutako integralak azkar aldatzen direnean (adibidez, entrenamenduaren hasieran). Aitzitik, hurbilpena orekatik gertu dagoenean eta inplikatutako integralak poliki-poliki aldatzen direnean, koefizienteak baxuak (memoria altua) izan behar dira. 106 5.4. Momentum metodoarekiko lotura 5.4. Momentum metodoarekiko lotura Momentum duen SGD optimizatzailea (SGDM)4[177,209] hurrengo bi urratseko metodo errekurtsiboaren arabera aurkezten da: vt+1 := βvt−g(θt),(5.13a) θt+1 := θt+λvt+1,(5.13b) non vtmomentum metagailu-tasa baita v0= 0 balioarekin hasiaraziz, eta 0 ≤ β < 1momentum-koefizientea izanik. β= 0 bada, SGD optimizatzaile arrunta berreskuratzen dugu (5.6). (5.13) erregela berridazten badugu θt+1 =θt−λgt eskemaren arabera, honako hau daukagu: gt=g(θt) + βgt−1= t X l=0 βt−lg(θl).(5.13c) SGDMaren bertsio sofistikatuago bat entrenamenduan zehar momentum-koefizientea aldatzean datza (ikus, adibidez, [50]); hau da, (5.13) araberako erregela jarraituz baina β-ren ordez βt∈[0,1) erabilita, {βt}t≥1hautatutako segidarako. Orduan, gtterminoak honako itxura hau du: gt=g(θt) + βtgt−1= t X l=0 t−l Y s=1 βl+s!g(θl).(5.14) Entrenamenduan zehar βthiperparametro egokiak hautatzea erronka bat da. Hala ere, SGDM optimizatzailearen ikaskuntza-tasa eta momentum-koefizientea berregokituz, γt-ren arabera t≥1 denean —gogoratu (5.12)—, honela: λt:= λt−1 γt γt−1 , λ0:= λ, (5.15a) βt:= γt−1 1−γt γt ,(5.15b) orduan memorian oinarritutako proposamena berreskuratzen dugu (5.12). Bi optimizazioek, (5.12)-k eta (5.13)–(5.14)-k, estokastikoki gradienteak metatzen dituzte parametro entrenagarriak berregokitzeko. Hala ere, (5.13)–(5.14)-k aurreko gradienteen batezbesteko geometrikoa hartzen du kontuan, gt-k ∂F ∂θ (θt)- ren antza hartzeko asmorik gabe. Gure proposamenak, (5.12)-k, ordea, egungo eta aurreko gradienteak berreskalatzen ditu iterazioz iterazio, gt-k ∂F ∂θ (θt) imita dezan —gogoratu 5.2. irudia—. 4Ingelesetik, SGD with Momentum (SGDM). 107 7. Lorpenak 7.2. Kongresuak Nazioko eta nazioarteko konferentzietan egindako ekarpenak. Aurkezlea agertzen da azpimarratuta. 2023 Carlos Uriarte, David Pardo, Ignacio Muga, Judit Mu˜noz-Matute. The Deep Double Ritz Method: a Deep Learning Residual Minimization Method for solving Partial Differential Equations. 2nd IACM Mechanistic Machine Learning and Digital Engineering for Computational Science and Technology Conference. University of Texas El Paso-n, AEB. 2023ko irailak 24-27. 2023 Carlos Uriarte, Jamie M. Taylor, David Pardo, Oscar A. Rodr´ıguez, Patrick Vega. Memory-based Monte Carlo integration for solving Partial Differential Equations using Neural Networks. XXXI COMCA – Congreso de Matem´atica Capricornio, Antofagasta, Txile. 2023ko abuztuak 2-4. 2023 Carlos Uriarte, Jamie M. Taylor, David Pardo, Oscar A. Rodr´ıguez, Patrick Vega. Memory-based Monte Carlo integration for solving Partial Differential Equations using Neural Networks. ICCS 2023 – 23rd International Conference on Computational Science. Praga, Txekiar Errepublika. 2023ko uztailak 3-5. 2023 Carlos Uriarte, David Pardo, Ignacio Muga. Goal-Oriented Deep Ritz and Least-Squares methods. ICCS 2023 – 23rd International Conference on Computational Science. Praga, Txekiar Errepublika. 2023ko uztailak 35. 2022 Carlos Uriarte, David Pardo, and ´ Angel J. Omella. A Finite Element based Deep Learning solver for parametric PDEs. SIAM MDS22 – Conference on Mathematics of Data Science. San Diego, Kalifornia, AEB. 2023ko irailak 26-30. 2022 Carlos Uriarte, David Pardo, Ignacio Muga, and Judit Mu˜noz-Matute. Solving Partial Differential Equations using Adversarial Neural Networks. CMN 2022 – Congress on Numerical Methods in Engineering. Las Palmas Kanaria Handikoa, Espainia. 2022ko irailak 12-14. 2022 Carlos Uriarte, David Pardo, Ignacio Muga, and Judit Mu˜noz-Matute. Adversarial Neural Networks for solving variationally formulated Partial Differential Equations. WCCM/APCOM 2022 – 15th World Congress on Computational Mechanics and 8th Asian Pacific Congress on Computational Mechanics. Yokohama, Japonia. 2022ko uztailak 31 – abuztuak 5. 114 7.3. Ikastaroak, mintegiak eta tailerrak 2022 Carlos Uriarte, David Pardo, Ignacio Muga, and Judit Mu˜noz-Matute. An Adversarial Networks approach for solving Partial Differential Equations. ICCS 2022 – 22nd International Conference on Computational Science. Londres, Erresuma Batua. 2022ko ekainak 21-23. 2022 Carlos Uriarte, David Pardo, Ignacio Muga, and Judit Mu˜noz-Matute. A Generative Adversarial Networks approach for solving Partial Differential Equations. ECCOMAS 2022 – 8th European Congress on Computational Methods in Applied Sciences and Engineering. Oslo, Noruega. 2022ko ekainak 5-9. 2021 Carlos Uriarte, ´ Angel J. Omella, David Pardo. A Finite Element based Deep Learning solver for parametric PDEs. ICCS 2021 – 21st International Conference on Computational Science. Krakovia, Polonia. 2021eko ekainak 16-18. 7.3. Ikastaroak, mintegiak eta tailerrak Ikastaro, mintegi eta tailerretan egindako ekarpenak. Azpimarratuta agertzen da ekarpen nagusia egin duen pertsona. Kasuren batean azpimarratutakorik ez egoteak esan nahi du ekarpen nagusia egin duen pertsona zehatzik ez dagoela. 2023 Carlos Uriarte, David Pardo. (Goal-Oriented) Deep Residual Minimization Methods. Mintegia Faculty of Computer Science, Electronics, and Telecommunications of the AGH University of Science and Technology-n, Krakovia, Polonia. 2023ko uztailak 11. 2023 Tom´as Teijeiro, ´ Angel J. Omella, Jamie M. Taylor, Carlos Uriarte, David Pardo. Parametric PDEs using Deep Learning. Ikastaroaren eta taldetailerraren antolatzailea Asturiasen, Espainia. 2023ko maiatzak 21-27. 2023 Carlos Uriarte, David Pardo, Ignacio Muga, Judit Mu˜noz-Matute. A Deep Double Ritz Method for solving Partial Differential Equations using Neural Networks. Ahozko aurkezpena Workshop Numerical Methods in Geophysics: Present, Future, and Applications-en, Valpara´ıso, Txile. 2023ko urtarrilak 12-13. 2023 ´ Angel J. Omella, Carlos Uriarte, and David Pardo. Coding Deep Neural Networks for PDEs. Irakaslea Pontificia Universidad Cat´olica de Valpara´ıso-n antolatutako eta Olmu´e-n garatutako ikastaroan, Txile. 2023ko urtarrilak 15-20. 115 7. Lorpenak 2022 Carlos Uriarte, David Pardo, Ignacio Muga, Judit Mu˜noz-Matute. A Deep Double Ritz Method for solving Partial Differential Equations. Ahozko aurkezpena XC Encuentro Anual de la Sociedad Matem´atica de Chile Punta de Tralca-n, Txile. 2022ko abenduak 8-10. 2022 Carlos Uriarte, David Pardo, Ignacio Muga, Judit Mu˜noz-Matute. A Deep Double Ritz Method for solving Partial Differential Equations. Mintegia Pontificia Universidad Cat´olica de Valpara´ıso-ko Matematikako Institutan, Txile. 2022ko azaroak 11. 2022 David Pardo, Magdalena Strugaru, Jamie M. Taylor, ´ Angel J. Omella, Jon A. Rivera, Carlos Uriarte, Ignacio Muga, Judit Mu˜noz-Matute. Deep Learning for Simulation and Inversion Problems. Hitzaldi nagusia Oden Institute for Computational Engineering and Sciences of the University of Texas-en, Austin, Texas, AEB. 2022ko urriak 21. 2022 David Pardo, ´ Angel J. Omella, Jamie M. Taylor, Carlos Uriarte, Jon A. Rivera, Magdalena Strugaru. Deep Learning for Simulation and Inversion Problems. Hitzaldi nagusia 9ºCongreso Metropolitano en Modelizaci´on y Simulaci´on Num´erica-n, Mexiko. 2022ko maiatzak 4-6. 2021 David Pardo, ´ Angel J. Omella, Jon Ander River, Carlos Uriarte, Ana Fern´adez-Navamuel. Solving Forward and Inverse Problems with Deep Learning. Ikastaroaren eta talde-tailerraren organizatzailea Kantabrian, Espainia. 2021ko maiatzak 10-18. 2021 Carlos Uriarte, David Pardo, ´ Angel J. Omella. A Finite Element based Deep Learning solver for parametric PDEs. Mintegia Pontificia Universidad Cat´olica de Valpara´ıso-ko Matematikako Institutuan. 2021ko martxoak 26. 2021 Carlos Uriarte, David Pardo, ´ Angel J. Omella. A Finite Element based Deep Learning solver for parametric PDEs. Mintegia Faculty of Computer Science, Electronics, and Telecommunications of the AGH University of Science and Technology-n, Krakovia, Polonia. 2021ko urtarrilak 21. 2020 David Pardo, ´ Angel J. Omella, Carlos Uriarte. Deep Learning for Solving Inverse Problems using TF2.0. Ikastaroaren eta talde-tailerraren antolatzailea Asturiasen, Espainia. 2020ko ekainak 28uztailak 4. 116 7.4. Ikerketa-egonaldiak 7.4. Ikerketa-egonaldiak Hogeita hamar egun baino gehiago iraun duten ikerketa-egonaldiak. 2023 Pontificia Universidad Cat´olica de Valpara´ıso-ko Matematikako Institutuan, Txile. Kudeatzaile-akademikoa: Ignacio Muga irakaslea (92 egun). 2022 Oden Institute for Computational Engineering and Sciences of the University of Texas at Austin-en, AEB. Kudeatzaile-akademikoa: Leszek F. Demkowicz irakaslea (37 egun). 7.5. Dibulgazio-jarduerak Dibulgazioko hitzaldiak. 2023 Elisabete Alberdi, Carlos Uriarte. Matematika eguneroko bizitzan. Ekuazio diferentzialak mundua azaltzeko. Dibulgazio hitzaldia Uhagon Kulturgunean, Markina-Xemein, Bizkaia, Espainia. 2023ko maiatzak 12. 2022 Carlos Uriarte. Ekuazio diferentzialak mundua azaltzeko eta sare-neuronal artifizialak horiek ebazteko. Dibulgazio hitzaldia Matematika Eguneroko Bizitzan 19. Edizioan, Bidebarrieta Kulturgunean, Bilbo, Espainia. 2022ko maiatzak 19. 117 Appendices/Eranskinak 119 A. Introduction (Chapter 1) A.1. Motivation Partial Differential Equations (PDEs)are of great value to society due to their broad applicability in modeling multiple biological, physical, or social phenomena [73,207,70,192]. Using derivatives in these equations allows us to describe complex relationships and their rates of change over time and/or space. For example, they allow to model: heat transfer [102,88], electromagnetic fields [106,221], fluid dynamics [23,14,223], population evolution [96,160,38], and financial [26,101] or health-care [193,141,165] forecasts. However, establishing an equation to describe a physical or social situation is only the first step. There arise different kinds of problems when modeling via PDEs, and the way we treat or approach them depends on the available computational resources. In forward problems, we are interested in determining the function that satisfies a PDE given some initial or boundary conditions. In this way, the solution to forward problems can provide insights into the modeled phenomena and can be used to make predictions, usually via post-processed measurements taken from the solution. In inverse problems, we have measurements of the solution, but we need to discover the value of specific parameters in the equation. This is an indirect problem, as we need to determine the parameters that provide a PDE solution that agrees with the given measurements [210]. Example A.1 illustrates a possible framework of forward and inverse problems in electromagnetic fields [72]. Example A.1 (Electromagnetic fields).Maxwell’s equations model electromagnetic fields according to four well-known physical laws [72,106]:            ∇×E=−jωµH−M,Faraday’s Law, ∇×H= (σ+jωε)E+J,Amper`e’s Law, ∇·(εE) = ρf,Gauss’ Law of Electricity, ∇·(µH)=0,Gauss’ Law of Magnetism. (A.1) Here, Eand Hdenote the electric and magnetic fields in the frequency domain, respectively; Jand Mare given source terms; σstands for the electrical conduc121 A. Introduction (Chapter 1) tivity of the media, εfor the electrical permittivity, µfor the magnetic permeability, ρffor the electric charge density, jfor the imaginary unit, and ωfor the angular frequency. Then, a possible forward-inverse scenario on electromagnetic fields could consist of parameters {σ,µ,ε}and measurements {Z}, where Z denotes the impedance tensor defined by E=ZH. For further details, see, e.g., [9]. Figure A.1 shows a graphic summary of this modeled framework. Parameters {σ,µ,ε} Solution {E,H} Measurements {Z} Solve PDE Post-processing Forward problem Inverse problem Figure A.1: Forward and inverse problems sketch in electromagnetic fields. Solving forward problems requires, in general, fewer resources than inverse problems. Given a choice of parameters, a forward problem typically consists of performing a single PDE simulation. In contrast, an inverse problem usually requires the iterative resolution of forward problems, conveniently readjusting the choice of parameters at each iteration based on the outcome of previous forward simulations [210,103]. Connecting both types of problems, we have parametric problems, where the aim is to study the behavior of the solution/measurements as a function of the parameters. The solution can depend on one or several parameters, and we can study how changes in these parameters affect the behavior of the system. Alternatively, the parametric approach can be viewed as a forward problem where the parameters of the equation take the role of variables of the solution. Typically, the parametric forward problem is nonlinear with respect to the introduced new variables, enormously escalating and hindering its resolvability. In this dissertation, we will only address the task of solving parametric and nonparametric linear PDEs, which is critical in many applications and frameworks. We will indistinctly call PDEs both one-dimensional—also known as Ordinary Differential Equations (ODEs)—and higher-dimensional differential equations. 122 A.2. Traditional Numerical Methods A.2. Traditional Numerical Methods There are various methods for solving PDEs, which can be broadly classified into two categories: analytic and numerical methods. Analytic methods involve using mathematical techniques such as separation of variables, integral transforms, or complex analysis to find an exact solution to the PDE [205,87,213]. However, these methods are often limited to relatively simple PDEs and geometries. Additionally, even if the PDE is analytically solvable, the resulting solution may be expressed in a complicated form (e.g., as an infinite power series), which can be difficult to interpret and apply in practical scenarios. On the other hand, numerical methods involve approximating the solution to the PDE using a computational algorithm. Generally, these methods propose approximations that are easily interpretable or numerically recoverable (e.g., as a finite linear combination of simple prescribed functions). Below, we briefly review the main aspects of three widely known and employed kinds of numerical methods for solving PDEs. The Finite Difference Method (FDM)discretizes the spatial domain and/or time interval of the PDE into a finite number of subdomains [203,178,130]. Then, the evaluation of the derivatives at these discrete points is approximated by solving algebraic equations containing finite differences and values from nearby points. The FDM converts a PDE into a system of linear equations that can be solved by matrix algebra techniques. It is conceptually simple but challenging in design for complex problems or irregular domains. In particular, it suffers from the so-called curse of dimensionality [29,30], where the size of the involved matrix grows exponentially with the dimension of the problem. Similarly, the Finite Element Method (FEM)[100,39,138,184] proposes a mesh-based scheme where the approximated solution is in the form of a finite linear combination of some prescribed functions (typically, piecewise polynomial) with local support. The support is purposely designed on specific regions of the domain, called elements, so the resulting matrix is sparse. The FEM is particularly suitable for dealing with PDEs in variational form but suffers from similar difficulties as the FDM. Spectral methods [81,36,46] are a class of numerical methods that approximate functions using a linear combination of orthogonal basis functions with global support, such as Chebyshev polynomials or Fourier series. These methods can achieve very high accuracy, and are particularly useful for smooth solutions with high oscillatory behavior, but may become computationally unfeasible for complex geometries. In essence, the numerical methods mentioned above are based on pre-establishing a finite-dimensional basis and parameterizing the approximation via the coefficients of the corresponding linear combination. When using a mesh-based basis 123