Solving Differential Equations with Constrained Learning
Authors: Viggo Moro
Luiz F. O. Chamon
Venue: ICLR 2025 Poster (2025)
Keywords: Constrained learning
partial differential equations
neural operators
physics-informed neural networks
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: (Partial) differential equations (PDEs) are fundamental tools for describing natural phenomena, making their solution crucial in science and engineering. While traditional methods, such as the finite element method, provide reliable solutions, their accuracy is often tied to the use of computationally intensive fine meshes. Moreover, they do not naturally account for measurements or prior solutions, and any change in the problem parameters requires results to be fully recomputed. Neural network-based approaches, such as physics-informed neural networks and neural operators, offer a mesh-free alternative by directly fitting those models to the PDE solution. They can also integrate prior knowledge and tackle entire families of PDEs by simply aggregating additional training losses. Nevertheless, they are highly sensitive to hyperparameters such as collocation points and the weights associated with each loss. This paper addresses these challenges by developing a science-constrained learning (SCL) framework. It demonstrates that finding a (weak) solution of a PDE is equivalent to solving a constrained learning problem with worst-case losses. This explains the limitations of previous methods that minimize the expected value of aggregated losses. SCL also organically integrates structural constraints (e.g., invariances) and (partial) measurements or known solutions. The resulting constrained learning problems can be tackled using a practical algorithm that yields accurate solutions across a variety of PDEs, neural network architectures, and prior knowledge levels without extensive hyperparameter tuning and sometimes even at a lower computational cost.

Published as a conference paper at ICLR 2025
SOLVING DIFFERENTIAL EQUATIONS WITH CONSTRAINED LEARNING
Viggo Moro
Luiz F. O. Chamon
University of Oxford
École polytechnique luiz. chanonn@oolyyeecmiiiuu.....
vi moromms www.sso.....o aw..cuukkk...c.......
ABSTRACT
(Partial) differential equations (PDEs) fundamental tools for describing natu- ral phenomena, making their solution crucial in science and engineering. While traditional methods, such the finite element method, provide reliable solutions, their accuracy is often tied to the use of computationally intensive fine meshes. Moreover, they do not naturally account for measurements prior solutions, and any change in the problem parameters requires results to be fully recomputed. Neural network-based approaches, such as physics-informed neural networks and neural operators, offer mesh-free alternative by directly fitting those models to the PDE solution. They can also integrate prior knowledge and tackle entire families of PDEs by simply aggregating additional training losses. Nevertheless, they are highly sensitive to hyperparameters such as collocation points and the weights associated with each loss. This paper addresses these challenges by developing a
learning (SCL) framework. demonstrates that finding a (weak) solution of PDE is equivalent to solving constrained learning problem with worst-case losses. This explains the limitations of previous methods that minimize the expected value of aggregated losses. SCL also organically integrates structural constraints (e.g., invariances) and (partial) measurements or known solu- tions. The resulting constrained learning problems can be tackled using practical algorithm that yields accurate solutions across variety of PDEs, neural network architectures, and prior knowledge levels without extensive hyperparameter tuning and sometimes even at a lower computational cost.
INTRODUCTION
(Partial) differential equations (PDEs) are key tools in science and engineering, playing a central role in the solution of inverse problems, systems engineering, and the description of natural phenom- ena (Lustig et al. 2008 Potter et al. 2010, Molesky et al., 2018 8 Evans 2010). As such, a variety of numerical methods have been developed to approximate their solutions, such as the well-known finite element method (FEM). Despite their celebrated precision and approximation guarantees, these methods provide solutions dition to mesh size, requires the solution to be recomputed. They are therefore unable to incorporate prior knowledge, such as real-world measurements known solutions similar equations (Brenner & Scott 2007 Vequee
single PDE at time. Any change to the problem, from boundary con-
2007 Katsikadelis 2016)..
Methods based on neural networks (NNs), such as physies-informed NNs (PINNs) ( (agarii al. 1998 Raissi et al.) 2019 Lu et al. 2021b) and neural operators (NOs) (Li et al. 2021 L et al. 2
Rahman et al.) 2023), have been developed with these challenges in mind. Rather than discretizing the PDE, they directly fit a NN to its solution. They can therefore be trained to simultancously solve entire families of PDEs and interpolate known solutions by simply incorporating additional losses
to their training objectives (Li et al., 2021] Cho et al. 2024| Li et al.) 2024). Yet, these methods are highly sensitive to hyperparameters such as the weights used to combine the training losses and the collocation points used to evaluate the PDE residuals, which often leads to low quality or
trivial solutions (Krishnapriyan et al. 2021 Wight & Zhao 2021 Wang et al.) 2022b a). This has prompted variety of heuristics to be proposed based on ad hoc weight updates WWagg al. 2021a
Maddu et al.) 2022) and adaptive or causal sampling Nabian et al. (202I); Krishnapriyan et al. (202))
McClenny & Braga-Neto (2023): Penwarden et al. (023) (see Appendix G for further related works).

---

Published as a conference paper at ICLR 2025
This paper shows that these limitations not methodological, but epistemological. It is not issue of how the problem is solved, but which problem is being solved. To do so, it
proves that obtaining (weak) solution of PDE requires solving constrained learning problem with worst-case losses (Prop. 3.Ip. i.c., is not enough to use either constrained formulations
et al. 2021b Basir & Senocak 2022) worst-case losses (Wang al., 2022a Daw al. 2023
incorporates prior scientific knowledge on the structure (e.g., invariance) and value (c.g., measure- ments, set points) of the solution without resorting to specialized models or data transforms (Sec. We therefore dub this approach science-constrained learning (SCL);
develops a practical algorithm that foregoes the careful selection of loss weights and collocation points. Contrary other methods, it explicitly approximates the (weak) solution of PDEs and yields reliability metrics that capture the difficulty of fitting specific PDE parameters and/or data points (Sec.
illustrates the effectiveness (accuracy and computational cost) of SCL (Sec. 5 for a diverse set of PDEs, NN architectures (e.g., MLPs and NOs), and problem types (solving single or parametric family of PDEs; interpolating known solutions; identifying settings that are difficult to fit).
2
PROBLEM FORMULATION
2.1 BOUNDARY VALUE PROBLEMS
Consider (bounded, connected, open) region Ω Rdd with (smooth) boundary 8Ω and a (partial) dif- ferential operator D. with coefficients т E Il C Ra defined on the domain D = Ω x (0, T'. Here, π captures (finite) set of parameters of the phenomenon, such as the diffusion rate or viscosity. Given a space F of functions mapping D to R, we define a boundary value problem (B VP) as find such that
u J
T(æ t), u(x, 0) = h(z, 0), u(x,1) h(x, l),
(z, [u[u[(z,
(a, t) E E x (0,7)
(PDE)
(BVP)
x x
(IC)
(x, E 0n x (0,
(BC)
where : D - R is a forcing function and h: B B describes the boundary (BC) and initial (IC) conditions over B (? x ( x (1) for f2 Quan. In what follows, we always consider the initial condition part of the boundary conditions and refer to them jointly as BC. Note that the developments in this paper also apply to formulations of (BVP) involving other BCs (e.g.. Neumann, periodic), parameterized BCs, and vector-valued PDEs. We showcase a variety of phenomena that can be described by BVP) in App. A and refer to, e.g., Evann 2010 for a more detailed treatment.
In general, a (strong) solution of (BVP) need not exist in any function space F. This motivates the rise of relaxations such as the weak formulation that replaces the pointwise equation (PDE) in (BVP) by the integral equation"
[ Dulu[(z, ()p( 1) dxdt = -
T(x, t)yp(x, t) dædt, for all yp PTT.
(I)
The are known as test functions and T is typically taken to be Sobolev space due to its natural compatibility with this setting (see App. from the fact that a solution of (BVP) (when it exists) is also solution of (D. although the converse not necessarily true. Indeed, [ allows for a wider range of solutions with less stringent regularity requirements, particularly with respect to continuity and differentiability (Evans) 2010). The BCs of (BV)) can often homogeneized (i.c., h 0) and imposed implicitly through T. thus fully describing its weak solution by (I) (Brenner & Scott 2007 Le VeVuue 2007 Katsikadelis 2016).
further details). The name weak formulation comes
2.2 SOLVING BOUNDARY VALUE PROBLEMS
There exists a wide range of numerical methods for solving BVPs, most of which rely on dis- cretizing (BVP) (e.g., finite difference method, FDM) its weak formulation (e.g., FEM). Their
Typicallly the weak formulation only considers the spatial domain x. handling t separately using time- stepping methods (Thomee 2013). Still, (T) is not uncommon and can be used to obtain weak formulations for parabolic PDEs (see, e.g., Knabner & Angerman 2003 Evans 2010 Steinbach & 7ank 2020).
2

---

Published as a conference paper at ICLR 2025
well-established approximation guarantees, accuracy, and stable implementations make them ubiqui- tous in scientific and engineering applications (Brenner & Scott 2007; Lequee 2007; Katsikadelis 2016). These classical methods, however, only tackle one VP at time and do not naturally incorporate prior knowledge, such as measurements or known (partial) solutions. These challenges can be addressed by directly parameterizing the solution of BVPs, most notably using NNs. While this approach may not achieve the precision of classical methods (Krishnapriyan et al. 2021] Wight & Zhao 202T] Grossmann al. 2024 McGreivy & Hakim 2024), they are able to simultaneously provide solutions for entire families of BVPs and extrapolate new solutions from existing ones. Generally speaking, these methods can be separated into unsupervised, that seek to solve (BVP) directly, and supervised, that leverage previously computed or measured solutions.
Unsupervised methods. These approaches train model ug D 11 R with parameters 0 E @ C RP (e.g., a multilayer perceptron, MLP) to fit (BVP). This is the case, for instance, of physics- informed neural networks (PINNs) that train up by solving, for fixed weights Upde: Mbc 0.
minimize Mpdee pxpee + (bc BEB
(PI)
M
(uo (1i) Emm T (z Em)) 4m)))
(0) -
(Im, E DD n. II,
M
m=1
N
(0) A
( (n) - M6 kn))
(Inz (n) E B.
N
n=1
We write ug (ar, () to emphasize that evaluate ug u which approximates the solution of (BVP) with coefficients π, at (a, E D. In practice, the input of the model ug is simply (I,t, x) E D x II. The losses fpde and foe promote the requirements in (PDE) and (IC)-(BC from (BVP). respectiely although the former is computed using automatic differentiation rather than discretization. Though the majority of PINNS target a single BVP, i.e., == 1 in PT), their extension to parameterized families of BVPs has been explored (Cho et al. 2024
Supervised methods. Rather than directly solving the BVP, these approaches fit the model us to a set of (partial) solutions uf, of (BVP). In this setting, ug is often a neural operator (NO) capable of handling infiite-dimensiona. (functional) inputs and outputs, such as forcing functions and ICs h(z, 0). Given, e.g., forcing-solutiong pairs (T,, u/). these NOs are trained by solving
minimize 6E
Σ
(T)(T)) u/] (
(PII)
In practice, the functions Tn and u. are discretized (in the time or spectral domain) to enable computations (Li et al. 202T Lu et al. 2021a Hao et al. 2023) Wei Zhang 2023). While not uncommon to combine PI and (PII), these semi- cemi-supereised methods typically rely MLPs Raissi ct al. 2019 Lu et al. 2021b), since it can be challenging to evaluate D(ue] for NOs (Li et al. 2024).
Limitations. Though effective in many applications, these methods are very sensitive to their hyper- parameters. Indeed, the choice of collocation points (xr, (), PDE coefficients n, and weights ('pdes /4h
affect both the quality and computational complexity of (PI) (Krishnapryan et al. 2021 Wightt &
Zhao 202T; Wang al. [22 la). The same holds for the discretization of the objective in (PI)) (Li
et al.] 2021 Hao et al.) 2023] Wei & Zhang 2023). Supervised methods face the additional challenge
that acquiring the PDE solutions } used in (PII) can be expensive (relying on, e.g., classical
methods) and it is challenging obtain good performance from small datasets. This issue is aggra- vated by the heterogeneous difficulty of fitting each solution. A variety of heuristics for collocation points (Nabian et al. (202T] Daw et al. 2023 Penwarden et al. (023)., weights (Wang et al.) 202 la Maddu et al., 2022 McClenny & Braga-Neto 2023), and PDE solutions (Pestouric et al. 2023
Musekamp al.] (2024) have been put forward to mitigate these challenges. Yet, they generally focus on specific "failure modes" or BVPs and seldom address the non-trivial interactions of these yperparameters, limiting their effectiveness.
3
LEARNING
In this section, we argue that the challenges faced by previous NN-based BVP solvers arise not because of how (PT] and (PII) are solved, but because they are not the appropriate problems solve in
3

---

Published as a conference paper at ICLR 2025
the first place. To do so, show that obtaining (weak) solution of (BVP) is equivalent to solving constrained learning problem with worst-case losses. Hence, it is not enough to use (approximations
of) worst-case losses in, e.g., (Wang et al. 2022a Daw et al., 2023), adapting loss weights as in, e.g., (Wang et al. 2021a Lu et al., 2021b, McClenny & Braga-Nert Neto 2023). Building on this result, we show how incorporate other forms of scientific knowledge that not mechanistic (i.e., PDEs) without resorting to specialized models, including structural information (e.g., invariances) and observations (measurements, simulations) of the solution. The resulting science-constrained learning (SCL) problem accommodates a variety of knowledge settings, from unsupervised to supervised, and is amenable to a practical algorithm capable of effectively tackling entire families of BVPs and extrapolating solutions from existing ones (Sec. 4 and 5
In the remainder of this paper, we use up to refer to any parameterized model (MIP, NO, etc.). For clarity, we derive our results for a single BVP instance, omitting the dependence on and/or т. We consider these extensions at the end of the section.
Mechanistic (PDE) knowledge. We begin by showing how weak solutions of (BVP) can be obtained using constrained learning. To do so, we relax the BCs of (BVP) to relate the weak formulation 1 a distributionally robust constraint. The BCs are then reintroduced using a constrained formulation. We start with the following proposition, where W refers to the (k, p)-th order Sobolev space (see App. B and (S) denotes the space of square-integrable probability distributions supported on S.
Proposition 3.1. Let ut E wk' (DD) where 2 1 is the degree of the differential operator D, be such that sup (x,,t)) (Dlu (xx, T(z,t)) If the dimension d of Q satisfies d s
4k' - 1. then ut satisfies ! with T wk (DD).
proof is provided in Appendix B The equality constraint suggested by Prop. 3.1 enforces ① but does not impose the BCs of (BVP). Since they must hold for all (Ir t) E B, it is more appropriate to incorporate them using a worst-case loss, namely, (a.(a. ( (ug (r, (u
h(x, t)) rather than an average loss in PI As long (I. - (ug(x, t) - h(zz, +)² is function in L this is equivalent to a distributionally robust loss similar to the one from Prop. 3.1 (see Appendix
B). We therefore
conclude that weak solution of BVP) is obtained by solving
minimize PEB
sup (BP)
(a.t)~ (uo(e,) (xx,)))
(PIII)
subject to
sup (EP)(D)
(D(uol(z, -(x,) ((x,t)) Sc.
where e 2 0 controls the trade-off between fitting the PDE and the BCs when us is not expressive enough to satisfy both.
Prop. 3.1 clucidates the challenges arising from the choice of collocation points in the unsupervised
approach most notably PINNs. Indeed, it is not enough to use a fixed distribution (e.g.. uniform): satisfying requires training against all distributions & E pp? The use of worst-case losses in constrained formulation what makes (PIII) considerably different from previous adaptive sampling methods and loss-weighting schemes. In fact, contrary to previous approaches, by (approximately) solving (IIl (as detailed in Sec. 4). we indeed (approximately) solve (BPP). At the same time, Prop. 3.1 establishes limitation of lcarning-based solvers by restricting the smoothness of their solutions (essentially, solutions in ((d+))44.2 This can be an issue for large-scale dynamical systems, such as those found smart grid applications, or when transforming higher-order PDEs in higher-dimensional first-order systems. While Prop. 3.. describes a sufficient condition for u' to satisfy (1). a necessary condition can be obtained by restricting pd to sufficiently smooth distributions, namely, those belonging to Sobolev space.
Structural knowledge. The constrained form of (PII) suggests that other information can be incorporated as long as they can be formulated as learning objectives, i.e., statistical losses. This is the case of certain forms of structural knowledge. Indeed, it is often possible to obtain information about the structure of the solution of a BVP, such as invariances or symmetrics, without explicitly solving it Olver 1979 Akhound-Sadegh al. 2023). While this structural information is already
4

---

Published as a conference paper at ICLR 2025
encoded in (PI.I.. using it explicitly can reduce training time well as the number of both collocation
points and/or observations [as (PII)) needed. It also helps ensure the physical validity of outcomes by explicitly avoiding degenerate solutions, common failure mode of unsupervised methods such as PINNs IKrishnapriyan et al., 2021) (we do not observe such issues with (PIII), see Sec. (5.1). Structural knowledge can also be used remove solution ambiguities, e.g., for the eikonal equation that is invariant to the sign of the solution (see App. A)
While structural knowledge can sometimes be incorporated into the model ue, e.g., using equivariant architectures (Cohen & Welling 2016, Batzner et al., 2022), it can also be imposed as a worst- case constraint. This is convenient for when such models are intricate
design. Consider, for example, (finite) invariance group G whose elements act on the domain (x, E D such that u (xx 1) ut [r ()), where ut a solution of (BVP). This invariance can be enforced by
sup PEP)(D)
B(xe))=a) (uo(z, uolm(e,,t)) Se, 7 EEg.
(2)
Notice that we use the same distributionally robust formulation as for the BCs in (PIII). Similar constraints can be constructed for other structures, such as equivariance. In contrast to (PI.). where
we want e 0, it can be beneficial use a larger values in (2) to accommodate models Ug that cannot fully capture the solution invariances.
Observational knowledge. In addition to mechanistic (i.c., PDEs) and structural knowledge, also consider (partial, noisy) observations of the BVP solution, obtained either via classical methods (e.g., FEM) or real-world measurements. Although this type of information is commonly associated with NOs, seen as they are typically trained in supervised manner as in (PII,, is not limited to that architecture. Given observations w/,, I....,J, of solutions of (BVP. we may formulate constraints of the kind
E(z,))=m ( Se. j= j=
(3)
where m some distribution (typically uniform) of points on D. Note that (3) simply a different
way of writing the Lo2-norm from the objective of (PII). However, rather than averaging I,2 losses, (3) constraints the maximum error across data points. By considering each sample individually, accounts for the heterogeneous difficulty of fitting them and enables the tolerance to be adjusted individually for each observation, e.g., using larger values for noisier samples.
Science-constrained learning. Combining (PII)P with 2 and general SCL problem accommodating all knowledge sources considered so far. Explicitly,
able formulate
me (uo(x,,) - 1xx, []]] (EP2(B)
minimize 868
sup
(M)
subject to
E (m.t)~u (EP {
sup
7(,, })
I Cyúe
(M)
(SCL)
sup (EP)(D)
E (m..)r (10(2, u0 t))) S c,, 7 Eg
(S)
E (z,t)) (uo(z, - u,(x, co. co. j=
(0)
Note that any subset of the constraints in (SCL) be used depending the available information. The objective of (SCL) is also not restricted to the BCs and can be replaced by any of the other terms.
In fact. (SCL) can be formulated without objective, i.e., feasibility problem.
is straightforward to extend (SCL) to simultancously solve a parameterized family of BVPs. However, rather than discretizing the parameter space in (PI), rely worst-case formulation
5

---

Published as a conference paper at ICLR 2025
that considers all of its possible values rather than only finite subset. Explicitly, rewrite SSCL) as
minimize MEB
sup B(x,tm)s,, rmpp
(uo(z,))(<)) -0x,,),)))))))) [))
(M)
PEP2
subject to sup (x,,,))),,)) PEPa
- T(x, tc) $ Epdo
(M)
sup ((atm))))rrp Tmp ( (1)( (2 t) - uo 7)x(()))) S
(S)
PEP2 B(x,))/m
S en. j=l....,
(0)
(SCL')
where m. p are fixed distributions (e.g., uniform), un solution of (BVP) with coefficients π and forcing function Tj. and the invariance "/ is now parametrized to account for the fact that its action may depend on π (e-g., translation invariance with different strides). Note that the to are now supported on B x II or D x II, which we omit in (SCL') for clarity. Hence, they target not only BVss (parameters л) that are hard to fit, but also the regions of the domain responsible for this difficulty, enabling performances that would require fine discretizations (see Sec. 5 Yet, this approach is not directly applicable to diment-aadmeniionll parameters (e.g.. 7) as it requires sampling from a function space. We leave this extension for future work, considering here fixed distribution p.
In the next section, we develop a practical algorithm to tackle (SCL) and SCL' by (i) leveraging non-convex duality results from constrained learning (Chamon & Ribeiro 2020) Chamon et al. 2023)
and (ii) deriving explicit approximations the suprema over &b from which we can sample efficintll..
4 ALGORITHM
To develop a practical algorithm for SCL [and (SCL')), we need to overcome the fact that is (i) a non-convex constrained optimization problem involving (ii) worst-case losses. A typical approach to (i) is to combine the all losses as penalties into a single training objective as in (PI). Though penalties and constraints
essentially equivalent in convex optimization (strong duality, (Bertsekas) 2009),, this is not the case in the non-convex setting of (SCL). Hence, regardless of how the weights
in (PI) are adapted (e.g., WWang et al. (2021a); Wight & Zhao (202)): Lu al. (2021b): Basir &
Senocak (2022), it need not provide a solution of (SCL).
We overcome this issue by first tackling (ii) using the following proposition:
Proposition 4.1. Let - f(=) E T.2. Then, for all 5 >> O there exists < sup, {(2) such that supsepa Ez [((2)] < Ezwo. [[(=)] +6, where p2 e (2(2)) x [((=)-a), [a]+ max(0, a).
A proof based on Robey* et al. (202T) can be found in App. B Prop. 4.1 shows that the worst-case losses in /SCL)(SCL') can be approximated arbitrarily well by an expectation with respect to ela. a distribution proportional to a truncation of the underlying loss. For clarity, consider (SCL) with only constraint (M), but similar manipulations hold for the constraints (S) and (O) as well as SCL') Explicitly, (M) can be written
minimize 9EB
A(x, h(2,t))² { ((e,)))a (uo(z, t) (hooz,)))((,,))))
(PIV)
E B( (Dlus)(o, 1) - T(x, 1)) [] Se.
subject
for whe (xx [(ua(x,t) - ((z, a] and widee (x,t) x [(Dluo((z,)) r(a, a]
supported on B and D respectively.
Observe that (PIV now has the form of a constrained learning problem. i.c., constrained optimiza- tion problem with statistical losses. We can therefore use non-convex duality results from (Chamon & Ribeiro 2020 Chamon et al., 2023 Elenter et al. 2024) to show that, under typical conditions from (unconstrained) learning theory and for rich enough parametrization, its solution can be approximated by solving the empirical dual problem
maximize A
min f A),
(DIV)
Bee

---

Published as a conference paper at ICLR 2025
where L is the empirical Lagrangian of (PIV) based on samples (a xhh) the and (ahht t viale namely.
Nic
L(0, A)
(uo(zn be (2))
1e(1r (2))
Nbc
nm1
Male
Dauolaahh (nrte) - (nn (nde)
N ==1
Contrary to previous approaches based on (PI. such as (Wang et al.
2021a Wight & Zhao 2021
Daw al. 2023), (DIV) truly approximates the solution of (B Indeed, (Chamon al. 2023
Thm. 1) and (Elenter et al. 2024] Thm. 3.1) guarantee that solutions of (DIV) are near-optimal
and near-feasible for (PIV) and, in view of Prop. [.1 and (4.1). (BV)) (see App. D for details),
is worth noting that this is only possible because (PIV) is statistical problem. Although similar
Lagrangian formulations have been used Lu et al. (202 (b); Basir & Senocak (202),, they deal with
deterministic constrained problem (fixed collocation points) for which this duality does not hold.
From a practical perspective, (DIV) does not require extensive hyperparameter tuning [such as μ in J. seen À is optimization variable. What is more, despite non-convexity, the duality between (PIV) and (DIV) allows the solution X" of (DIV) to be interpreted as a sensitivity of the objective (BC residuals) to small relaxations of (Chamon & Ribeiro 2020 Chamon al. 2023 Hounie et al. 2023b). This information can be used to evaluate the fit of noisy measurements
in (O) the reliability of solutions for different parameters π in (SCL') (see Sec. Finally,
DIV) is amenable to practical algorithms such as dual ascent, which updates λo = 0 as
Npte
Ax+1 = An +
(Dluo (aahh (odo) - (x:ble (22e))
for of E argmin Ż(0, Ax). Aee
pdc
(5)
Even if the empirical Lagrangian minimizer of is only computed approximately, (5) can be shown to
converge to a neighborhood of solution of (DIV) (Chamon et al. 2023 Elenter et al. 2024).
Still, to turn (5) into a practical algorithm, we need to obtain samples from tha which is only known implicitly [see (PTV)) Nevertheless, since the losses in (PIV) are non-negative, the vi, are smooth, dully-supported be sampled using Markov Chain Monte Carlo (MCMC) techniques (Robert Casella (2004). We use the Metropolis-Hasting (MH) algorithm in our experiments since it avoids additional backward passes and higher-order derivatives resulting from differentiating D2 while still providing strong empirical results (see App-[Cand Sec.5) Exploring first-order methods (e.g.. Langevin Monte Carlo) and algorithms adapted to discontinuous distributions (c.g.. (Nishimura et al. 2020 to enable faster convergence (reduce N) and better approximations (increase a) is left for future work.
eqquare-itterrale distributions for a = 0. They are therefore amenable to
possible to replace v/a by fixed distribution (e.g., uniform) or even fixed points for some of the constraints in (SCL) at negligible accuracy costs. We find this is consistently the case for the BC objective, although we emphasize that this is architectures, such FNOs (Li et al. 2021), cannot make predictions at arbitrary points of the domain. It is then more appropriate to take Va to be a uniform distribution over equispaced points.
Additionally, certain
The resulting method for solving (SCL' summarized in Alg. [ Note that rather than obtaining Lagrangian minimizers as in (5). 1 alternates between optimizing for 0 (step 8) and À (step 9). Such gradient descent-ascent schemes are commonly used in convex optimization (Arrow et al. 1958
Bertsekas 2015), although their convergence is less well understood in non-convex settings Lin et al. (2020); Fiez et al. (202I): Yang al. (2022). The convergence guarantees for (5) can be recovered by repeating step 8 (6-update) multiple times before updating A. Note that Alg. [ does not rely on training heuristics, such as adaptive or causal sampling Krishnapriyan et al. (202T) b); McCleeny & Braga-Neto (2023): Daw et al. (2023); Penwarden et al. (2023): Wang et al. (2024). ad hoc weight updates (Wang et al. 2021a; Maddu et al., 2022), and conditional updates (Lu et al. 2021b) Basir & Senocak 2022). Indeed, steps 4-7 are empirical estimates of expectations with respect to vo that themselves approximate the worst-case losses in SCL' (Prop. 4.1.. Steps 8-9 describe a traditional
Fiez
primal-dual (gradient descent-ascent algorithm for solving problems such as (DIV). which itself
yields approximate solutions of (SCL') ) (Chamon & Ribeiro 2020 Chamon et al. 2023 Elenter et al.
2024) and, consequently, (BVP) (Prop. 3.1
7

---

Published as a conference paper at ICLR 2025
Algorithm 1 Primal-dual method for (SCL)
1: Inputs: Differential operator DR. invariant transformations E G observations set (Ty, Tyt uz)) parameterized model ud,, and xp = Xo = =100 =
2: for k 1.
K
Nuc
3:
1
Ub (mTh (x (hr) - h( ent) )
tbc (xhee in mhh
v0
NBC
n=1
Nuk
4:
pode k
D a/ [uo (TT (oihle (alle)
I [xp (pplee )
pde
yndo pile) ~ 2/ PDB
Npc
nml
N.
5:
un (1 (1) uo, (1
(14 1h mp)) 1/0
N
1
6:
uoe (Ty, Ty (x (2))) u (1 t
(xx (. m
No
7: 0k+1 = 0k
Up Volk + p Vol + Σ No Vole + = x Volk
i=1
j=1
8:
A = [2pe +
- -pde)] t
1a+1 = Ni + nnall
- -6)] 12t21 = [2 + nnllfh -co]
9:
end
EXPERIMENTS
In this section, we showcase the use of SCL by training MLPs and FNOs (Li et al. 2021 to solve six PDEs (convection, reaction-diffusion, eikonal, Burgers', diffusion-sorption, and Navier- Stokes). We consider different subsets of constraints from /SCL) illustrate a variety of knowledge settings, but train only the most suitable model in each case since our goal is to illustrate the natural uses of SCL rather than exhaust potential. Detailed descriptions are provided in the appendices, including BVPs (App.(A), training procedures (App. Ep. and further results (App. F. Code to reproduce these experiments is available at thttps: : //g1 thub. . com vmoro sequel, we use fixed points (z, for the BC objective rather than wit to illustrate how computational complexity can reduced without significantly affecting the results. We still use v for π.
In the
5.1 SOLVING A SPECIFIC BVP
We begin by solving (VP) for fixed parameters æ, forcing function r, and BCs. Though this may not be the best application for NN-based solvers [see, e.g., McGreivy & Hakim (2024)), it remains valid demonstration. We use (M) train MLPs to solve convection, reaction-diffusion. and eikonal problems, comparing the results to PINNs fi.e., (PI with weights chosen as in (Daw al. [202)]]. All experiments use 1000 collocation points per epoch obtained by (PINN) sampling uniformly (we find this performs better than using fixed points in, e.g., (Raissi et al.) 2019 Lu et al. 202Tb)), (R3) using the adaptive heuristic from (Daw et al. 2023), and (SCL) using MH after 4000 burn-in steps. Table shows that SCL matches and often outperforms other methods, particularly
Table Relative L2 error for solving specific BVPs (average standard deviation across seeds).
PINN PI 1.17 + 0.65 %
R3 Daw et al. 2023
SSL (M)
A = 30
0.999 0.53 %
0.971 0.30 %
Convection
B 50
56.5 20 %
29.0 = 33 %
3.74 = 0.87 %
(v. p) (3,3)
0.745 = 0.014 %
0.736 = 0.068 %
1.82 = 0.74
React.-Diff.
(v. (3,5) 5)
79.6 = 0.27 %
0.665 = 0.046 %
0.762 + 0.11%
Eikonal
87.1 = 39 %
85.5 = 27 %
9.95 + 1.9%
8

---

Published as a conference paper at ICLR 2025
(PI) 4k
(PI) 30k (SCL'XMM)
Epoch:
0.5k-1k
30k-30.5k
Last 0.5k
(PI) 7k
1.0
Relative L2 error
2
0.06
Density
0.04
0.4
0.02
0.2
0.0
0
0.00
0.0
0.5
1.0
10
30
20
25
30
B
B
(b)
(a)
Figure 1: Solving parametric convection BVPs: (a) relative L2 error vs. (legend reports number of differential operator evaluations per epoch), (b) Samples from vit at different training stages.
challenging scenarios (e.g., convection with / jointly adapts the weight of each loss and the points used to evaluate them during training (see Alg. [ Although R3 also targets points with high PDE residuals, does not approximate the worst-case loss needed to guarantee (weak) solution (Prop. B.T.. The improved performance of SCL comes at higher computational cost (= 5x), although this is largely offset when solving parametric problems.
non-linear eikonal). This is because
5.2 SOLVING PARAMETRIC FAMILIES OF BVPs
key advantages of based approaches is their ability to solve entire families of BVPs at once.
Consider fitting MLP (2r, 4, 9, to solve convection problems for 3 [1,30] using (SCL) (M)
and (PI) as in Cho et al. 2024).. We do not use their tailored "P"INN" architecture, although it would be compatible with SCT. The π 3, taken be 4, 7, and 30 equispaced values in [1330]. Fig.
shows that PP can only achieve the error of (SCL') hM) for the finest discretization, at which point evaluates the PDE loss times more per epoch. The same pattern holds for reaction-diffusion (4-7 times) and 2D Helmholtz (4-5 times) problems (App. F). This effectiveness is due to the worst-case loss of (SCL') jointly selecting (Ir, t. .3) target challenging coefficients as well as the domain regions responsible for that difficulty. In fact, inspecting viDD at different stages of training (Fig. (b) shows that SCL first fits the solution "causally." focusing on smaller values of While this has been proposed in (Krishnapriyan et al., [021] Penwarden al.) 2023 Wang et al.) 2024), it arises naturally by solving (SCL' (M). As training advances, however, Alg. shifts focus to fitting higher convection speeds B. This occurs without any prior knowledge of the problems manual tuning.
As we have argued, this is a use case for which SCL is particularly well-suited. Indeed, consider solving 2D Helmholtz equation for parameters (01, a2) [1, using FEM solver with mesh size chosen to obtain a similar accuracy as SCL. Across 100 experiments, the average relative L2 for the FEM solver was 0.036 for an average runtime of 4.3 minutes per solution (see App. FF) Hence, in the time it took to train the SCL model (31.1 hours), we could evaluate less than 440 parameters combinations using the FEM solver. This 20 times less than the 10* combinations used to evaluate the error of (SCL') (M) (average eиror 0.013 for a runtime of 4 minutes). It is worth noting that these numbers are for highly-optimized FEM implementation (Baratta al. (023).
5.3 LEVERAGING INVARIANCE WHEN SOLVING BVPs
Next, we showcase how SCL can be used to overcome computational limitations or searce mechanistic knowledge. Consider training an MLP to solve convection BVP with 30 and periodic initial condition h(0,, x) = ..). This problem is commonly used to showcase a "failure mode" (PI) since it yields degenerate solutions when using fixed collocation points (Fig. 2h) (Krishnapriyan et al. 2021). Note that the constrained (MC also fails when using fixed collocations points (Fig.2p), though its stochastic version finds accurate solutions (Table IT. One way overcome this limitation is by leveraging additional knowledge. this case, we know from the BVP structure that its solution must be periodic with period π 7/15. By incorporating this invariance using (SCL

---

Published as a conference paper at ICLR 2025
(PI)
(M)
(SCL)(M+S)
Analytical Solution
Magnitude
5.0
x 2.5
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
t
t
Figure 2: Using invariance in convection BVPs (BC and PDE losses
PI and (SCL) evaluated
using a fixed set of collocation points).
Table 2:
Relative Lg error test set (average across 10 seeds, see App. for standard deviation).
Dual Variable (A)
0.07
0.06
1
(PI))
SCC(()
Burgers'
10 10-3
0.0540 %
0.0444 %
10
4.29 %
3.31 %
0.04
Navier-Stokes
10 10-4
32.2 %
29.9 %
0.00
0.25
0.50
0.75
1.00
1.25 1.50
10 10-5
27.6 %
26.0% %
Magnitude of Initial Condition
Diffusion-Sorption
0.274 %
0.218 %
Figure 3: Final value of dual variables ())))) vs. magnitude of IC for Burgers' equation.
can obtain accurate solutions despite our use of fixed collocation points for (SCL (Fig. 2b). Note that the worst-case loss in (S) is fundamental to avoid degenerate solutions.
5.4 SUPERVISED SOLUTION OF BVPs
conclude by exploring supervised methods that rely only on pairs of initial conditions hy (ar, 0) and corresponding (weak) solutions u.. To showease the versatility of SCL, we train FNOs rather than MLPs, fixing m to a uniform distribution over regular grid accommodate their predictions. We also cast the SCL problem using only the observational constraints (S).)OO), i.e., without any objective. In contrast to (PII), which minimizes the average error, this formulation enforces maximum еиrог of across samples. Though apparently minor, this leads to better prediction quality, as shown Table 2 This occurs because the difficulty of fitting PDE solutions varies across ICs. While the average tends to emphasize the majority of "easy-to-fit samples," enforcing maximum error gives more weight to challenging data points. This becomes clear in Fig- which compares the magnitude of each IC in the training set with its final dual variables OB; for the Burgers' equation. Immediately, we notice a trend where ICs with either small or large magnitudes appear harder to fit. This information can be leveraged guide the collection of additional data points improve the NO architecture. Indeed, any model improvement can immediately take advantage of SCL since it is (virtually) independent of the choice of ug. The benefits of having O are clear when predicting which IC is hard to fit is intricate, as is the case for the Navier-Stokes equation (see App. F).
6 CONCLUSION
This paper developed SCL, technique for solving BVPs based on constrained learning. It demon- strated that finding (weak) solutions of PDEs is equivalent to solving constrained learning problems with worst-case losses, which also allows prior knowledge to be naturally incorporated to the solution of the BVP, e.g., structural constraints (e.g., invariances), real-world measurements, and previously known solutions. It developed practical algorithm tackle SCL problems and showcased its performance across a variety of PDEs and NN architectures. SCL not only yields accurate solutions, but also tackles many challenges faced by previous methods, such as extensive hyperparameter tuning, degenerate solutions, and fine discretizations domain and/or coefficients. Future work includes exploring worst-case losses for functional parameters such as and the of SCL for active learning.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
This work was partly funded by the Deutsche Forschungsgemeinnscaatt (DFG, German Research Foundation) under Germany's Excellence Strategy (EXC 2075-390740016), It was performed in part on the HoreKa supercomputer funded by the Ministry of Science, Research and the Arts Baden- Württemberg and by the Federal Ministry of Education and Research. Viggo Moro thanks the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for their support and G-Rescarch for supporting the travel costs.
REFERENCES
Tara Akhound-Sadegh, Laurence Perreault-Levasserr, Johannes Brandstetter, Max Welling, and Siamak Ravanbakhsh. Lie point symmetry and physics informed networks, 2023.
K. Arrow, L. Hurwicz, and H. Uzawa. Studies in linear and non-linear programming. Stanford University Press, 1958.
Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth N. Wells. DOLFINx: the next generation FEniCS problem solving environment, 2023.
Shamsulhaq Basir and Inane Senocak. Physics and equality constrained artificial neural networks: Ap- plication forward and inverse problems with multi-fidelity data fusion. Journal of Computational Physics, 463: 1 1301, 2022.
Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P Mailoa, Mordechai Korblluth, Nicola Molinari, Tess E Smidt, and Boris Kozinsky. E(3)-cquivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature communications, 13(1):2453, 2022.
D. Bertsekas. Convex optimization theory. Athena Scientific, 2009.
D. P. Bertsekas. Convex optimization algorithms. Athena Scientific, 2015.
Morteza Boroun, Zeinab Alizadeh, and Afrooz Jalilzadeh. Accelerated primal-dual scheme for class of stochastic nonconvex-concave saddle point problems. In American Control Conference, pp. 204-209, 2023.
S. Brenner and R. Scott. The Mathematical Theory of Finite Element Methods. Texts in Applied Mathematics. Springer New York, 2007.
Souvik Chakraborty. Transfer learning based multi-fidelity physics informed deep neural network. Journal of Computational Physics, 426: 109942, 2021.
Nithin Chalapathi, Yiheng Du, and Aditi S. Krishnapriyan. Scaling physics physics-informen informed hard constraints
with mixtuue-offexpprtts International Conference on Learning Representations, 2024.
L. F. O. Chamon and A. Ribeiro. Probably approximately correct constrained learning, In Advances in Neural Information Processing. 2020.
L. F. O. Chamon, S. Paternain, M. Calvo-Fullana, and A. Ribeiro. Constrained learning with non-convex losses. JEEE Trans. on Inf. Theory, 69[3):1739-116: 2023.
Yuyao Chen, Lu Lu, George Em Karniadakis, and Luca Dal Negro. Physics-informen neural networks for inverse problems in nano-optics and metamaterials. Opt. Express, 28(8): 1611-1166633 2020.
Minhao Cheng, Qi Lei, Pin- Yu Chen, Inderjit Dhillon, and Cho-Jui Hsieh. CAT: Customized adversarial training for improved robustness. In International Joint Conference on Artificial Intelligence, pp. 673-679, 2022.
Pao-Hsiung Chiu, Jian Cheng Wong, Chinchun Ooi, My Ha Dao, and Yew-Soon Ong. Can-pinn: A fast physics-iniee informed neural network based on coupled- ielhommtti-nuumetiii differentiation method. Computer Methody in Applied Mechanics and Engimeering, 395:1 14909, 2022.

---

Published as a conference paper at ICLR 2025
Woojin Cho, Minju Jo, Haksoo Lim, Kookjin Lee, Dongeun Lee, Sanghyun Hong, and Noseong Park. Parameterized physics- inffomed neural networks for parameterized PDEs. In International Conference Machine Learning, 2024.
Taco Cohen and Max Welling. Group equivariant convolutional networks. International Conference on Machine Learning, pp. 2990-2999, 2016.
Andrew Cotter, Heinrich Jiang, Maya Gupta, Serena Wang, Taman Narayan, Seungil You, and Karthik Sridharan. Optimization with non- con-differentiable constraints with applications to faimess, recall, chur,, and other goals. Journal of Machine Learning Research, 20(172):1-59) 2019.
Arka Daw, Jie Bu, Sifan Wang, Paris Perdikaris, and Anuj Karpatne. Mitigating propagation failures in physics-informed neural networks using retain-resampll--eleass (R3) sampling. In International Conference on Machine Learning, pp. 7264-7302, 2023.
Shaan Desai, Marios Mattheakis, Hayden Joy, Pavlos Protopapas, and Stephen J. Roberts. One-shot transfer learning of physics-infonnn informed neural networks. In for Science Workshop (ICML), 2022.
Guneet S. Dhillon, Kamyar Azizzadenesheli, Jeremy D. Bernstein, Jean Kossaifi, Aran Khanna, Zachary C. Lipton, and Animashree Anandkumar. Stochastic activation pruning for robust adver- sarial defense. In International Conference Learning Representations, 2018.
Mahesh Dissanayake and Nhan Phan-Thien. Neural-network-bassee approximations for solving partial differential equations. Communications in Numerical Methods in Engineering, 10:195-201, 1994.
J. Elenter, L. F. O. Chamon, and A. Ribeiro. Near-optimal solutions of constrained learning problems. In International Conference on Learning Representations, 2024.
L.C. Evans. Partial Differential Equations. Graduate studies mathematics. American Mathematical Society, 2010.
Rizal Fathony, Anit Kumar Sahu, Devin Willmott, and Zico Kolter. Multiplicative filter networks. In International Conference on Learning Representations, 2021.
Tanner Fiez, Lillian Ratliff, Eric Mazumdar, Evan Faulkner, and Adhyyan Narang. Global conver- gence to local minmax equilibrium in classes of nonconvex zero-sum games. In Advances in Neural Information Processing Systems, pp. 29049-29063 2021.
Alberto Fiorenza, Maria Rosaria Formica, Tomáš G. Roskovee, and Filip Soudsky. Detailed proof of classical Gagliardo-Nirenberg interpolation inequality with historical remarks. Zeitschrift für Analysis und ihre Anwendungen, 40(2):217-236, 2021.
Han Gao, Luning Sun, and Jian-Xun Wang. Phygeonet: Physics-informed geometry- convolutional neural networks for solving parameterized steady-stat pdes on irregular domain. Journal of Computational Physics, 428:1 10079, 2021.
Ian J. Goodfellow, Jonathon Shlens, and Christian Szegedy. Explaining and harnessing adversarial examples. In International Conference on Learning Representations, 2015.
Somdatta Goswami, Cosmin Anitescu, Souvik Chakraborty, and Timon Rabezuk. Transfer learning enhanced physics informed neural network for phase- modeling of fracture. Theoretical and Applied Fracture Mechanics, 106:10244 102447, 2020.
Tamara G Grossmann, Urszula Julia Komorowska, Jonas Latz, and Carola-Bibiane Schönlieb. Can physics-informed neural networks beat the finite element method? IMA Journal of Applied Mathematics, 89(1)1143-774 2024.
Gauray Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-bused operator learning for differential equations. Advances neural information processing systems, 34:24048-22 -24622, 2021.
Jayesh K Gupta and Johannes Brandstetter. Towards multi cullii-pattiotmmpral-sssse generalized PDE
modeling. Transactions Machine Learning Research, 2023.
12

---

Published as a conference paper at ICLR 2025
Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying. Yinpeng Dong, Songming Liu, Ze Cheng, Jian Song, and Jun Zhu. GNOT: A general neural operator transformer for operator learning. In International Conference Machine Learning, pp. 12556- 12569, 2023.
Jacob Helwig, Xuan Zhang, Cong Fu, Jerry Kurtin, Stephan Wojtowytsch, and Shuiwang Ji. Group equivariant Fourier neural operators for partial differential equations. arXiv preprint arXiv:2306.050m 05697, 2023.
Kurt Hornik. Approximation capabilities of multilayer feedforward networks. Neural networks, 4(2): 251-257, 1991
I. Hounie, L. O. Chamon, and A. Ribeiro. Automatic data augmentation via invariance-constrained learning. In International Conference Machine Learning, 2023a.
Hounie, A. Ribeiro, and L. F. O. Chamon. Resilient constrained learning. In Advances in Neural Information Processing, 2023b.
Ameya D. Jagtap, Kenji Kawaguchi, and George Em Karniadakis. Adaptive activation functions accelerate convergence in deep and physics-informed neural networks. Journal of Computational Physies, 404: 109136, 2020.
Spren Fiig Jarner and Ernst Hansen. Geometric ergodicity of metropolis algorithms. Stochastic Processes and their Applications, 85(2):341-361, 2000.
Namgyu Kang, Byeonghyeon Lee, Youngjoon Hong, Scok-Bae Yun, and Eunbyung Park. Pixel: Physics-informed cell representations for fast and accurate pde solvers. In AAAI Conference Artificial Intelligence, pp. 3186-8194, 2023.
J.T. Katsikadelis. The Boundary Element Method for Engineers and Scientists: Theory and Applica- tions. Elsevier Science, 2016
Michael Kearns, Seth Neel, Aaron Roth, and Zhiwei Steven Wu. Preventing faimess Gerrymandering: Auditing and learning for subgroup fairness. In International Conference Machine Learning. pp. 2564-2572, 2018.
Ehsan Kharazmi, Zhongqiang Zhang. and George E.M. Karniadakis. hp-vpinns: Variational physics- informed neural networks with domain decomposition. Computer Methods in Applied Mechanics and Engineering. 374:1 13547, 2021.
Diederik P. Kingma and Jimmy Ba. arXiv:1412.698009.
Adam: A method for stochastic optimization, 2017.
Knabner and L. Angerman. Numerical Methods for Elliptic and Parabolic Partial Differential Equations. Texts in Applied Mathematics. Springer New York, 2003.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research, 24(89):1-97, 2023.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael Mahoney. Char- acterizing possible failure modes in pnysics-informed neural networks. In Advances in Neural Information Processing, volume 34, pp. 26548-24222 265660, 2021.
I.E. Lagaris, A. Likas, and D.I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Trans. Neural Netw, 95):987711000 1998.
Randall J. Lequee Finite Difference Methods for Ordinary and Partial Differential Equations. Society for Industrial and Applied Mathematics, 2007.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. In Advances in Neural Information Processing Systems, pp. 6755-6766. 2020.
13

---

Published as a conference paper at ICLR 2025
Zongyi Li. Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference Learning Representations, 2021.
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. Journal of Machine Learning Research, 24(388):1-22 2023.
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. differential equations. ACM/ / IMS J. Data Sci., 1(3), 2024.
informed neural operator for learning partial
Tianyi Lin, Chi Jin, and Michael Jordan. Near-optimal algorithms for minimax optimization. In Conference Learning Theory, pp. 9738-2779, 2020.
Lu Lu, Pengzhan Jin, Guofei Pang. Zhonggiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218-229, 2021a.
Lu Lu, Raphaël Pestourie. Wenjie Yao, Zhicheng Wang, Francese Verdugo, and Steven G. Johnson. Physics-informee neural networks with hard constraints for inverse design. SIAM Journal on Scientific Computing. 43(6):B 105-1110--1111111 132,, 2021b.
Michael Lustig, David L Donoho, Juan M Santos, and John M Pauly, Compressed sensing MRI. IEEE signal processing magazine, 25(2):72-85 0277722-22222 2008.
Suryanarayana Maddu, Dominik Sturm, Christian L Müller, and Ivo F Sbalzarini. Inverse Dirichlet weighting enables reliable training of physics informed neural networks. Machine Learning: Science and Technology, 3(1):015026, 2022.
Stefano Markidis. The old and the new: Can physics-informed deep-learning replace traditional linear solvers? Frontiers in Big Data, 4, 2021.
Levi D. McClenny and Ulisses M. Braga-Neto. Self-adaptive physics-informed neural networks. Journal of Computational Physics, 474:111722 2023.
Nick McGreivy and Ammar Hakim. Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations. Nature Machine Intelligence, 6 110):1256-1269, 2024.
Aleksander Mądry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. Towards deep learning models resistant to adversarial attacks. International Conference on Learning Representations, 2018.
Sean Molesky, Zin Lin, Alexander Y. Piggott, Weiliang Jin, Jelena Vucković, and Alejandro W. Rodriguez. Inverse design nanophotonics. Nature Photonics, 12(11):659-670, 2018.
Ben Moseley, Andrew Markham, and Tarje Nissen-Meyer. Finite basis physics-ifformed neural networks (FBPINNs): A scalable domain decomposition approach for solving differential equations. Advances in Computational Mathematics, 49(4);62, 2023.
Daniel Musekamp, Marimuthu Kalimuthu, David Holzmüller, Makoto Takamoto, and Mathias Niepert. Active learning for neural PDE solvers. In Advances in Neural Information Processing. Workshop on Data-driven and Differentiable Simulations, Surrogates, and Solvers., 2024.
Mohammad Amin Nabian, Rini Jasmine Gladstone, and Hadi Meidani. Efficient training of physics- informed neural networks via importance sampling. Computer-Aided Civil and Infrastructure Engineering, 36(8))662299)) 2021.
Louis Nirenberg. On elliptic partial differential equations. Annali della Scuola Normale Superiore di Pisa-Scienze Fisiche Matematiche, 13(2): (1(2):115-40 162, 1959.
Akihiko Nishimura, David B Dunson, and Jianfeng Lu. Discontinuous Hamiltonian Monte Carlo for discrete parameters and discontinuous likelihoods. Biometrika, 107(2):365-380 2020.
14

---

Published as a conference paper at ICLR 2025
Peter J. Olver. Symmetry groups and group invariant solutions of partial differential equations. J. Differential Geometry, 14:497-542, 1979.
Ravi G. Patel, Indu Manickam, Nathaniel A. Trask, Mitchell A. Wood, Myoungkyu Lee, Ignacio Tomas, and Eric C. Cyr. Thermodynamically consistent physies-informed neural networks for hyperbolic systems. Journal of Computational Physics, 449:110755 2022.
Michael Penwarden, Ameya D. Jagtap, Shandian Zhe, George Em Karniadakis, and Robert M. Kirby. unified scalable framework for causal sweeping strategies for Physics-Informed Neural Networks (PINNs) and their temporal decompositions. Journal of Computational Physics, 493, 2023.
Raphaël Pestourie, Youssef Mrouch, Chris Rackauckas, Payel Das, and Steven G. Johnson. Physics- enhanced deep surrogates for partial differential equations. Nature Machine Intelligence. 5[12]: 2023.
Lee C.
Potter, Emre Ertin, Jason T. Parker, and Müjdat Cetin. Sparsity and compressed sensing in radar imaging. Proceedings of the IEEE, 98(6): 100--1020, 2010.
Dimitris C. Psichogios and Lyle H. Ungar. A hybrid neural network-firs principles approach to process modeling. Aiche Journal, 38: 199--111,,, 1992.
Md Ashiqur Rahman, Zachary E Ross, and Kamyar Avizzadenesheli. U-NO: U-shaped neural operators. Transactions on Machine Learning Research, 2023.
M. Raissi, Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
C.P. Robert and G. Casella. Monte Carlo statistical methods. Springer Verlag, 2004.
A. Robey*, L. F, O. Chamon*, G. J. Pappas, H. Hassani, and A. Ribeiro. Adversarial robustness with semi-infinite constrained learning. In Advances in Neural Information Processing, 2021.
R. T. Rockafellar and R. J-B Wets. Variational Analysis, volume 317. Springer Science & Business Media, 2004.
Hwijae Son, Jin Woo Jang, Woo Jin Han, and Hyung Ju Hwang. Sobolev training for physics informed neural networks. (00i::20100099::: 2021.
Olaf Steinbach and Marco Zank. Coercive space-time finite element methods for initial boundary value problems. Electronic Transactions on Numerical Analysis, 52:154-194 2020.
Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Dan MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. PDEBench: An extensive benchmark for scientific machine learning. In Advances in Neural Information Processing Track on Datasets and Benchmarks, 2022.
V. Thomee. Galerkin Finite Element Methods for Parabolic Problems. Springer Series in Computa- tional Mathematics. Springer Berlin Heidelberg, 2013.
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized Fourier neural operators. In International Conference on Learning Representations, 2023.
Chuwei Wang, Shanda Li, Di He, and Liwei Wang. Is L physics informed loss always suitable for training physics informed neural network? Advances in Neural Information Processing, 2022a.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks, SIAM Journal Scientific Computing, 4365))3305-A3381 2021a.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. On the cigenvector bias of Fourier feature networks:
From regression to solving multi-scale PDEs with physics-iffrmmed neural networks. Computer Methods in Applied Mechanics and Engineering, 384:1 13938, 2021b.
Sifan Wang. Xinling Yu, and Paris Perdikaris. When and why PINNs fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449:1 10768, 2022b.
15

---

Published as a conference paper at ICLR 2025
Sifan Wang, Shyam Sankaran, and Paris Perdikaris. Respecting causality for training physics- informed neural networks. Computer Methods in Applied Mechanics and Engineering. 421: 116813, 2024.
Min Wei and Xuesong Zhang. Super-resolution neural operator. In Proceedings the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 18247-18256, 2023.
Colby L. Wight and Jia Zhao. Solving Allen-Cahn and Cahn-Hilliard equations using the adaptive physies informed neural networks. Communications Computational Physics, 29(3):930 -54,, 2021.
Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non- adaptive and residual-based adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 403: 15671, 2023.
Dongxian Wu, Shu-Tao Xia, and Yisen Wang. Adversarial weight perturbation helps robust general- ization. In Advances in Neural Information Processing, pp. 2958-2969, 2020.
Chen Xu, Ba Trung Cao, Yong Yuan, and Günther Meschke. Transfer learning based physics- informed neural networks for solving inverse problems in engineering structures under different loading scenarios. Computer Methods in Applied Mechanics and Engineering, 405:115852, 2023.
Junchi Yang, Negar Kiyavash, and Niao He. Global convergence and variance reduction for class of nonconvex-nnonoocaae minimax problems. In Advances in Neural Information Processing Systems, pp- 1153-1165, 2020.
Junchi Yang. Xiang Li, and Niao He. Nest your adaptive algorithm for parameter-agnostee nonconvex minimax optimization. In Advances in Neural Information Processing Systems, 2022.
Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems. Computer Methods in Applied Mechanics and Engineering. 393, 2022.
16

---

Published as a conference paper at ICLR 2025
A APPLICATIONS OF (BP))
We showcase here variety of phenomena that can be described using (BVP). The PDEs detailed in this section are used throughout experiments to illustrate the performance of SCL.
A.1 CONVECTION EQUATION
The one-dimensional convection equation models the transport of a scalar quantity u(æ; 1), such as temperature or concentration, along the spatial dimension We consider the convection problem Ou(x, 1) Ōt + 3
Öu(r, 1)
= 0,
(sr, 1) E (0, (m, x (0, 1]
(6a)
dr
u(z, 0) - sin(x),
a E [0, 2T]
(6b)
u(0, u(2 t);
E (0,
(6c)
where the coefficient 3 denotes the convection speed. Despite its use of periodic BCs, namely, (6c).
this problem (and its variations) can be cast straightforward extension of (BVP).
A.2 REACTION-DIFFUSION EQUATION
The one me-dimensional reaction-diffie diffusion equation describes variety of phenomena, including chemi- cal reactions, population dynamics, and heat transfer, depending on the form of its reaction term. In this paper, we consider the reaction-diffusion problem with periodic BC and Gaussian IC. Explicitly, Qu(x,t) t)
Ou(x, t)
pu(x, t) u(xxtt)) f)),
(x, E (0, 2T) x (0,
(7a)
0t
Ox
1 1(---T)
u(I, 0) = exp 2 n7/
x E [0, 2xl
(7b)
(0,tt) t) = u(2TT, t),
E (0, 1]
(7c)
where v > 0 and p are the diffusion and reaction coefficients, respectively.
A.3 EIKONAL EQUATION
The eikonal equation is encountered in many applications involving wave propagation, e.g., elec- tromagnetism. It also describes the (signed) distance between any point x S2 and some fixed boundary äs, hence its usage in vision applications. In this case, we consider the BVP
Vu(x, v)ll = 1, (x, v) E
(8a)
u(z, y) 0, (x, as
(8b)
where Ω 1)2 and as is complex shape, in case, the gears figure from (Daw et al.
2023). In order to ensure that negative (positive) distances are assigned to the interior (exterior) of the shape, we add structural constraint to the solution which enforces that is non-negative on the boundary of Ω. Explicitly, we enforce that
u(z, v) 2 O, (z, v) 00.
This can be done by adding a structural constraint to (SCL). i.c., by replacing (SCL with a loss that induces 9 explicitly
E (( (EP
sup
[[ uo(x,))) S.
(10)
In our experiments, we find that this constraint is not particularly difficult to enforce and that we can obtain good results by replacing the worst-case vaa by fixed, equispaced points on 0Ω.
A.4 HELMHOLTZ EQUATION
The two-dimensional Helmholtz equation models wave propagation and vibration phenomena in various physical contexts. Here, consider the problem
Vuu(x, y) + u2u(x,, y) = Txx, y, T);
(x, y) E interion(S)
(1
u(x, v) sin(Taa1)) sin(maay))
(a, v) E as
(11b)
17

---

Published as a conference paper at ICLR 2025
where S 10,, ko 0 is the number; coefficients (a1. 02) represent the spa- tial frequencies in the x and y directions, respectively; and the forcing function T(Ø, y, T) (N22 n at x²)3) (sin(raja) sin(Ta2y))
10:
A.5 BURGERS' EQUATION
The one-dimensional Burgers' equation models the behavior of scalar field u(x,() () under the combined effects of nonlinear convection and diffusion. is commonly used in fluid dynamics and traffic flow to describe shock waves and turbulence. In particular, we consider the following BVP with periodic BCs
@ulx, 1) 1 (,,(1r 1)
1 ()
+
(x, t) E (0, x (0,
(12a)
δt
2
dx
Ox
u(2, 0) ho(z),
E [0, 1]
(12b)
(0,,t) t) u(1,t): t);
$ E (0,
(12c)
where is the viscosity coefficient, which governs the strength of the diffusion term. We consider ICs hg from the same distribution as et al. (002))
A.6 DIFFUSION-SORTIION EQUATION
The Inffusion-sorptiion equation models the transport of scalar field u(x, t) (e.g.. contaminant concentration) in porous medium. It is commonly used in environmental science and chemical engineering. terms of this PDE, we consider the following BVP
Ōu(x,t) t)
S Ou(x, t)
(2, (0, 1) (0,500)
(13a)
Ōt R(u(x,t)) ())
Ox
u(x,0) 0) = ho(r),
E [0, 1]
(13b)
0u(1,t)
u(0, t) = 1, u(1,t)
t E (0, 500]
(13c)
Ox
4 is the diffusion coefficient and R is the retardation factor defined where == x 1
R(u) = + ②) pak
where P3 = 2880 is the bulk density, 6 = 0.29 is the medium porosity, 3.5 10- Freundlich's sorption parameter, and n 0.874 is Freundlich's exponent. We consider the distribution of ICs ho from Takamoto et al. (2022).
A notable aspect of 13 is its nonlinearity due the dependence u(x, t) of the effective diffusion coefficient P R(a). What is more, it can become singular when O, making this equation particularly
challenging to solve.
A.7 SNVIER-STOEE EQUATION
The two-dimensional, incompressible Navier-Stokes equation in vorticity form removes the pressure term and focuses the dynamics of rotational flow. It is used to describe the local rotation of fluid, i.e., the vorticity w(r, t) = V x u(x: t), where is the two-dimensional velocity field and V f denotes the curl of f. Although vorticity is more challenging model than velocity, it offers deeper understanding of the flow dynamics. The velocity field can be recovered from the vorticity using Poisson's 's equation.
Explicitly, for x = [21, [2] we consider the BVP Ow(x, t) + u(I, 1)
Vw(a, t) 2 + TTx),
(x, E (0, x (0,7)
(14a)
01.
w(I, 0) Wi,
E [0, 1)2, (a,t) (0, 1)² x (0.
(14b)
u(x, t) O,
(14c)
18

---

Published as a conference paper at ICLR 2025
where v >> 0 is the viscosity coefficient and V f denotes the divergence of f. The forcing function taken to be
7(x) 0.1 sin (2 + z2)) + cos (2 + 22))
and the ICs wp are taken from the same distribution as in Li et al. 202 We consider three settings, namely, v 10-3 with 50, 10-4 with T 30, and
and 7 20.
19

---

Published as a conference paper at ICLR 2025
B WEAK SOLUTIONS AND ROBUST LEARNING
The space T of test functions plays a fundamental role when defining the weak formulations I Typically, it chosen to be a Sobolev space due to its natural compatibility with BVPs and the fact that it leads to less stringent differentiability requirements issues with the strong formulation (BVP). Sobolev space consists of functions in some Lebesgue space L whose weak derivafives also in L.. Recall that I,P is the space of p-integrable functions. To define a Sobolev space, we therefore need to start by defining a weak derivative.
It therefore overcomes the main
A locally integrable function f defined open set 2 C Rdd is weakly differentiable with respect to z; if there exists Df also locally integrable (i.e., f, Df E (2)) such that
J D((+)(()d== - (x)) dE(2) Oz:
dz, for all a (Z),
(15)
where C (2) is the space of infinitely differentiable, compactly supported functions. We say Df is the weak derivative of To generalize (15) to higher-order derivatives, consider the multi-
index = (ai,..., .e)) to be a d-tuple of non-negative integers and let |a| -Ei ag. We then define the a- weak derivatives of f, denoted Da f. as
[ D f(z)((2)dz = / (xz) Ozi ... Ozg dz, for all & (Z).
0 (()))
(16)
We can now define what we mean by Sobolev space Definition B.1 (Sobolev space). For an integer k 2 0 and 1 S S co, we define the Sobolev space (Z) = (S E (2( D f (Z) for all multi-indices a with |a| S k).
Evans 2010).
We write Wr. whenever the set 2 clear from the context. Note that since LP T for p 2 1, Sobolev spaces impose more restrictions than weak differentiability. Also, while W is in general Banach space, 5 is Hilbert space (Evans 2010).
Having set the groundwork, we can now proceed with the proof of Prop. 3.1
B. PROOF OF PROP. 3.1
The proof follows by constructing a measure of the deviation from the weak formulation | and showing that it dominated by the proposed worst-case statistical loss. This immediately implies Prop. 33.1. Explicitly, note that the weak formulation (1) can equivalently be expressed as
Up (Diulee, T(x t)) (0(x, t)dzdt = 0, for all yp E Wrk'2
2
(17)
where we omitted the dependence on the coefficients ı for conciseness. Using Jensen's inequality, we can upper bound (17) for any as in
(Dlulo,, - T(x, 1)) (x) so(ar, ()dodl
(cr, -(a,)))-((r)) S / (1x) (x,,,, t)dadl. (18)
Since to E Wk' c I??, the partition function
Z = Jp pxxx t) dxdt = 11 < oo is dell-de
We can therefore consider the normalized v = 222 /2p in (18) to get
-(x, U/, (D(,(,,))- p(z,, t)drdt t)ddd)
2
s Zo /
( "a t)dxdt.
(19)
Since va is a non-negative, normalized function, is the density of a probability measure, i.c., v E P. The following technical lemma shows that it is in fact in pa, ie., it is gquare-integrable probability density.
Lemma B.2. Let up E Wk. 2 and 1p 2² be probability distribution. Then, t p2
20

---

Published as a conference paper at ICLR 2025
Before proving Lemma
B.2
let us conclude the proof. The hypothesis u² implies that
J (Dlu)(z, r(x, 1x(x, t)drdt = 0 PEP2
sup
and since Zy is bounded, we obtain from
that
1² U/ (Diw))a,))) - r(z,t) pl t)drdt
S
Noticing that this holds for all s Wk' recover 17 which concludes the proof.
Proof of Lemma B.2 show P?, we must show that Φ E Lt. Indeed,
J (22,1))
1
||1|| | 1
72 plx,, 1) dadi
(20)
Since yp E W C L?, suffices it to show that the numerator finite. To do so, we can use the Gagliardo-Nirenberg interpolation inequality (Nirenberg 1959 Fiorenza al.) (021) write
| S "(( 1 -111l12)) ), for a [0, 1],
(21)
as long as (i) 2 m and (ii) 4am = d 1. additional condition must hold in particular cases:
(iii) if m d ++ 1 is a non-negative integer, then < 1.
Note that d 1 is the dimension of the space-time domain D. Since & W the right-hand of (21) is finite. We therefore only need to show that (i)-(iii) are satisfied under the hypothesis of the theorem. To do so, consider three cases:
(a) d = 0: in this case, d + is odd so that condition (iii) does not apply. Immediately, for and 0.25.
21 holds
(b) d = 1: since (ii requires
now have that d+1 is non-negative integer. Hence, condition (iii) applies. Nevertheless, we can still take m = and a 0.5 (II).
(c) i we can now consider other cases by taking
d and a
d +
m
4m
Indeed, (ii) holds by construction. Additionally, from the hypothesis of the theorem, have + 1 < 4k', which by the monotonicity of the ceiling operation implies that (i) holds. suffices to show that (iii) never applies. For 2 < d < 3, we have m and (d + 1)/2 > 1. For 3, we have
d ++1
d+1
d d+1
d+1
<
+1
<.
(22)
4
2
2
Thus, under the hypothesis of the theorem, implies that p2.
21 implies that ||s/|12. < co. From 20 this in turn
B.2 PROP.
For the sake of completeness, we provide a short discussion of the preliminary material needed to prove this proposition.
21

---

Published as a conference paper at ICLR 2025
B.2.1 PRELIMINARIES
We begin by showing that the supremum can be written as distributionally robust optimiza- tion problem. Indeed, consider the worst-case loss {(0) = (e.()))) [(ug(z, t)) and assume
that (z, t) E(us t)) is function in L This loss can be written in epigraph form
(())
inf t subject to [uu(2 t)) S t, for all (z, D.
(PV)
Writing PVV in Lagragian form (Bertsekas 2009 Ch. 4), we obtain
[()) - inf LE
sup EEL
IPV 0,4,))
(PVI)
where L3 denotes the subspace of almost everywhere non-negative functions of L3, Here, the Lagrangian Lpv(O, 0, L, 1) is defined as
LPV 1p (,) ++ / (1x (((uu() 8)) t dxdt
(23)
=t =t(1- - S (1x, t)dwdt + ()),oot)((o dx)dt,
Since 23) is linear function of t, ((0) the optimal value a linear program parametrized by Hence, strong duality holds (Bertsekas
2009 Ch. 4) and we obtain that
(0) sup WE
where
min tEl
Lpy 0,t, v).
(24)
Since is unconstrained and PV is linear in t, the dual function diverges -00 un- less fo (2x t)drdt = 1. From
23 and (4). we thus obtain that
{(0) sup / (ar ()E(ua(s, ))dddt +EPa
(25)
We also quickly review the necessary variational results for normal integrands. The majority of this exposition is adapted from Rockafellar & Wets (2004). Throughout, we let the tuple (T, denote measurable space, where T is nonempty set and A is a rlllebr of measurable sets belonging to T.
Definition B.3 (Carathéodory integrand). function T R R is called Carathéodory integrand if it is measurable in for each x and continuous in x for each t.
Definition B.4 (Decomposable space). A space FF of measurable functions g : T + R" is decomposable in association with measure μ on A if for every function go E F, for every set E A with (A)) < co, and for every bounded, measurable function 91 A + IR.,, the space F contains the function g R" defined by
(t(()) 99 (t)
for E T,A, for E A.
g(t)
(26)
Note that Lebesgue spaces L are decomposable for all E [1, oo] (see, e.g., (Rockafellar & Wets 2004 Ch. 14)). We can now state a crucial result concerning the interchangability maximization and integration.
Theorem B.5 (Thm. 14.60 in Rockafellar & Wets 2004)), Let F be decomposable space of measur- able functions and F T R be Carathéodory integrand. Then, as long as f(F, ((T) # for all é E F,
inf [ 6EF
f(T,
inf f(r,0)) dr. HERR
(27)
Moreover, as long as this common value is not oc, one has that
J f(r. ((T)) = Ō(T) E argmin f(T, for almost every e DEF
argmin
(28)
#ER"
22

---

Published as a conference paper at ICLR 2025
B.3 PROOF OF PROPOSITION 3.2
Proof. Consider a sequence ii E p converging to {(O) in L2, i.e., a solution of 25 In other words, for every 8 > O, there exists No such that
(0) / 2n (x, t)e(uo(( t))d 2 {()) - 6, for all n 2 Ng.
(29)
Consider now Gs minlssN | Since Vi E p C L²,
I 1 is finite for all n. And
since N₂ is finite, is cg. We can therefore rewrite
29 as
Jo b(z t)e(u))) t)) (0) = sup
vEta
(PVII)
S (b(x, ()dadt =1, 11vb1222 S C6.
subject to
where we rewrote 22 as {v L3 SD (((xt)dzdd 1}. Notice that (PVII is a convex quadratic
program in tb. Furthermore, note that a zero-mean normal distribution with variance /2/2 is strictly feasible for (PVII) (it belongs p and strictly satisfies the Lー-norm constraint). Hence, Slater's condition holds and we find that (PVII) is strongly dual (Bertsekas 2009 Ch. 4). We therefore conclude that for every cø > O, there exists μs IR and S 76 do such that
/ (1x ())(u ())drdt ou [/o (xx, t)dxdt 1 78 [1bbl cs)
(x( sup VEL2
/ (,(a,))((ua ()) rab(,, f)² + asb(x, f) dxdt 7808 05. wer2,
sup
To conclude, we use the fact that 2 is decomposable and since that the integrand is Carathéodory to exchange the supremum and the integral obtain that
Eo(
sup vil(uo(z,)) f)) + 7/8b/ asrb drdt 78C5 0x Sek
straightforward calculation the inner maximization problem shown above yields that the solution to (PVII) is given by
(z t) [e(uo (x, t))
-a]+
(30)
2y
where [2]+ max(0, 2) denotes the projection onto the non-negative orthant and a, ? are chosen so that
/ va (1r, t)drdt = and
|| I.= S C5-
Hence, for any 5 > 0 in 129 we can find such that 30 approximates f.
23

---

Published as a conference paper at ICLR 2025
C SAMPLING WITH THE METROPOLIS-HASTING ALGORITHM
Algorithm uses samples from one or more the in order to compute the losses in steps 4-7. This is, however, straightforward unless those distributions are fixed to, c.g., uniform (as we typically do for the BCs). That is because we only know vo up to normalization factor. We overcome this issue using MCMC techniques, more specifically, the -asttings algorithm (Robert & Casella 2004). In Alg. 2
we consider the general case of sampling from to x L, where É non-negative scalar-valued loss. We denote by En the desired samples and R. their support. In Alg. for instance, we would take in = (xn pile Eplee an
plee) R = D x 11, and
((zn) (Dnme [uer (TPple)
(xo :nhle) Txxfa epalle )²
(31)
in step and on = (2 (n) R D, and
((zn) = (noo (9) -
uf ((3)),)))
in step 7.
Typically, the covariance Σ is taken to be diagonal (independent proposals), e.g., 1²1. The choice of parameter o² of the proposal (step 3) affects the mixing rate, i.e. how fast the samples converge to the desired distribution. Smaller values of 0² will lead to slower mixing chains since the algorithm will not explore the space efficiently. On the other hand, large values will cause the acceptance probability in step 4 to be too small, so that the algorithm will remain stuck (step 5). Oftentimes, the parameter a² adapted during a burn-in phase hit specific acceptance rate, around 30% as rule-of-thumi (Robert & Casella, 2004). In experiments,
find reasonable value for o² and keep it fixed throughout training. We also consider a burn-in period by using only the last No samples generated by Alg. 2 For the single BVP experiments in Sec. [5.1 we use No 1000 and for the parameterized BVPs in Sec. 5.2 we use No = 2500.
Given that sample only from bounded domains (i.e., some subset of D II), the target dis- tribution th has finite tails for any a, satisfying sufficient conditions for uniform ergodicity [see, e.g., (Jarner & Hansen 20000] The law of the samples obtained from Alg. therefore converge (in Kullback-Leibler divergence) to via- Additionally, since we prove in Prop. 3.1 and vo (and more generally, v are square-integrable, alternative sampling technique with faster mixing rates can be used. That is the case, for instance, of Langevin Monte Carlo (LMC) (Robert & Casella 2004). Yet, the LMC algorithm uses first-order information of For the PDE loss in (31), this means higher-order space-time derivatives of ug and thus, additional backward passes. It is not clear that the benefits of faster mixing outweigh the increase in computational complexity, especially given that good results can be obtained using Alg. 2
Algorithm 2 Metropolis- Hastings algorithm with Gaussian proposal
1: ~ Uniform(R)
b Sample initial state
2: for O,..., N
3: 8 Gaussian (Zn, Σ) {(2)
b> Draw proposal
4: Pu = min 1, {(2)))
T(2 E R)
D Evaluate acceptance probability
5: 2n+1 2
with probability Pr.
Update state
7n+1 Zn, with probability Pn
6: end
7: return {21,..., NN}
24

---

Published as a conference paper at ICLR 2025
D GENERALIZATION RESULTS
Here, we formalize the generalization guarantces for when solutions of the empirical dual problem (probably approximately) near-optimal and near-feasible for the statistical primal problem. These were first detailed in [Chamon & Ribeiro (2020); [Chamon et al.] (2023). As done in Sec. 4 for simplicity and clarity, consider (M) Specifically, we are interested when solutions of DIV
are (probably approximately) near-optimal and near-feasible for (PIV). Generalization guarantees for the for the complete (SCL) as well as its parametric extension (SCL) follow in the same way.
We begin with the essential (non-convex) duality assumptions, In particular, we assume that the }} is sufficiently expressive (Assumption D.1) and that there exists a function ug E Fo that is strictly feasible (Assumption D.2). For NNs in particular, universal approximation theorems indicate that these assumptions are satisfied for large enough models (see, e.g.. (Homiik) 1991)) Finally, we impose a learning theoretic limit on the complexity of the hypothesis space Fe in order to ensure that our empirical approximations are well-posed (Assumption D.3). The main theorem our results are based namely (Chamon et al. 2023 Thm. 1), also require the losses to be convex, M -Lpischhitz continuous, and [0, B bounded. Since we only consider quadratic losses on the bounded domain D, these assumptions hold immediately.
hypothesis space Feg {ug
Assumption D.1. The parametrization us is rich enough that for each 02, 02 E θ and | E [0, there exists 0 @ such that (2)(), [8uo, (2r, () + (1 B)us, (x, 1) uo(a, 1)l v. Assumption D.2. There exist 6' such that ugr strictly feasible for PIV i.e., such that
(2.1)~ (Diuol(z, t) - r(x, c- Mv.
Assumption D.3. There exist C(N, 8) monotonically decreasing with N such that with probability 1 8 over samples (a 2 whe and (2ntle inde)) viplle holds for all E θ that
N
E E(a,t)-e (uo(x,))) - Mx,)))
(mo (h)) 1 )2w)))
S C(N,8)
N
E (0
(D)uo)(2, ((x,))))
Σ
(D(u))
(pde) T (a-fole ,
< <((N,
N
Under these assumptions, we can bound the empirical duality gap between (IV and = p - D*, where
DIV i.c..
p* minimize llee)
E(z,t)~aft (uo(x, t) h(z,t)))
subject to
E (2.))mott (D/uo)(z f) T(x, tx) sc
and
D* = max min L(0, A). X Ag0
Proposition D.4. Under Assumptions
D.I
D.3 it holds with probability (3m + 2) that
< 0 (Mu+
where X* is a solution of DIV.
Prop. D.4 is obtained directly from (Chamon et al. 2023) Thm. 1). This duality gap bound is enough to guarantee that the dual ascent algorithm in (5) provides near-optimal and near-feasible randomized solution of (PIV). Since all our losses are strongly convex (quadratic), we can further show that randomization is not necessary using the last iterate guarantees from (Elenter et al. 2024 Prop. 4... This is in spite of the fact that (PIV) a non-convex optimization problem.
On the other hand, the convergence of primal-dual methods such Alg. iin non-convex settings is the subject of active research, see, e.g., (Yang et al. 2020 Lin et al. 2020 Fiez et al., 2022]] Boroun et al. 2023). Transferring the guarantees from (5) to Alg. requires additional conditions, e.g., step size separation as in (Yang et al.) 2020)). Such convergence guarantees are, however, beyond the scope of this paper and left for future work.
25

---

Published as a conference paper at ICLR 2025
E
EXPERIMENTAL DETAILS
E.
HYPERPARAMETERS AND IMPLEMENTATION DETAILS
Throughout experiments, the relative L2 enror performance metric, which we define as
Σm us [T (Zn, (n) ul (T,, (1) bn)
Crel ( (T,
(32)
Σ
[ut ( (I,, kn)]]
where ut is the solution of VBP obtained either analytically by using classical numerical methods. For MLPs. the collocation points { } are taken from dense regular grid of points (see exact numbers below), and for FNOs, they are determined by the test sets from (Li al. 2021 Takamoto
et al. 2022). For parametrized problems, we report the average error
Σ (TI hz),
Erel J
evaluated either on a dense regular grid of points (for coefficients see exact numbers below) or based on the test sets from (Li et
2021 Takamoto et al. 2022).
To provide sensitivity measures, we run all experiments for 10 different seeds and report average and standard deviations of the results. We find that for certain difficult problem (e.g., diffusion with 3 = 50 or
with (v,, p) (3, 5)) the hyperparameters of SCL and R3 sometimes need to adjusted for certain seeds. This occurs rarely, but shows that there may not one-size- fits-all hyperparameter settings. For PINNs in (P]) we were unable to find any hyperparameters that solved those problems.
E.1.1 SOLVING SPECIFIC BVP (SEC. 5.1
In this section, we formulated SCL problems of the form (CL in order to use the same information that PINNs traditionally rely on. Recall that do use the worst-case distribution do (or even random points) for the BCs, but instead consider fixed, regularly distributed points. This reduces the overall computational complexity of the problem essentially performance cost. Explicitly, we consider the following problem
minimize GEB
uul (bc) 1i(2 be (12)
n=1
subject to
E (g,
(Diuel(z, - Txx, t)) Epde.
To compute the objective for the convection and reaction-diffusinn PDEs, we use 256 points (2 0), be [0, 2Tr). for the IC and 100 points equally spaced in E (0, 1] to evaluate the period BC. For the eikonal PDE, recall from (10) we use an additional structural constraint. In this case, we therefore formulate the SCL problem
M
minimize BEe
Σ
[Jo ymm)
M
subject to
E
(a [Dlue](s, ((x,)) (1))
Epde
Σ n=1
uo(ons 3n)] S.s
N
where fixed collocation points for the BCs and structural constraint, namely, M 2234 points on as (the gears figure from (Daw et al. 2023)) and 40 points on On. We use the exact same points for (PI)
26

---

Published as a conference paper at ICLR 2025
Problem hyperparameters. For SCL, the tolerance Epde selected by starting with small value (e.g., 10-1) and increasing it when the dual variables became too large during training to accommodate difficult problems. After a coarse hyperparameter search, we kept the weights j in (PI) used in (Daw et al. 2023). Note that eikonal PDE, since in this case the BCs play less critical role. All values are displayed Table B When solving the Eikonal equation with SCL, used 6, 10 the tolerance for the structural constraint.
used different weights for the BC and IC to solve the
We used MLPs with 4 hidden layers for ug each with 50 neurons for the convection and reaction-diffusion equations or 128 neurons for the eikonal equation and hyperbolic tangent activation function.
Model.
Training. To evaluate the PDE loss, all methods used 1000 collocation points sampled uniformly at random at the beginning of each epoch (PINN), obtained using the R3 from (Daw al., 2023) (R3), or using Alg. (SCL). For R3, we use the hyperparameters from (Daw et al. 2023). For Alg. 2 we use Σ diag 0.25, 0.01) for drawing proposals for x and respectively for both convection and reaction-diffusion. For the eikonal PDE, we use E = 0.04 T. In both cases, we run the algorithm for = 5000 and use only the last 1000 samples. All methods were trained using Adam with the default parameters from (Kingma & Ba 2017) and learning rates described in Table 4 Note that the baselines only use learning rate Tpp, since they do not use dual methods.
PDEs were tested on a dense regular grid of 256 x 100 points (x, t) E D against their analytical solutions. The solution of the eikonal PDE was tested on dense regular grid of 384 384 points (æx, y) E Ω against the ground truth predictions from (Daw et al. 2023).
Testing. The solution of the convection and
1.1.2 SOLVING PARAMETRIC FAMILIES OF BVPS (SEC. 5.2
The SCL problem we formulate here is similar the previous section, although we used the parameterized version (SL')(M). Once again, replace the (x,() marginals of the worst-case distribution v by fixed, uniform distribution. Note, however, that we keep the worst-case formulation for the coefficients Explicitly, we consider the SCL problem
No.
2
minimize 6E0
(o((T) (2)) h()x((()))))( (2))
Noc
subject to
(a.,,,))- FDE
E
[D=[uo(T)](:; (x) 1)) S Epde
Once again, we compute the objective for the convection and reaction-diffusion PDEs using 256 points (xr 0),
2rr) for the IC and 100 points equally spaced in t E (0, 1] evaluate the period BC. For the Helmholtz PDE, we use 4 x 256 points equally space around at to evaluate the BC. We the exact same points for (PI). Note that we include the coefficients in the forcing function to account for the Helmholtz BVP (see Sec.
Problem hyperparameters. Once again, the tolerance épde were selected by starting with a small value (e.g., 10 10-4) and increasing when the dual variables achieved too large value during training to accommodate difficult problems. The weights p in (PI) for the baselines were taken from (Daw Exact values are displayed in Table
al. 2023
Table 3: Problem hyperparameters for solving specific BVP
μD
μBBC
μАC
Epde
Convection: B =00
100
100
10-3
Convection: B 50
100
100
5 10-3
Reaction-diffusion: (v. p) (3,3)
100
100
10-2
Reaction- (v. = (3,5)
1
100
100
5 x 10-3
Eikonal
10
500
5 x 10-1
27

---

Published as a conference paper at ICLR 2025
Table 4: Training hyperparameters for solving specific BVP
Id
Learning decay
Iterations
(only SCL)
Convection: 30
10
10
0.9m every 5 000 iter.
175 000
Convection: / == 50
10
10
0.91 every 5 000 iter.
200 000
(v, p) = (3.3)
10 10-3
10
200.000
Reaction-diffusiin (v, p) (3.5)
10 10-3
10
200
Eikonal
10-3
10-4
0.9r every 5 000 iter.
60000
Model. We used MLPs with 4 hidden layers of 50 neurons.
we used 1000 collocation points per coefficient value, sampled uniformly at random at the beginning of each epoch. For (SCL. (M), we used Alg. [ For the PDE loss, used Σ = diag(0.25,, 0.01, o to sample from (Ir, L, π for both the convection (a,
Training. When training using PI..
for coefficient 3) and reaction-diffusion [0,,2 (1,1) for coefficients (V, equations. For the Helmholtz PDE, we used = 0.04 to sample from (x, t, T) for coefficients π = 02). The same variances of were used sample worst-case coefficients for the BC (recall that the distribution over collocation points is fixed). In all cases, SCL uses the last 2500 out of 5000 samples generated by the MH algorithm, except for the Helmholtz PDE with (01, a2) E [1, where we all 5, 000 samples to account for the additional difficulty of the problem.
All models were trained for 200, 000 using Adam with the default parameters from (Kingma Ba 2017) and learning rate of 10-3 for (PI). For SCL, the dual learning rate was T = 10 In all cases, we decayed the learning rates by factor of 0.9 every 5000 epochs, with the exception of the ceaction-diffusion PDE where we found it better to keep the learning rate constant.
Testing. The solution of the convection and reaction-diffusion PDEs were tested on a dense regular grid of 256 x 1000 points (z, t, m) E D II and 256 x00 x00x 100 points (æ, t, p) D II, respectively, against their analytical solutions. The solution the Helmholtz PDE was tested on dense regular grid of 256 256 x 100 x 100 points (r, y, d1, a2) E Ω II against its analytical solution.
E.1.3 LEVERAGING INVARIANCE WHEN SOLVING BVPs (SEC. 5.3
The SCL problems formulated in this section are of the form SCL (M+1). To showcase the advantages of integrating additional knowledge, such as the structure the BVP solution, we consider fixed collocation points for the constraints (M). This is in fact not uncommon for PINNs, see, e.g., (Raissi et al. 2019: Lu al. 2021b). These points are sampled uniformly at random once and then kept constant throughout training. Recall that for our convection BVP (Sec. A), the solution is
Table 5: Problem hyperparameters for solving parametric family of BVPs
PBC
Epdc
Convection
1
100
10 10-3
Reaction-diffusionn: (u,p) E 50,
1
100
5 10-3
Reaction-diifusion: (v.p) E [0, 10
100
10-2
Reaction-diffusion: (v,p) E [0, 20
100
10-1
Helmholtz: (a1, a2) E [1,2]
100
5 10
Helmholtz: (01, 02) [1,3]
100
5
28

---

Published as a conference paper at ICLR 2025
Table Problem hyperparameters for supervised solutions
# training samples
# validation samples
# test samples
FNO architecture
Burgers'
10-3
200
200
16 modes, layers
Diffusion sorption
10
1000
500
500
8 modes, 5 layers
Navier-Stokes: = 10 10-3
10-2
1000
500
500
8 modes, 8 layers
Navier-Stokess = 10
5 x 10 2
1000
500
500
8 modes, 8 layers
Navier-Stokes: 10-1
10-2
800
200
200
8 modes, 8 layers
periodic with period 2r/3. We therefore use the problem
minimize Ag
uo (.)) (()n)tt))
M
subject to
(iuol lmd) - T (Frre tmm)
I ppde
M
(mo(o,;)) t)
2T B
27
x +
S
For both (P and (SCL') (M) we use a total of N 456 collocation points, namely 256 points (z. 0),
[0, 2, for the IC and 100 points equally spaced in (0, 1] to evaluate the period BC. We use M = 100 collocation points sampled uniformly at random in the beginning of training and kept fixed throughout for the PDE loss.
Problem hyperparameters. For SCL, we take Cpde = 10 and = 10-3 For Pl). we use the weights y from (Daw et al. 2023), namely, μiD = 1, /JBC = 100, and uuc = 100.
Model. We used MIPs with 4 hidden layers of 50 neurons.
For the invariance loss, we used diag(0.5, 0.1) to sample from x..)). All models were trained for 200 000 epochs using Adam with the default parameters from Kingma & Ba 2017) and learning rate of 10-3 for (PIp. For SCL, the dual learning
Training. For SCL' (I), we used Alg. 2
rate 77 10 We decayed the learning rates by factor of 0.9 every 5000 epochs.
Testing. The solution was tested on dense regular grid of 256 x 100 points (z, t) E against its analytical solution.
E.1.4 SUPERVISED SOLUTION OF BVPs (SEC. 5.4
For supervised experiments, we formulate an SCL without objective using only data constraints (ob- servational knowledge). Since we FNOs, that can only make predictions uniform grids, replace Di in (SCL) with uniform distribution over fixed regular grid. The problem the FNOs tackle is that of predicting the solution ut of BVP given its IC f(z,, 0). Hence, the training data is composed of pairs (u,, h2) describing ICs and their corresponding solution. We therefore pose the SCL problem
minimize 9ee
N N
subject to
(En) (x) u (I bn)) eor j j= -.....,
For SCL, the tolerance was chosen as before, using a coarse hyperpa- rameter search. The final values reported in T
Problem hyperparameters.
Model. We used the ENO architecture from (Li et al. 2021 with 64 hidden channels, 128 projection channels, and no lifting channels. The number of modes and layers reported in Table
29

---

Published as a conference paper at ICLR 2025
Training and Testing. The datasets from (Li et al. equation, whereas the diffusion sorption dataset was taken from (Takamoto et al. 2022). All models
2021 were used for Burgers' and Navier-Stoken
were trained for 500 epochs using Adam with the default settings from (Kingma & Ba 2017) with
learning rate 10-3 and batch size of 20. For SCL., the dual learning rate was T) = 10-1 All learning
rates were decreased by factor of 0.5 every 100 epochs. All test errors reported for the model that achieved the lowest validation error during training. The sizes of the training, validation and test sets are reported in Table 6
30

---

Published as a conference paper at ICLR 2025
F ADDITIONAL EXPERIMENTS
F.1 SOLVING PARAMETRIC FAMILIES OF BVPS
We begin by presenting additional experiments focused on solving parametric families of BVPs (Sec. 5.2) and show how the samples from MH can be used to gain insights into the PDE and the training process.
In what follows, we report the "relative (computational) complexity" of (SCL in terms of differential operator evaluations per epoch. Explicitly,
# differential operator evaluations per epoch for (SCL)
Relative complexity
100%
# differential operator evaluations per epoch for (PI)
Recall that in order to evaluate the PDE loss, PI] uses 1000 collocation points per discretized
coefficient ny whereas (SCL.) takes 5000 steps of Alg. 2
Convection equation. Table considers simultaneously solving all BVPs corresponding to the convection equation with 3 I, 3.3 and compares (SCL') with (PI). We see that (SCL') outperforms or matches (P while being more efficient. In particular, (SL))(M significantly outperforms (Pī) in terms of relative L2 error for all but the finest discretization where they perform similarly. However, for that discretization, (SCL.) (M) is much more efficient that (PI). In that sense, it strikes better compromise between error and computational cost. This is even clearer from Fig. 4 particularly when we normalize the z-axis in terms of differential operator evaluations.
Table 7: Relative L2 error and computational efficiency for the parametric convection problem.
Average Relative Le Error
Relative complexity
Discretization for (PI
(Pi
(SCL')
SCL') ) + (PI)
{1.0, 10.0, 20.0, 30.0) (1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0)
0.365
125% 71 %
0.220
0.0110
2.0, 4.0, 5. 303
0.0476
16%
B = 10
B 20
B = 30
Relative La error
109
(PD)
(PI)
(SCL)(M)
(SCL)(M)
10
(PI)
(SCL)(M)
10
50k
100k
150k
200k
50k
100k
150k
200k
50k
100k
150k
200k
Epoch
Epoch
Epoch
β 10
β = 20
ß = 30
Relative L2 error
10a
(PD
(M)
10-1
(P1) (SCL)(M)
(PI) (SCL)(M)
10
0
1e9
1e9
1e9
Differential operator evaluations
Figure 4: Relative L2 error as function of training epoch and differential operator evaluations for the parametric convection problem.
31

---

Published as a conference paper at ICLR 2025
Reaction-diffusion equation. Additional results for the parametric reaction-diffusion BVP shown in 8 Same as for the convection PDE, (SCL") (M) is more efficient than (PT), striking a better compromise between computational complexity and performance. Indeed, in order to achieve the same error as ISCL"(M) (PI) requires between 5 and 7 times more evaluations of the PDE loss, i.e., of the differential operator D, per epoch. This is once again clear when looking at the evolution of the error during training (Fig. D). especially when the r-axis is displayed is terms of PDE evaluations. The distribution of errors across parameters is also more homogeneous for the SCL solution (Fig. 7
Finally, we can once again inspect the samples from vo throughout training to understand where the advantage of SCL comes from (Fig. o). First, we do not note any interesting behavior over the z-marginal (the samples are mostly uniform and the histogram is therefore omitted). Once again, we see that (M starts by focusing more on earlier times t, fitting the solution of the PDE "causally." Additionally, since the diffusion term tends to make the solution more homogeneous for larger times, it is clear that these are regions that are easier to fit and therefore require less attention. Once again, this behavior is not manually encouraged, but arises naturally from Alg. T As for the (v, p). we see that the distributions shift during training, indicating the change difficulty of PDE. In the end, the samples for p quite uniform,
fitting the solution of the while we notice that there remains a strong focus on smaller values of v. Note that these distributions reflect the error patterns of the final solution (Fig.
Table 8: Relative L2 error and computational efficiency for the parametric reaction-diffusion problem.
Coefficients range
Average relative L2 error
Relative complexity SSCL
Discretization for
(PI)
SCL (M)
5.0)
0.0793
55.6%
V/ E 0.
3.3,
0.0190
31.3%
0.0126
P E 0,5
{0.0, 1.25, 2.5, 2.5,, 5.03²
0.0119
20%
{0.0, 1.0, 2.0, 4.0.
0.0105
13.9%
{0.0, 5.0, 10.0)²
0.636
55.6%
E 00, 10 p E [0, 10
2.0, 7.5, 10.0)2
0.0228
0.0133
20%
{0.0, 2.0, 4.0, 6.0, 8.0, 10.03
0.0131
13.9%
E [1
10.0, 10.0, 20.0)²
0.841
55.6%
0.0204
E [1,20
(1.0, 5.0, 10.0, 15.0. 20.03
0.0128
20%
32

---

Published as a conference paper at ICLR 2025
(v, p) (1, 10)
(v, p) (5, 10)
(v, p) (10, 10)
Relative La error
10"
(PD)
(PI) (SCL))()
(PI)
(SCL)(M)
(SCL)()
10-1
11:8
50k
100k
150k
200k
50k
100k
150k
200k
50k
100k
150k
200k
Epoch
Epoch
Epoch
(v. p) (1, 10)
(v, p) = (5, 10)
(v, p) (10, 10)
Relative La error
109
1110
(PD
(PI)
(PI)
(MCL)(
(SCL)(()
(SCL)(M)
10
10-2
2
4
1e9
1e9
1e9
Differential operator evaluations
Figure 5: Relative L2 error function of training epoch and differential operator eval- uations for the parametric reaction-diffusion problem. (0.0, 6.0, 10.0)²
(PT) uses the discretization (u, E
Sampled y Epoch: 0.5k-1k
Sampled p
Sampled Times
Epoch: 0.5k- 1k
Epoch: 0.5k-1k 1k
0.100
0.075
0.050
0.025
0.000
Epoch: 30k-30.5k
Epoch: 30k-30.5k
Epoch: 30k-30.5k
Relative Frequency
0.10
0.05
0.00
Epoch: Last 500
Epoch: Last 500
Epoch: Last 500
0.05
0.00
0.0
0.2
0.4
0.8
1.0
Figure 6: Histogram of (marginalized)) MH samples of dh for the parametric reaction reaction-diffusion equation.
33

---

Published as a conference paper at ICLR 2025
(PI)
(SCL)(M)
4
10-1
Relative 4g error
3
10
V
(a)
(PI)
(SCL'((M)
10
Relative L2 error
10
4
2
(b)
Figure 7: Relative Lg error for reaction-diffusion solutions trained using SCL') and PI with
discretization (a) (v. p) E (0.0, 2.5, 5.0)2 and (b) (v. E (0.0, 2.0, 4.0, 6.0, 8.0, 10.0)
34

---

Published as a conference paper at ICLR 2025
Helmholtz equation. Once again, see from Table g that SCL' (M) makes efficient of computations than (P.. Indeed, in order to achieve the same error as (SCL') (M), (PI) requires between 3 and 4 times more evaluations of the PDE loss (i.e., of the differential operator D) per epoch. This is clear by looking at the evolution of the error during training after normalizing the r-axis in terms of computational complexity (Fig-5). Naturally, taking finer discretizations eventually leads to lower errors (Fig. D but the computational cost associated also rises considerably. On the other hand, keep the computational cost of (SCL') fixed throughout all experiments, showcasing its good performance across scenarios with little to no manipulation.
We can also inspect the samples from the throughout training to gain better understanding of the difficulties perceived by the MLP to fit solutions of this problem (Fig. TO). We display only the x and 01 marginals, seen as they display the same behaviors as y and 02 respectively due the symmetry of the Helmholtz equation. Here, we notice that the distribution of has an alternating pattern initially. This makes sense seen as the solution the Helmholtz equation is periodic. SCL clearly picks up this pattern, focusing the modes of the solution. As training continues, the sampling becomes more uniform, although with a focus on the boundaries of the domain where the notice
MLP clearly has difficulties fitting the solution. With respect to the problem coefficients, that th concentrates on larger values a1, especially in the beginning of training. These are indeed coefficients for which the solution of the problem is harder to fit (as evidenced by Fig. 1).
Table 9: Relative L2 error and computational efficiency for the parametric reaction-diffusion problem.
Coefficients range
Average relative L2 error
Relative complexity SCL') + (I])
Discretization for (PI)
(PI)
SCL (M)
1,5,,.5.33
0.0307
55.6%
d1 E [ E [
{1.0, 1.25, 1.5, 2.0
0.00593
0.0125
20%
1,, 1.6, 1.8,
0.00463
13.9%
(1.0,200, 3.0)²
1.34
55.6%
a1 E [1
3.0,
0.00943
0.0549
20%
02 E [1.3]
[1.0, 1.4, 1.8, 2.2, 2.2, 2.0,
0.00953
13.9%
(Ø1, a2) (1.0,
(α1, a2) = (1.5, 2.0)
(01, az) = (2.0, 2.0)
Relative L2 error
109
(PD)
(PI)
(PI)
(SCL)(M)
(SCL)(()
(SCL'XM)
10-1
10 2
50k
100k
150k
200k
50k
100k
150k
200k
50k
100k
150k
200k
Epoch
Epoch
Epoch
(a), 02) (1.0, 2.0)
(81, a2)
(1.5,
(81, a2) = (2.0,2.0)
Relative L2 error
10"
(PI)
(P
(M)
(SCL))(M)
(SCL)(M)
10-1
10
0
2
2
4
6
4
6
1e9
1e9
1e9
Differential operator evaluations
Figure 8: Relative L2 function of training epoch and differential operator evaluations for the
parametric Helmholtz problem. PI uses the discretization (2) E (1.0, 1.2, 1.4, 1.6, 1.8, 2.0
35

---

Published as a conference paper at ICLR 2025
(PI)
(SCL))(M)
2.0
0.05
0.04
Relative L2 error
0.03
1.4
0.02
1.2
0.01
0.00
1.0
1.2
1.4
1.6
1.8
2.0
1.0
1.2
1.4
1.6
1.8
2.0
a1
a1
(a)
(PI)
(SCL)(M))
Relative L2 error
1.6
10-2
82
1.4
1.2
1.0
1.2
1.4
1.6
1.8
2.0
1.0
1.2
1.4
1.6
1.8
2.0
A1
A1
(b)
Figure 9: Relative L2 error for Helmholtz solutions trained using (SCL') and PI with discretization
(a) (a1 a2) E {1.0, 1.5, 2.0}² and (b) (a1, E {1.0, 1.2, 1.4, 1.6, 1.8, 2.0
36

---

Published as a conference paper at ICLR 2025
Beginning (Epoch: 0.5k- 1k)
Beginning (Epoch: 0.5k- 1k)
0.10
Relative Frequency
Relative Frequency
0.08
0.08
0.06
0.06
0.04
0.04
0.02
0.02
0.00
0.00
Middle (Epoch: 30k-30.5k)
Middle (Epoch: 30k-30.5k)
Relative Frequency
0.08
Relative Frequency
0.10
0.06
0.08
0.06
0.04
0.04
0.02
0.02
0.00
0.00
End (Epoch: Last 0.5k)
End (Epoch: Last 0.5k)
0.125
Relative Frequency
0.10
Relative Frequency
0.100
0.08
0.075
0.06
0.050
0.04
0.02
0.025
0.00
0.000
-1.0
005
1.0
1.0
1.4
1.8
2.0
X
81
Figure 10: Histogram of (marginalized) MH samples of s/o for the parametric Helmholtz equation.
37

---

Published as a conference paper at ICLR 2025
F.2 SUPERVISED SOLUTION OF BVPs
Burgers' equation. Fig 11 shows the box plots for the relative L2 error across samples. They show that not only the average error distribution is shifted down. This is due to the variety of weights given to different samples, weights that in fact vary during the training process (Fig [2). The few large dual variables are related to ICs that are harder to fit and can provide important information for data collection or architecture improvements. We have already shown which ICs are harder for FNOs to fit Fig.
samples smaller when using (SCL) but in fact the whole
Box Plot for Train Errors
Box Plot for Test Errors
0.0010
Relative L2 error
0.0008
0.0006
0.0004
Eo
(PII)
(SCL)()
(PII)
(SCL)
Figure 11: Distribution of train and test errors (across data points) for the Burgers' equation (orange line indicates the median).
10a
Dual Variable (A)
0.06
Frequency
0.04
101
Maximum Minimum
0.02
Average
0.00
10a
100
200
300
400
500
0.07
Epoch
Dual Variable (A)
(a)
(b)
Figure 12: Dual variables obtained by Alg. 1 for the Burgers' equation: (a) evolution during training and (b) distribution the dual variables after training.
38

---

Published as a conference paper at ICLR 2025
Diffusion-sorption equation. We once again display the distribution of the relative L errors across training and test data points for models trained using (PII) and (S (Fig. [3p. Here, we clearly see that by bounding the maximum emror rather than minimizing its average leads to more homogeneous fit across samples. This once again due to the different weight assigned to each data
sample, weights that also evolve throughout training (Fig 14b). By inspecting the ICs with large and small values of X, we notice the pattern showcased in Fig 14p, where IC with either large or small magnitude are more challenging to fit than those with moderate ones.
Box Plot for Train Errors
Box Plot for Test Errors
0.012
0.001
Relative La error
0.010
0.008
0.006
0.004
(PII)
(SCL)()
(PII)
(SCLXO)
Figure 13: Distribution of train and test (across data points) for the diffusion-sorption equa-
tion (orange line indicates the median).
Dual Variable (A)
Dual Dual Variable (A)
0.07
0.04
0.06
Maximum
0.02
Minimum
0.05
Average
0.00
100
200
300
400
500
0.00
0.05
0.10
0.15
0.20
Epoch
Magnitude Initial Condition
(a)
(b)
Figure 14: Dual variables obtained by Alg. 1 for the diffusion diffusion-sorptionn equation: (a) evolution during training and (b) value function of the IC magnitude.
39

---

Published as a conference paper at ICLR 2025
Table 10: Relative L2 error on test set (mean standard deviation).
V
(PI
(S)
Burgers'
10
0.0540 + 0.0027 %
0.0444 + 0.0020 %
10-3
4.29 = 0.40 %
3.31 = 0.16 %
Navier-Stokes
10-4
32.2 + 0.87
29.9 = 0.54
10
27.6 + 0.63 %
26.0 = 0.33 %
Diftusion-Sorption
0.274 + 0.049 %
0.218 + 0.036 %
Navier-Stokes equation. We start by displaying an extended version of Table 2 including standard deviations of results over runs to show the consistency of our results across random seeds (Table(TO). We then turn to Fig 15 which shows that (SCL)(O) not only improves the average relative L2 error, but its entire distribution across train and test data points. For the Navier-Stoke equations, however, is harder find relation between IC properties and the difficulty fitting the solution. That because, we show in Fig [6] the dual variables do not have such extremely different values. This is certainly due to the fact ICs similarly difficult to fit. Still, some ICs have outlier values of À (Fig [6p), which does point to the fact that the FNO does struggle more to certain conditions. That being said, not easy to identify what in those conditions make them hard (Fig 17 Nevertheless, this is not issue
the tolerance ca set very loose (0.01) and that these situations,
need not know beforehand which ICs are challenging: suffices to run Alg. T to solve (O).
Box Plot for Train Errors
Box Plot for Test Errors
0.0060
Eo 0.01
0.10
Relative La error
0.0055
0.08
0.0050
0.06
0.0045
0.0040
0.02
0.0035
(PII)
(SCL)
(PII)
(SCL)
Figure 15: Distribution of train and test errors (across data points) for the Auvier-Stocke equation
with 10 (orange line indicates the median).
Dual Variable (
10a
0.06
Frequency
Maximum Minimum Average A
0.04
10
0.02
10a
100
200
300
400
500
0.00
0.01
0.02
0.03
0.04
0.05
0.06
Epoch
Dual Variable
(a)
(b)
Figure Dual variables obtained by Alg. for the Navier-Stokes Stokes equation with 10-3, (a) evolution during training and (b) distribution of the dual variables after training.
40

---

Published as a conference paper at ICLR 2025
IC for Smallest Dual Variables (1)
IC for Largest Dual Variables (1)
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
0.0
0.6
0.8
0.0
0.2
0.6
0.8
X
X
0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
Magnitude
Figure 17:
Initial conditions corresponding the smallest and largest final dual variables for Navier- Stokes equation with v = 10
41

---

Published as a conference paper at ICLR 2025
G ADDITIONAL RELATED WORK
G.1 PIYSICS-INFORMED NEURAL NETWORKS
Fitting MLPs to the solution of BVPs goes back to (Psichogios Ungar 1992, Dissanayake & Phan-Thien 1994. Lagaris al. 1998). The advent of differentiable programming and automatic differentiation, however, lead to an increased interest in this approach, which has since been used
to tackle both forward and inverse problems involving a variety of PDEs (see, e.g., (Raissi et al.
2019:
Wight & Zhao 2021 Chen et al., 2020; Lu et al., 2021b Basir & Senocak, 2022; Yu et al.
2022
Xu et al. (023)) ). This led to new architectures tailored for PDEs (Raissi
2019 Fathony
et al.
2021 Gao et al. 2021 Wang et al. 2021a Kang et al. [203] Moseley et al. 2023 Cho et al.
2024
Chalapathi et al. 2024) leveraging positional embedding (Wang et al. 202 lb) and adaptive
activation functions (Jagtap et al. 2020).
Training PINNs. PINNS tend to be very sensitive to training hyperparameters, particularly the choice of collocation points and loss weights. Many works have theoretically and empirically
investigated the origins of these issues (Krishnapriyan et al. 2021 Markidis [021] Wiigtt & Zhao
2021 Wang et al. 2021a 2022b a Grossmann al. (024). Based these observations, adaptive heuristics have been proposed to select the collocation points based on importance sampling (Nabian
al. 2021 Wu al. 2023), adversarial training Wang al. (2022a), rejection sampling (Daw et al.
2023). and iassaeety-inppired rules (Penwarden et al.) 2023 Wang et al. 2024). Similarly, empirical rules for determining the loss weights [/ in (Ph)) have been developed using the magnitude of the gradients (Wang et al. (021a), cigenvalues of the neural tangent kernel (Wang et al. [2022b), inverse- Dirichlet weighting (Maddu et al., 2022), soft attention mechanisms (McClenny & Braga-Neto 2023), and (augmented) lagrangian formulations (L.u et al. 2b22bb] [Basir & (annocak]2022) Other
works have addressed these challenges by changing the problem formulation inspired by traditional
numerical methods (Kharazmi al. 2021 Chiu et al. 2022, Patel et al. 2022), proposing different
objective functions (Yu et al. 2022; Son et al. 2022], and using sequential (Wight & Zhao 2021
Krishnapriyan et al. (021) and transfer learning (Goswami et al. 2020 Chakraborty 2021
Desai
ct al. 2022) techniques.
In contrast these approaches, we address these issues jointly by using worst-case losses and constrained learning to obviate these hyperparameters. In fact, we prove that the constrained learning problems we pose yield (weak) solutions of BVPs. Hence, it is not enough to either adversarial training to estimate the worst-case loss as in (Wang et al. (022a) or constrained learning to manipulate the loss weights as in (Lu al. 2021b Basir& & Senocak 2022). Both are required simultaneously. Leveraging findings from adversarially robust learning. we also replace the gradient methods used in (Wang et al. (0022)) i.e., the technique from (Mądry et al. 2018), by the sampling- -baee approach
in (Robey et al. 2021).
G.2 NEURAL OPERATORS
In contrast to the MLPs and convolutional NNs (CNNs) typically used in PINNs, NOs are NNs capable of handling BVP solutions for different IC or forcing functions. Many different architectures have been proposed. such as DeepONets (Lu et al., (021a). FNOs (Li et al., 202T), and NO based on U-Nets (Gupta & Brandstetter 2023). FNOs particular have become quite popular and garnered many efforts towards addressing their limitations, such as improving memory efficiency (Rahman et al. 2023). designing equivariant FNOs (Helwig et al. 2023), extending FNOs to general geometries (Li ct al.
inputs and outputs. They can therefore be trained to find
2023). ) factorizing the Fourier transform (Tran et al. 2023). and leveraging multiwavelets (Gupta
al.
2021)
Training NOs. Regardless of these improvements, the vast majority of NOs are trained in supervised manner by minimizing their average error across samples as in (PII (see, e.g., (Lu et al., 202 Li et al. 2020 Kovachki et al. 2023)). Unless substantial domain knowledge has been used during data collection, challenging cases may be underrepresented in the dataset, which could hinder the accuracy of the NO. Although semi-super techniques involving PDE losses have also been used (Li et al. 2024), computing the space-time derivatives needed to evaluate D, [u] is challenging for NOs.
42

---

Published as a conference paper at ICLR 2025
In this paper, do not develop new NO architectures, but focus the problem of training them. Explicitly, rather than targeting the average eтrоr, we target the maximum error across samples. This much better suited to handle the heterogeneous difficulty in fitting the data. We also incorporate structure in the solution during training, without the need to design new architectures.
G.3 CONSTRAINED AND ADVERSARIALLY ROBUST LEARNING
The main tool used in the development of SCL is constrained learning, or more specifically, robustness- constrained learning. Constrained learning is technique to train MI systems under requirements,
such as fairness (Kearns et al. 2018 Cotter al. 2019) Chamon Ribeiro
2020 Chamon et
2023) and robustness (Chamon & Ribeiro 2020 Robey* et al. [021] Hounie al. 2023a Chamon ct al. 2023), or to handle applications in which we want to attain good performance with respect to than metric. As in unconstrained learning, is formulated statistical risk minimization problem, albeit with constraints. Despite its non-convexity in virtually every modern ML. task, certain duality properties hold when using sufficiently expressive parametrizations, leading practical learning rule with generalization guarantees (Chamon & Ribeiro 2020) Chamon et al. (202)). These duality results have also been exploited automatically adapt each constraint specification their underlying difficulty, striking better compromises between objective and requirements Hounie et al. (2025b).
Aside from dealing with the nominal accuracy vs. robustness trade-off typical in ML systems, constrained learning has itself been used optimize robust losses. Typically, this is done using some combination of gradient ascent and random initialization, restart, and pruning heuristics (Goodfellow (et al. 2015 Mądry et al.) 2018 Dhillon al., 2018: Wu al.) 20200 Cheng et al. 2022). Using semi-infinite optimization techniques, however, these deterministic methods can be replaced by sampling approach that has been successful variety of domains (Robey* et al. 2021). Different
MCMC methods (Robert & Casella 2004) have been used in this context, including Langevin Monte
Carlo (LMC) (Robey et al. 2021) and Hettroooli--Hastinge (MH) (Hounie et al. 2023a).
43


Reviews:
[{"review_id": "RIpTc7j2cK", "reviewer": "Reviewer_JwZx", "summary": "This paper aims to adapt ideas from adversarial robustness and use them to find neural network solutions to partial differential equations with improved guarantees over existing neural solutions to PDEs.  The scheme, as I understand it, is to observe that the weak form of the solution to a PDE can be achieved by finding the minimum of a kind of \"worst-case PINN loss\" in which the distribution over collocation points is chosen to maximize the expected PINN loss.  This framing allows a particular result from adversarial robustness to apply: the worst-case (wrt distributions on inputs) expected loss on a region can be approximated by an expectation under a distribution proportional to the shifted-and-truncated loss on that region.  The amount of shifting governs the accuracy of the approximation.  I understand this to be a kind of continuity argument: if the worst-case distribution is concentrated on the points with the highest loss, then smoothing those out a bit doesn't change the value all that much.\n\nThe result is that you can get a reframing of the weak problem to look sort of like a PINN with a particular distribution on the collocation points.  The paper discusses how various kinds of constraints, e.g., symmetries, can then be framed in this way.  The actual sampling itself of the \\psi_\\alpha distribution is done with Metropolis-Hastings, or it is discarded in favor of a uniform distribution.  The approach is empirically evaluated on a couple of small problems and compared to some neural network baselines. No comparisons to conventional PDE solvers are performed.", "questions": "What did you use to compute the ground truth solutions?  How long did it take to compute those ground truth solutions relative to the methods compared in the paper?  What I'd like to see here is a thoughtful comparison with conventional methods and points along the Pareto curve trading off compute time against accuracy (here, that means the fineness of discretization for, e.g., FEM).", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 4, "contribution": 2, "strengths": "I think the main strength of this paper is that it allows one to reason more formally about the PINN-type setup and connect it directly to the weak form solution.  It's a creative idea to adapt the result from adversarial robustness. The paper is mostly well written.", "weaknesses": "The main weakness of the paper is that we don't get any sense of how well this actually works for solving PDEs.  Like many \"deep learning for PDEs\" papers, it does not compare to conventional methods for solving PDEs, but only other neural networks.  It is currently unclear whether PINNs and related ideas are actually ever a good idea.  Very few papers examine the actual Pareto curve of computational cost versus accuracy with respect to, e.g., FEM.  A recent paper reviewing the area has shed light on just how bad this situation is from a scientific perspective:\n\nMcGreivy, Nick, and Ammar Hakim. \"Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations.\" Nature Machine Intelligence (2024): 1-14.\n\nThe MH24 paper criticizes weak standard numerical baselines, but the present paper does not compare to any conventional numerical methods at all.  Do we conclude from this paper that I should stop using my PDE solver and start using this method?  The abstract claims \"accurate solutions\" but then only demonstrates that relative to deep learning baselines on tiny problems.\n\nAnd to be clear, this situation is not really the authors' fault.  They are just following the trend and applying the standard that the neural PDE solver community has created.  Nevertheless, the for the larger scientific community to find this work valuable, we must apply a standard that makes sense for people who want to actually solve PDEs.\n\nSo we might instead take the contribution of the paper to be the guarantee, except:\n\n1) As stated in the paper, existing methods already have guarantees.\n\n2) The guarantee itself centers on sampling from \\psi_\\alpha which, as mentioned in the paper, becomes difficult exactly when the error bound becomes tight.  Using Metropolis--Hastings essentially throws that guarantee away unless you can prove, e.g., uniform ergodicity of the resulting Markov chain. So in the end we are faced with what looks like just another difficult integral, just with a different form and a lack of a guarantee.  Note that one could raise the same criticism of the original adversarial robustness result: sampling from little regions around the worst case inputs is just as difficult as finding them in the first place.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "rQrImM1wYn", "reviewer": "Reviewer_x438", "summary": "The paper develops a technique called ``Science constrained learning'' in order to address some of the issues with Physics Informed Neural Networks (PINNs) by allowing structural constraints and known solutions to be taken into account in addition to the usual constraints considered by PINNs.", "questions": "I would be very pleased to hear the authors' response to my comments above, especially with respect to the proof of the reduction in conditioning number in Appendix B.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Overall the paper is well written and formatted and the development of Science Constrained Learning (SCL) is well described. There are several different background techniques described in the paper, in addition to the new technique. I think this work is worthy of inclusion in the conference and I am leaning towards an acceptance.", "weaknesses": "There are several different background techniques described in the paper, in addition to the new technique. Therefore, I think it would be beneficial to the readability of the paper to include a clear list of the scientific contributions of SCL to help the reader. I note that the paper is a full 10 pages long, however I think that some of the background could be shortened somewhat (especially in the introductory sections) in order to give the description of the new technique a bit more room to breathe. \n\nOne of the main claims for the new method is that it is not sensitive to the balancing of weights for the different components of the loss function. This is well justified by the use of the dual formulation. However this is not mentioned until the end of the 4th page. Due to the importance of this aspect of the solution, I would recommend mentioning it earlier in the paper. Otherwise the reader is left wondering what the difference is to e.g.\\ the work in [1], which uses derivative constraints and empirical observations. \n\nThe use of the Lagrangian does indeed support the idea that the method addresses the issue of balancing the weights in the loss function which can affect PINNs. However the paper also makes the statement that this new methodology addresses sensitivity to the number of training points. I cannot really see a clear explanation in the paper for why this might be the case. Either this explanation should be brought out more and clarified in the paper or, alternatively, it should be explained that this is an empirical observation.\n\n[1] Kentaro Hoshisashi, Carolyn E Phelan, and Paolo Barucca. No-arbitrage deep calibration for volatility smile and skewness. arXiv preprint arXiv:2310.16703, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "rYAJimiVM6", "reviewer": "Reviewer_4CoR", "summary": "Neural network-based approaches, such as physics-informed neural networks and neural operators, offer a mesh-free alternative by directly fitting those models to the PDE solution. They are known to be highly sensitive to hyperparameters such as collocation points and the weights associated with each loss. This paper addresses these challenges by developing a science-constrained learning (SCL) framework. It demonstrates that finding a (weak) solution of a PDE is equivalent to solving a constrained learning problem\nwith worst-case losses.", "questions": "The specific questions (associated with the items above) are:\n\n+ Can the authors clarify how their approach (for instance Algorithm 1 updating schedule) compares to the update schedule in Paris' paper:   https://epubs.siam.org/doi/pdf/10.1137/20M1318043    Section 2.4 (not the specific updates, but the general updates as given by the spirit of the logic presented)?  In others of Paris' papers, he points to :   https://arxiv.org/pdf/2007.04542.  How does the new method compare to the \"adding weights to the loss function\" section 2.2.1 (Which is analogous to Paris' work also)?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The strength of the paper is their proposed solution to the hyperparameter issue:  They develop a science-constrained learning (SCL) framework. It demonstrates that finding a (weak) solution of a PDE is equivalent to solving a constrained learning problem\nwith worst-case losses.\n\nThe paper provides ways of incorporating structural knowledge, observational knowledge, and science-contrained learning.  The algorithm is easy to follow, and the mathematical proofs are helpful in terms of understanding the benefits of the method.", "weaknesses": "The reviewer is having a difficult time distinguishing the contribution from the work of Paris Perdikaris and his group.  Examples are:\n\nA. Daw, J. Bu, S. Wang, P. Perdikaris, A. Karpatne, Rethinking the Importance of Sampling in Physics-informed Neural\nNetworks, arXiv preprint arXiv:2207.02338\n\n S. Wang, S. Sankaran, P. Perdikaris, Respecting causality is all you need for training physics-informed neural networks,\narXiv preprint arXiv:2203.07404\n\nParis references in one of his talks a paper by this former advisor and colleagues \"A unified scalable framework for causal sweeping strategies for ...\" that has a categorization (and references).   \n\nThe important point is that the algorithm in the current paper (Algorithm 1) has many of the same characteristics you find in Paris' work which you incorporate learning the hyperparameters into the learning.\n\nAt some point, the reviewer is wondering how common this is now considered when application papers now reference the constrained paper (NIPS) as say that they use PINNs + Constrained learning:\n\nhttps://www.sciencedirect.com/science/article/pii/S1385894724012919#b78\n\nIs the current paper really 'novel' (i.e. a contribution) if others now consider things commonplace.  The reviewer is open to being convinced of the contributions if the literature search were broader and the comparisons where against the various adaptive sampling methods mentioned by Paris and others.\n\nIn terms of concrete actions, can the authors:\n+  Provide a more comprehensive literature review that clearly positions their work relative to existing constrained learning approaches for PINNs?\n\n+  Explicitly discuss how their method differs from or improves upon other adaptive sampling techniques, particularly those by Perdikaris and others (as highlighted by the question below)?\n\n+ Highlight any novel theoretical results or empirical improvements over state-of-the-art methods", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IU0KNu5fFe", "reviewer": "Reviewer_RnPx", "summary": "This paper developed SCL, a technique for solving BVPs based on constrained learning. \nIt then developed a practical algorithm to tackle these problems and showcased its performance across a variety of PDEs.\nSCL not only yields accurate BVP solutions, but tackles many of the challenges faced by previous methods, such as extensive hyperparameter tuning and computational costs.", "questions": "- The proposed formulation seems to be a rephrasing of PDE boundary value problems with invariance, observation data into a single constraint problem. The resulting Lagrangian problem is then solved straightforwardly in Algorithm 1, reminiscent of standard multi-task learning algorithms. Is there more to it than this? \n\n- (Lines 214-) \"Compared to penalty methods, (Dˆ -CSL) does not require extensive tuning of coefficients [such as μ in (PI)] and guarantees generalization individually for each constraint (near-optimality and near-feasibility) rather than for their aggregated value.\": Could you elaborate on this? It seems to be a key statement regarding the paper's contribution. Could you provide some references to support this?\n\n- (Lines 58-) \"Instead, we use semi-infinite constrained learning techniques to develop a hybrid sampling-optimization algorithm that tackles both problems jointly without extensive hyperparameter tuning and training heuristics (Sec. 4).\":  I am not fully convinced by this claim. Algorithm 1 still includes hyperparameters, such as $\\eta_p$ and $\\eta_d$, and training heuristics are still necessary.\n\n- (Lines 472-) \"it (causality) arises naturally by solving (SCL′)(M). As training advances, however, Alg. 1 shifts focus to fitting higher convection speeds β. Note that this occurs without any prior knowledge of the problems or manual tuning.\": Could you clarify why causality naturally arises in this context?\n\n- Under what conditions, does Algorithm 1 converge?  (Chamon et al., 2023; Elenter et al., 2024) could be valuable references, to my understanding.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The authors tackle important problems in physics-informed learning: reducing computational costs and avoiding hyperparameter tuning.\n\n- The experimental results show that the proposed formulation outperforms baselines in some cases.", "weaknesses": "- The key contribution, Algorithm 1, can be seen as a variant of DNN training with learnable loss weights, which can often be found in multi-task learning, lacking novelty.\n\n- The dependence on random seeds is unknown, potentially causing reproducibility issues. Error bars should be added.\n\n- The paper is difficult to follow. The paper benefits from paragraph writing, and please present information in relevant paragraphs and sections. It was challenging to discern whether certain sentences conveyed key ideas or were merely supplementary notes.\n\n- I have several critical questions. Please see Questions below. I am open to discussion.\n\n- (Minor) (Lines 1730, 1733, etc.) Fig -> Fig.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)