Solving Differential Equations with Constrained Learning

## Abstract

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.

## 1 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 phenomena (Lustig et al., ; Potter et al., ; Molesky et al., ; Evans, ). As such, a variety of numerical methods have been used to approximate their solutions, such as the celebrated finite element method (FEM). Despite their approximation guarantees, the accuracy of these methods is often in proportion to the computational complexity associated with finer discretizations. Additionally, they do not naturally incorporate prior knowledge, such as real-world measurements or previous solutions to similar equations. Lastly, any change to the PDE problem, such as initial condition or mesh size, requires the solution to be recomputed (Brenner & Scott, ; LeVeque, ; Katsikadelis, ).

Methods based on neural network (NN) architectures,
such as physics-informed NNs (PINNs) (Lagaris et al., ; Raissi et al., ; Lu et al., )
and neural operators (NOs) (Li et al., ; Lu et al., ; Rahman et al., ),
have been proposed in an attempt to address these challenges. Rather than discretizing the PDE, these mesh-free approaches directly fit NNs to its solution. What is more, they can be trained to simultaneously solve families of PDEs and leverage data from measurements or known solutions by incorporating additional losses in their training objectives (Li et al., ; Cho et al., ; Li et al., ). Yet, the effectiveness of these methods is highly sensitive to hyperparameters such as the number of points used to evaluate the PDE and the weights used to combine the different training losses. These issues hinder their accuracy, especially for PDEs with non-smooth solutions (Krishnapriyan et al., ; Wight & Zhao, ; Wang et al., ; ). Though heuristics have been proposed to mitigate these issues, e.g., Nabian et al. (); Maddu et al. (); Wu et al. (); Daw et al. (), they are often *ad hoc*, derived for specific PDEs, and may strike poor compromises between performance and computational complexity.

In this paper, we address these drawbacks using *constrained learning*. In particular, we prove that obtaining a (weak) solution of a PDE is equivalent to solving a constrained learning problem involving worst-case losses (Sec. 3).
Hence, it is not enough to use adversarial training to estimate the worst-case losses (Wang et al., ) *or* constrained formulations to manipulate the loss weights (Lu et al., ; Basir & Senocak, ).
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).
This approach also naturally incorporates prior scientific knowledge on the structure (e.g., invariance) and observations (e.g., measurements) of the solution, which is why we dub it *science-constrained learning* (SCL).
We illustrate the effectiveness of SCL for a diverse set of PDEs, NN architectures (e.g., multilayer perceptrons, MLPs, and NOs), and prior knowledge levels. Our experiments show that it consistently yields accurate solutions for entire families of PDEs, sometimes at a lower computational cost (Sec. 5).

## 2 Problem Formulation

### 2.1 Boundary Value Problems

Consider a (bounded, connected, open) region $\Omega\subset\mathbb{R}^{d}$ with (smooth) boundary $\partial\Omega$ Consider a (partial) differential operator $D_{\pi}$ with coefficients $\pi\in\Pi\subset\mathbb{R}^{q}$ defined on the domain $\mathcal{D}=\Omega\times(0,T]$. Here, $\pi$ captures a (finite) set of parameters for the phenomenon, such as the diffusion rate or viscosity. Given a space $\mathcal{F}$ of functions mapping $\mathcal{D}$ to $\mathbb{R}$, we define a boundary value problem (BVP) as

$$ $\displaystyle\operatorname{find}$ $\displaystyle u\in\mathcal{F}$ (BVP) $\displaystyle\operatorname{such\ that}$ $\displaystyle D_{\pi}[u](x,t)=\tau(x,t),$ $\displaystyle(x,t)\in\Omega\times(0,T]$ (PDE) $\displaystyle u(x,0)=h(x,0),$ $\displaystyle x\in\bar{\Omega}$ (IC) $\displaystyle u(x,t)=h(x,t),$ $\displaystyle(x,t)\in\partial\Omega\times(0,T]$ (BC) $$

where $\tau:\mathcal{D}\rightarrow\mathbb{R}$ is the *forcing function* and $h:\mathcal{B}\rightarrow\mathbb{R}$ describes the boundary (BC) and initial (IC) conditions over $\mathcal{B}=\bar{\Omega}\times\{0\}\cup\partial\Omega\times(0,T]$ for $\bar{\Omega}=\Omega\cup\partial\Omega$. In what follows, we always consider the initial condition as part of the boundary conditions and refer to them jointly as BC. Note that the results 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., (Evans, ), for a detailed treatment of BVPs.

In general, a (strong) solution of (BVP) need not exist in any function space $\mathcal{F}$. This motivates the rise of relaxations such as the *weak formulation* that replaces the pointwise equation (PDE) in (BVP) by the integral equation(^1^11Typically, the weak formulation only considers the spatial domain $x$, handling $t$ separately using time-stepping methods (Thomee, ). Still, (1) is not uncommon and can be used to obtain weak formulations for parabolic PDEs (see, e.g., (Knabner & Angerman, ; Evans, ; Steinbach & Zank, )).)

$$ $\int_{\mathcal{D}}D_{\pi}[u](x,t)\varphi(x,t)\,dxdt=\int_{\mathcal{D}}\tau(x,t )\varphi(x,t)\,dxdt,\quad\text{for all }\varphi\in\mathcal{T}\text{.}$ (1) $$

The $\varphi$ are known as *test functions* and $\mathcal{T}$ is typically taken to be a *Sobolev space* due to its natural compatibility with this setting (see App. B for further details). The name *weak formulation* comes from the fact that a solution of (BVP) (when it exists) is also a solution of (1), although the converse is not necessarily true. Indeed, (1) allows for a wider range of solutions with less stringent regularity requirements, particularly with respect to continuity and differentiability (Evans, ). The BCs of (BVP) can often be homogeneized (i.e., $h\equiv 0$) and imposed implicitly through $\mathcal{T}$, thus fully describing its weak solution by (1) (Brenner & Scott, ; LeVeque, ; Katsikadelis, ).

### 2.2 Solving boundary value problems

There exists a wide range of numerical methods for solving BVPs, most of which rely on discretizing BVP (e.g., finite difference method, FDM) or its weak formulation (e.g., finite element method, FEM). Some approaches, such as the boundary element method (BEM), also take advantage of the PDE structure (Green’s function) to reduce the computational complexity associated with handling fine discretizations  (Brenner & Scott, ; LeVeque, ; Katsikadelis, ). Though ubiquitous in scientific and engineering applications, these classical methods only tackle one BVP at a time and do not naturally incorporate prior knowledge, such as measurements or known (partial) solutions.

To address these challenges, methods based on parameterizing the solution $u$ of BVPs have been proposed, most notably using NNs. While they may not achieve the precision of classical methods (Krishnapriyan et al., ; Wight & Zhao, ; Grossmann et al., ), they are able to provide solutions for whole families of BVPs, extrapolating new solutions from existing ones, and leverage extrinsic information, such as real-world measurements. Generally speaking, they can be separated into *unsupervised methods*, that seek to solve (BVP) directly, and *supervised methods*, that leverage previously computed or measured solutions.

##### Unsupervised methods.

This approach directly fits a model $u_{\theta}:\mathcal{D}\times\Pi\to\mathbb{R}$ with parameters $\theta\in\Theta\subset\mathbb{R}^{p}$ to (BVP) (e.g., a multilayer perceptron, MLP). An eminent example of such methods are physics-informed neural networks (PINNs), which train $u_{\theta}$ by solving

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}\ \mu_{\text{D}}\ell_{ \text{D}}(\theta)+\mu_{\text{BC}}\ell_{\text{BC}}(\theta)$ (PI) $$

for fixed weights $\mu_{\text{D}},\mu_{\text{BC}}\geq 0$, where

$$ $\displaystyle\ell_{\text{D}}(\theta)$ $\displaystyle\triangleq\sum_{i=1}^{I}\Bigg{[}\frac{1}{M}\sum_{m=1}^{M}\Big{(}D _{\pi_{i}}[u_{\theta}(\pi_{i})](x_{m},t_{m})-\tau(x_{m},t_{m})\Big{)}^{2}\Bigg {]}\text{,}\quad$ $\displaystyle(x_{m},t_{m})\in\mathcal{D}\text{,}\ \ \pi_{i}\in\Pi\text{,}$ $\displaystyle\ell_{\text{BC}}(\theta)$ $\displaystyle\triangleq\frac{1}{N}\sum_{n=1}^{N}\Big{(}f_{\theta}(x_{n},t_{n}) -h(x_{n},t_{n})\Big{)}^{2}\text{,}$ $\displaystyle(x_{n},t_{n})\in\mathcal{B}\text{.}$ $$

Note that we write $u_{\theta}(\pi)(x,t)$ to emphasize that we evaluate $u_{\theta}(\pi)$, that approximates the solution of (BVP) with coefficients $\pi$, at $(x,t)\in\mathcal{D}$. In practice, however, the input of the model $u_{\theta}$ is actually $(x,t,\pi)\in\mathcal{D}\times\Pi$.
The losses $\ell_{\text{D}}$ and $\ell_{\text{BC}}$ promote the requirements in (PDE) and (IC)–(BC) from (BVP), respectively. However, in contrast to the numerical techniques used by classical methods (e.g., FEM), the former is computed via automatic differentiation. Some formulations also include an additional data term used to incorporate known solutions (Raissi et al., ; Lu et al., ). Though the vast majority of PINNs target a single BVP, i.e., $I=1$ in (PI), their extension to parameterized families of BVPs has also been explored (Cho et al., ).

Though effective in many applications, (PI) is very sensitive to its hyperparameters. Indeed, the number of collocation points $(x,t)$ (typically sampled uniformly at random), the choice of PDE coefficients $\pi_{i}$ (typically, a uniform grid), and the weights $\mu_{\text{D}},\mu_{\text{BC}}$ (often hand-tuned) affect both the quality and computational complexity of (PI). These issues are exacerbated if the BVP has non-smooth solutions even for a few values of $\pi$ (Krishnapriyan et al., ; Wight & Zhao, ; Wang et al., ). Adaptive heuristics for collocation points (Nabian et al., ; Daw et al., ; Wu et al., ) and weighting rules (Wang et al., ; Maddu et al., ; McClenny & Braga-Neto, ) have been put forward to mitigate these challenges. Yet, they generally lack theoretical grounding, focusing on specific “failure modes,” and seldom address the non-trivial interactions of these hyperparameters. Hence, they may strike poor compromises between performance and computational complexity and their effectiveness often depends on the BVP.

##### Supervised methods.

Rather than directly solving the BVP, supervised approaches fit the model $u_{\theta}$ to a set of (partial) solutions $u^{\dagger}_{n}$ of (BVP). In this setting, however, $u_{\theta}$ is typically a neural operators (NOs) capable of handling infinite-dimensional (functional) inputs and outputs. As such, they can model the solutions of a BVP for, e.g., different forcing functions $\tau$ or IC $h(x,0)$. Given, e.g., forcing-solution pairs $(\tau_{n},u^{\dagger}_{n})$, these NOs are trained by solving

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}\ \frac{1}{N}\sum_{n=1} ^{N}\big{\|}u_{\theta}(\tau_{n})-u^{\dagger}_{n}\big{\|}_{L_{2}(\mathcal{D})}^ {2}\text{.}$ (PII) $$

In practice, the functions $\tau_{n}$ and $u^{\dagger}_{n}$ are discretized (in the time or spectral domain) to enable computations (Li et al., ; Lu et al., ; Hao et al., ; Wei & Zhang, ). Note that MLPs can also be trained using supervised methods (e.g., (Raissi et al., )), but they cannot directly handle infinite-dimensional inputs and are therefore restricted to BVPs with finite-dimensional coefficients $\pi$ (as in, e.g., (Cho et al., )). Similarly, training NOs by directly using (BVP) as in (PI) is challenging as it is not straightforward to compute $D[u_{\theta}]$ (Li et al., ).

While supervised methods can operate without explicit knowledge of the BVP, acquiring the (partial) solutions they are trained on using classical methods (e.g., FEM) or real-world measurements can be expensive. What is more, obtaining good performance from small datasets can be challenging since some BVPs are highly sensitive to their coefficients and IC (see, e.g., Navier-Stokes in App. A). This issue is aggravated by the heterogeneous difficulty of fitting different solutions, which suggests the uniform weighting in (PII) is suboptimal. As with (PI), the performance of (PII) is also affected by the discretization used to evaluate its objective, though empirical evidence suggests this is less of an issue than with unsupervised methods (Li et al., ; Hao et al., ; Wei & Zhang, ).

## 3 Science-Constrained Learning

In this section, we address the drawbacks of the previously introduced NN-based BVP solvers using *constrained learning* Chamon & Ribeiro (); Chamon et al. (). Seen as (BVP) is itself a constrained problem, this approach yields a more natural formulation. Although constrained formulations of (PI) have been proposed (Lu et al., ; Basir & Senocak, ), they fall short of guaranteeing a (weak) solution. Indeed, we prove that the weak form (1) leads to a *different* constrained problem. We then demonstrate how other scientific knowledge can be integrated into this mechanistic description, such as the structure of the solution (e.g., invariances) or (partial) observations obtained from measurements or classical numerical methods. In Sec. 4 and 5, we show that this approach leads to a straightforward algorithm that performs well without training heuristics.

We next introduce the framework used to cast these SCL problems (Sec. 3.1) and then show how mechanistic, structural, and observational knowledge can be formulated as learning objectives, i.e., statistical losses (Sec. 3.2). In the remainder of this paper, we use $u_{\theta}$ to refer to any parameterized model, regardless of its architecture (e.g., MLP, NO, etc.).

### 3.1 Constrained and robust learning

Consider distributions $\mathfrak{D}_{i}$ over data pairs $(z,y)\in\mathcal{Z}\times\mathcal{Y}$ together with losses $\ell_{i}:\mathbb{R}^{k}\times\mathcal{Y}\to\mathbb{R}_{+}$, $i=0,\ldots,m$. Constrained learning seeks a model $\phi_{\theta}:\mathcal{Z}\to\mathbb{R}^{k}$ parameterized by $\theta\in\Theta$ that solves

$$ $\displaystyle P^{\star}=\min_{\theta\in\Theta}$ $\displaystyle\operatorname{\mathbb{E}}_{(z,y)\sim\mathfrak{D}_{0}}\!\Big{[} \ell_{0}\big{(}\phi_{\theta}(z),y\big{)}\Big{]}$ (P-CSL) $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(z,y)\sim\mathfrak{D}_{i}}\!\Big{[} \ell_{i}\big{(}\phi_{\theta}(z),y\big{)}\Big{]}\leq\epsilon_{i},\quad i=1, \ldots,I.$ $$

Typically, we think of the objective $\ell_{0}$ as a top-line metric (e.g., accuracy), whereas the other $\ell_{i}$ together with $\epsilon_{i}$ define requirements of the solution (e.g, robustness or invariance) Cotter et al. (); Chamon & Ribeiro (); Chamon et al. (). A common alternative to (P-CSL) is to combine the $\{\ell_{i}\}$ into a single training objective [as in, e.g., (PI)], leading to the *penalty methods* ubiquitous in ML. Though the penalty and constrained approaches are sometimes related (e.g., they are in fact equivalent in convex optimization (Bertsekas, )), this is not the case in the non-convex settings of ML. Still, under typical conditions from (unconstrained) learning theory and for rich enough parametrizations, the solution of (P-CSL) can be approximated by solving the *empirical dual problem* (Chamon & Ribeiro, ; Chamon et al., ; Elenter et al., )

$$ $\displaystyle\hat{D}^{\star}=\max_{\lambda_{i}\geq 0}\ \min_{\theta\in\Theta} \ \hat{L}(\theta,\lambda_{i})\text{,}$ (D̂-CSL) $$

where $\hat{L}$ is the *empirical Lagrangian* of (P-CSL) based on samples $(z_{n_{i}},y_{n_{i}})\sim\mathfrak{D}_{i}$, namely,

$$ $\hat{L}(\theta,\lambda_{i})\triangleq\dfrac{1}{N_{0}}\sum_{n_{0}=1}^{N_{0}} \ell_{0}\big{(}\phi_{\theta}(z_{n_{0}}),y_{n_{0}}\big{)}+\sum_{i=1}^{I}\lambda _{i}\Bigg{[}\dfrac{1}{N_{i}}\sum_{n_{i}=1}^{N_{i}}\ell_{i}\big{(}\phi_{\theta} (z_{n_{i}}),y_{n_{i}}\big{)}-\epsilon_{i}\Bigg{]}\text{.}$ (2) $$

Compared to penalty methods, (D̂-CSL) does not require extensive tuning of coefficients [such as $\mu$ in (PI)] and guarantees generalization individually for each constraint (near-optimality and near-feasibility) rather than for their aggregated value. What is more, the duality between (P-CSL) and (D̂-CSL) endows the solution $\lambda^{\star}_{i}$ of (D̂-CSL) with a sensitivity interpretation, despite their non-convexity (Chamon & Ribeiro, ; Chamon et al., ; Hounie et al., ). Explicitly, $\lambda^{\star}_{i}$ describes how much the objective $P^{\star}$ would decrease for small relaxations of $\epsilon_{i}$. Finally, (D̂-CSL) is amenable to practical algorithms such as dual ascent, which updates $\lambda_{i}^{(0)}=0$ as

$$ $\lambda_{i}^{k+1}=\lambda_{i}^{k}+\eta\Bigg{[}\dfrac{1}{N_{i}}\sum_{n_{i}=1}^{ N_{i}}\ell_{i}\big{(}\phi_{\theta^{\dagger}_{k}}(z_{n_{i}}),y_{n_{i}}\big{)}- \epsilon_{i}\Bigg{]}_{+}\text{,}\quad\text{for}\ \ \theta^{\dagger}_{k}\in\min _{\theta\in\Theta}\ \hat{L}(\theta,\lambda_{i}^{k})\text{.}$ (3) $$

Even if the empirical Lagrangian minimizer $\theta^{\dagger}_{k}$ is only computed approximately, (3) can be shown to converge to a neighborhood of a solution of (D̂-CSL) (Chamon et al., ; Elenter et al., ).

The “in expectation” constraints from (P-CSL) may be too lenient for certain applications. This is the case, for instance, for (adversarially) robust learning problems, where requirements are cast in terms of the worst-case loss

$$ $\bar{\ell}(z)=\max_{\tilde{z}\in\mathcal{R}(z)}\ \ell(\phi_{\theta}(\tilde{z}) )\text{,}$ (4) $$

where $\mathcal{R}(z)$ is a compact *robustness set*. For simplicity, we omit the label $y$ in (4) as it can be incorporated in the loss $\ell$. Although many algorithms have been proposed to tackle learning problems involving $\bar{\ell}$ (e.g., (Goodfellow et al., ; Mądry et al., )), it turns out to be closely related to the expected constraints in (P-CSL). This fact is captured by the next proposition adapted from (Robey* et al., , Prop. 3.2), where $\mathcal{P}^{2}$ denotes the space of square-integrable probability distributions.

###### Proposition 3.1 .

Suppose that $\tilde{z}\mapsto\ell(\phi_{\theta}(\tilde{z}))\in L^{2}$ for $\tilde{z}\in\mathcal{R}(z)$. Then, for all $\delta>0$ there exists $\alpha<\bar{\ell}(z)$ such that
$\bar{\ell}(z)=\sup_{\psi\in\mathcal{P}^{2}}\ \operatorname{\mathbb{E}}_{\tilde
{z}\sim\psi}\!\big{[}\ell(\phi_{\theta}(\tilde{z}))\big{]}\leq\operatorname{
\mathbb{E}}_{\tilde{z}\sim\psi_{\alpha}}\!\big{[}\ell(\phi_{\theta}(\tilde{z})
)\big{]}+\delta$,
where $\psi_{\alpha}\in\mathcal{P}^{2}$ is supported on $\mathcal{R}(z)$ and $\psi_{\alpha}(\tilde{z})\propto\big{[}\ell(\phi_{\theta}(\tilde{z}))-\alpha
\big{]}_{+}$ for $[a]_{+}=\max(0,a)$.

A proof is provided in (Robey* et al., , App. C). Prop. (3.1) shows that the worst-case loss (4) can be arbitrarily well approximated by an expectation with respect to $\psi_{\alpha}$, a distribution proportional to a truncation of the underlying loss $\ell$. Hence, by sampling from $\psi_{\alpha}$, we can use (D̂-CSL) to also tackle robust learning problem. In Section 4, we investigate when and how this can be done efficiently. But first, we show how to cast (BVP) as a constrained learning problem that can also incorporate structural and observational prior knowledge.

###### Proposition 3.1 .

### 3.2 The science-constrained learning problem

In this section, we explore how mechanistic, structural, and observational knowledge can be combined to find solutions of BVPs. For clarity, we first focus on a single BVP instance, omitting the dependence on $\pi$ and/or $\tau$. We consider these extensions at the end of the section.

##### Mechanistic 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) to a robust learning constraint. This constraint is then written in terms of an expected loss using Prop. 3.1. Finally, the BCs are reintroduced using (P-CSL). We start with the following proposition, where $W^{k,p}$ refers to the $(k,p)$-th order Sobolev space (see App. B).

###### Proposition 3.2 .

Let $u^{\dagger}\in W^{k^{\prime},2}(\mathcal{D})$, where $k^{\prime}\geq 1$ is the degree of the differential operator $D$, be such that
$\sup_{\psi\in\mathcal{P}^{2}}\operatorname{\mathbb{E}}_{(x,t)\sim\psi}\!\Big{[
}\big{(}D[u^{\dagger}](x,t)-\tau(x,t)\big{)}^{2}\Big{]}=0$.
If the dimension $d$ of $\Omega$ satisfies $d\leq 4k^{\prime}-1$ , then $u^{\dagger}$ satisfies (1) with $\mathcal{T}=W^{k^{\prime},2}(\mathcal{D})$.

A proof can be found in Appendix B. Prop. 3.2 elucidates the challenges arising from the choice of collocation points in unsupervised approaches such as (PI), most notably PINNs. Indeed, it is not enough to use a fixed distribution (e.g., uniform): obtaining a weak solution requires training against an infinite number of them. Using Prop. 3.1, however, the infinite problem of satisfying the weak form (1) can be reduced to a single statistical constraint, namely,

$$ $\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}}\!\Big{[}\big{(}D[u_{\theta} ](x,t)-\tau(x,t)\big{)}^{2}\Big{]}\leq\epsilon\text{,}$ (5) $$

for $\psi_{\alpha}$ as in Prop. 3.1 and $\epsilon>0$. Note that Prop. 3.2 describes a sufficient condition for $u^{\dagger}$ to satisfy (1). A necessary condition can also be obtained by restricting the set $\mathcal{P}^{2}$ to sufficiently smooth distributions, namely, those belonging to a Sobolev space.

Yet, it is not enough to satisfy the integral equality (1): a weak solution of (BVP) must also satisfy its BCs. Since Prop. 3.1 imposes no restrictions on $\mathcal{T}$ to ensure this is the case, we combine the constraint (5) with a worst-case loss describing the BC satisfaction. By once again applying Prop. 3.1, we can find a weak solution of (BVP) by solving the constrained learning problem

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{BC}}}\! \Big{[}\big{(}u_{\theta}(x,t)-h(x,t)\big{)}^{2}\Big{]}$ (PIII) $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{PDE}}}\! \Big{[}\big{(}D[u_{\theta}](x,t)-\tau(x,t)\big{)}^{2}\Big{]}\leq\epsilon\text{.}$ $$

We write $\psi_{\alpha}^{\text{BC}}$ and $\psi_{\alpha}^{\text{PDE}}$ to emphasize that they are different distributions, i.e.,
$\psi_{\alpha}^{\text{BC}}(x,t)\propto\big{[}\big{(}u_{\theta}(x,t)-h(x,t)\big{
)}^{2}-\alpha\big{]}_{+}$
supported on $\mathcal{B}$
and $\psi_{\alpha}^{\text{PDE}}(x,t)\propto\big{[}\big{(}D[u_{\theta}](x,t)-\tau(x,
t)\big{)}^{2}-\alpha\big{]}_{+}$
supported on $\mathcal{D}$. The presence of these distributions is what make (PIII) considerably different from other constrained formulations such as (Lu et al., ; Basir & Senocak, ). Indeed, only by using $\psi_{\alpha}^{\text{BC}}$ and $\psi_{\alpha}^{\text{PDE}}$ can we claim that (PIII) approximates (weak) solutions of (BVP) (Prop. 3.1).

###### Proposition 3.2 .

##### Structural knowledge.

It is often possible to obtain information about the structure of the solution of a BVP, such as invariances or symmetries, without explicitly solving it (Olver, ; Akhound-Sadegh et al., ). While this structural information is already encoded in (PIII), using it explicitly can reduce training time as well as the number of both collocation points and/or observations [as in (PII)] needed. It also helps ensure the physical validity of outcomes by explicitly avoiding degenerate solutions, a common failure mode of unsupervised methods such as PINNs (Krishnapriyan et al., ) (we do not observe such issues with (PIII), see Sec. 5.1).

While it is sometimes possible to incorporate structural knowledge into the model $u_{\theta}$, e.g., using equivariant architectures (Cohen & Welling, ; Batzner et al., ), it can also be imposed as a statistical constraint. This is convenient for when such models are intricate to design. Consider, for example, a (finite) invariance group $\mathcal{G}$ whose elements $\gamma_{i}$ act on the domain $(x,t)\in\mathcal{D}$ such that $u^{\dagger}(x,t)=u^{\dagger}[\gamma_{i}(x,t)]$, where $u^{\dagger}$ is a solution of (BVP). This invariance can be enforced by

$$ $\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{ST}}}\!\Big{[}\big{(} u_{\theta}(x,t)-u_{\theta}[\gamma_{i}(x,t)]\big{)}^{2}\Big{]}\leq\epsilon\text {,}\quad\gamma_{i}\in\mathcal{G}\text{,}$ (6) $$

where $\psi_{\alpha}^{\text{ST}}$ supported on $\mathcal{D}\times\mathcal{G}$. Similar constraints can be constructed around other structures, such as equivariance. In contrast to (PIII), where we want $\epsilon\approx 0$, it can be beneficial to use a larger values in (6) to accommodate models $u_{\theta}$ that cannot fully capture the solution invariances.

##### Observational knowledge.

The last type of scientific knowledge we consider are (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 a supervised manner as in (PII), it is not limited to that architecture. Given observations $u^{\dagger}_{j}$, $j=1,\dots,J$, of solutions of (BVP), we may formulate constraints of the kind

$$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{OB}}}\! \Big{[}\big{(}u_{\theta}(x,t)-u_{j}^{\dagger}(x,t)\big{)}^{2}\Big{]}$ $\displaystyle\leq\epsilon\text{,}\quad j=1,\dots,J\text{.}$ (7) $$

If $u^{\dagger}_{j}$ is only known for a subset of $(x_{q},t_{q})$, then $\psi_{\alpha}^{\text{OB}}$ becomes a discrete distribution over those points. This is the case, for instance, of Fourier NOs (FNOs) that only make predictions on a regular grid (Li et al., ). By considering each data point individually instead of averaging them as in (PII), (7) accounts for the heterogeneous difficulty of fitting them. It also enables the tolerance $\epsilon$ to be adjusted individually for each sample, e.g., using larger values for noisier observations.

##### Science-constrained learning.

We can now combine (PIII) with (6) and (7) to formulate a constrained learning problem accommodating all knowledge sources considered so far. Explicitly,

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{BC}}}\! \Big{[}\big{(}u_{\theta}(x,t)-h(x,t)\big{)}^{2}\Big{]}$ (M) (SCL) $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{PDE}}}\! \Big{[}\big{(}D[u_{\theta}](x,t)-\tau(x,t)\big{)}^{2}\Big{]}\leq\epsilon_{ \text{pde}}$ (M) $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{ST}}}\! \Big{[}\big{(}u_{\theta}(x,t)-u_{\theta}[\gamma_{i}(x,t)]\big{)}^{2}\Big{]} \leq\epsilon_{\text{s}}\text{,}\quad\gamma_{i}\in\mathcal{G}$ (S) $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{OB}}}\! \Big{[}\big{(}u_{\theta}(x,t)-u_{j}^{\dagger}(x,t)\big{)}^{2}\Big{]}\leq \epsilon_{\text{o}}\text{,}\quad j=1,\dots,J\text{.}$ (O) $$

Note that only a subset of the constraints in (SCL) can be used depending on the information available. 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 an objective, i.e., as a feasibility problem.

Naturally, (SCL) can be extended to simultaneously fit the solutions of a parameterized family of BVPs. However, rather than discretizing the parameter space as in (PI), we rely on the robust, worst-case formulation introduced in Sec. 3.1. Explicitly, we rewrite (SCL) as

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t,\pi)\sim\psi_{\alpha}^{\text{BC}} ,\,\tau\sim\mathfrak{m}}\!\Big{[}\big{(}u_{\theta}(\pi,\tau)(x,t)-h(x,t)\big{) }^{2}\Big{]}$ (M) (SCL ′ ) $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t,\pi)\sim\psi_{\alpha}^{\text{PDE} },\,\tau\sim\mathfrak{m}}\!\Big{[}\big{(}D_{\pi}[u_{\theta}(\pi,\tau)](x,t)- \tau(x,t)\big{)}^{2}\Big{]}\leq\epsilon_{\text{pde}}$ (M) $\displaystyle\operatorname{\mathbb{E}}_{(x,t,\pi)\sim\psi_{\alpha}^{\text{ST}} ,\,\tau\sim\mathfrak{m}}\!\bigg{[}\Big{(}u_{\theta}(\pi,\tau)(x,t)-u_{\theta}( \pi,\tau)\big{[}\gamma_{i}(\pi)(x,t)\big{]}\Big{)}^{2}\bigg{]}\leq\epsilon_{ \text{s}}\text{,}\quad\gamma_{i}\in\mathcal{G}$ (S) $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{OB}}}\! \Big{[}\big{(}u_{\theta}(\pi_{j},\tau_{j})(x,t)-u_{j}^{\dagger}(x,t)\big{)}^{2 }\Big{]}\leq\epsilon_{\text{o}}\text{,}\quad j=1,\dots,J\text{,}$ (O) $$

where $\mathfrak{m}$ is a fixed distribution (e.g., uniform), $u_{j}^{\dagger}$ is a solution of (BVP) with coefficients $\pi_{j}$ and forcing function $\tau_{j}$, and the invariance $\gamma_{i}$ is now parametrized to account for the fact that its action may depend on $\pi$ (e.g., translation invariance with different strides). By incorporating the coefficients $\pi$ into the worst-case distribution $\psi_{\alpha}$, we consider all of its possible values rather than only a finite subset, achieving performances that would require fine discretizations (see Sec. 5). Additionally, since the $\psi_{\alpha}$ are now supported on $\mathcal{B}\times\Pi$ or $\mathcal{D}\times\Pi$, it targets not only difficult BVPs but also the region of the domain responsible for this difficulty. In contrast, this approach is not directly applicable to infinite-dimensional parameters such as $\tau$, since it would require sampling from a function space. We therefore leave this extension for future work, considering instead a fixed distribution $\mathfrak{m}$.

Since (SCL) and (SCL^′) are constrained learning problems, they can be solved using their empirical dual form (D̂-CSL). Doing so, however, requires samples from the distributions $\psi_{\alpha}$ whose closed-form we do not know (see Prop. 3.1). Next, we address this challenge and develop a practical algorithm for solving SCL problems.

## 4 Algorithm

The first step in developing a practical algorithm for (SCL)/(SCL^′) is to obtain samples from the different $\psi_{\alpha}$. Yet, Prop. 3.1 only describes these distributions implicitly. For instance, the PDE constraint in (5) uses $\psi_{\alpha}^{\text{PDE}}\propto\big{[}(D[u_{\theta}](x,t)-\tau(x,t))^{2}-
\alpha\big{]}_{+}$. Without access to its partition function, sampling $\psi_{\alpha}^{\text{PDE}}$ can be intricate. Additionally, the truncation threshold $\alpha$ controls not only the approximation of the robust loss, but also the discontinuity and sparsity of $\psi_{\alpha}^{\text{PDE}}$. Indeed, as the value of $\alpha$ grows, $\psi_{\alpha}^{\text{PDE}}$ displays increasingly isolated modes. Sampling from such distributions is well-known to be challenging (Nishimura et al., ; Robert & Casella, ).

On the other hand, since the losses in (SCL) are non-negative, all the $\psi_{\alpha}$ are smooth and fully supported for $\alpha=0$. They are therefore amenable to be sampled using Markov Chain Monte Carlo (MCMC) techniques (Robert & Casella, ). We use the gradient-free Metropolis-Hastings (MH) algorithm (see App. C) in our experiments to avoid additional backward passes and the higher-order derivatives resulting from differentiating $D_{\pi}$ as in Wang et al. (). Still, both $\psi_{0}$ and MH turn out to provide strong practical results (see Sec. 5). A similar pattern has been observed in other ML applications (Robey* et al., ; Hounie et al., ).
Naturally, sampling from $\psi_{0}$ remains computationally demanding. It is, however, possible to replace $\psi_{\alpha}$ by a fixed distribution (e.g., uniform) or even fixed points for some of the constraints in (SCL)/(SCL^′) at negligible accuracy costs. We find this is consistently the case, e.g., for the BC objective, although we emphasize that this is application-dependent. Additionally, certain architectures cannot make predictions on arbitrary domain points, such as FNOs that only consider regular grids (Li et al., ). In those cases, it is more appropriate to take $\psi_{\alpha}$ to be a uniform distribution over equispaced points.

Having obtained samples from the $\psi_{\alpha}$ in (SCL)/(SCL^′), we can now construct its empirical dual problem (D̂-CSL) and proceed to solving it. Since obtaining Lagrangian minimizers can be time consuming, we replace the dual ascent method from (3) by the alternating scheme in Alg. 1. Such gradient descent-ascent schemes are commonly used in convex optimization (Bertsekas, ), although their convergence is less well understood in non-convex settings Lin et al. (); Fiez et al. (); Yang et al. ().
Note that contrary to previous constrained approaches (Lu et al., ; Basir & Senocak, ), Alg. 1 does not include training heuristics, such as conditional updates.
Next, we showcase the results of SCL for a variety of applications.

Figure: Algorithm 1 Primal-dual method for (SCL)

## 5 Experiments

In this section, we showcase the use of SCL to train both MLPs and FNOs (Li et al., ) to solve six different PDE problems (convection, reaction-diffusion, eikonal, Burgers’, diffusion-sorption, and Navier-Stokes). Throughout these experiments, we consider different subsets of constraints from (SCL)/(SCL^′) to illustrate a variety of knowledge settings. We denote by (SCL)(M) the problem (SCL) with only constraints (M).
For each setting, we consider only the most suitable model since our goal is not to exhaust the potential of SCL, but to illustrate some of its natural uses. Details of the experimental setups are provided in the appendices, including descriptions of the PDEs (App. A), training procedures, hyperparameters (App. D), and additional results (App. E).
In the sequel, we use fixed points $(x,t)$ for the BC objective rather than $\psi_{0}^{\text{BC}}$ to illustrate how computational complexity can be reduced without significantly affecting the results. We still use $\psi_{0}^{\text{BC}}$ for the marginal $\pi$.

### 5.1 Solving a specific BVP

**Table 1: Relative $L_{2}$ error for solving specific BVPs.**
|  | Convection | Reaction-Diffusion | Eikonal |  |  |
| --- | --- | --- | --- | --- | --- |
|  | $\beta=30$ | $\beta=50$ | $(\nu,\rho)=(3,3)$ | $(\nu,\rho)=(3,5)$ |  |
| PINN (PI) | $0.662\,\%$ | $64.8\,\%$ | $0.755\,\%$ | $80.0\,\%$ | $106\,\%$ |
| R3 (Daw et al., ) | $0.995\,\%$ | $2.32\,\%$ | $0.854\,\%$ | $0.746\,\%$ | $100\,\%$ |
| (SCL)(M) | $0.591\,\%$ | $2.50\,\%$ | $0.578\,\%$ | $0.493\,\%$ | $7.40\,\%$ |

We begin with a traditional use-case: solving (BVP) for a fixed coefficient $\pi$, forcing function $\tau$, and BCs. We use (SCL)(M) to train MLPs to solve convection, reaction-diffusion, and eikonal problems, comparing the results with PINNs [i.e., (PI) with weights chosen as in (Daw et al., )]. All experiments use $1000$ collocation points per epoch obtained by (PINN) sampling uniformly (we find this performs better than using fixed points as in, e.g., (Raissi et al., ; Lu et al., )), (R3) using the adaptive heuristic from (Daw et al., ), and (SCL) using MH after $4000$ burn-in steps.
Table 1 shows that SCL matches and often outperforms other methods, particularly in challenging scenarios (e.g., convection with $\beta=50$ or the non-linear eikonal). This is because it jointly adapts the weight of each loss and the points used to evaluate them during training (see Alg. 1). Although R3 also targets points with high PDE residuals, it does not approximate the worst-case loss needed to guarantee a (weak) solution (Prop. 3.2). The improved performance of SCL comes at a higher computational cost ($\approx 5$x), although this is largely offset when solving parametric problems.

### 5.2 Solving parametric families of BVPs

One of the key advantages of NN-based approaches is their ability to solve entire families of BVPs without retraining. Consider fitting an MLP $u_{\theta}(x,t,\beta)$ to simultaneously solve convection problems for $\beta\in[1,30]$ using (SCL^′)(M) and (PI) as in Cho et al. (). Note that we do not use their tailored “P^2INN” architecture, although it would be compatible with SCL. The $\pi_{j}=\beta_{j}$ are taken to be 4, 7, and 30 equispaced values in $[1,30]$. The results in Fig. 1a show that (PI) can only reach the error of (SCL^′)(M) for the finest discretization, at which point it evaluates the PDE loss 6 times more per epoch. The same pattern holds for reaction-diffusion (4–7 times) and 2D Helmholtz (4–5 times) problems (App. E).
The effectiveness of (SCL^′)(M) is due to its worst-case loss jointly selecting $(x,t,\beta)$ to target both challenging coefficients and the domain region responsible for that difficulty. In fact, by inspecting $\psi_{0}^{\text{PDE}}$ at different stages of training (Fig. 1b) it becomes clear that SCL begins by fitting the solution “causally,” focusing first on smaller values of $t$. While doing so has been proposed to improve training (Krishnapriyan et al., ; Wang et al., ), it arises naturally by solving (SCL^′)(M). As training advances, however, Alg. 1 shifts focus to fitting higher convection speeds $\beta$. Note that this occurs *without* any prior knowledge of the problems or manual tuning.

Figure: Figure 1: Solving parametric convection BVPs: (a) relative $L_{2}$ error vs. $\beta$ (legend reports number of differential operator evaluations per epoch). (b) Samples from $\psi_{0}^{\text{PDE}}$ at different training stages.
Refer to caption: x1.png

### 5.3 Leveraging invariance when solving BVPs

Next, we showcase how SCL can be used to overcome computational limitations or scarce mechanistic knowledge. Consider training an MLP to solve a convection BVP with $\beta=30$ and periodic initial condition $h(0,x)=\sin(x)$. This problem is commonly used to showcase a “failure mode” of (PI) since it yields degenerate solutions when using fixed collocation points (Fig. 2a) (Krishnapriyan et al., ).
Note that the constrained (SCL)(M) also fails when using fixed collocations points (Fig. 2b), even though its stochastic version finds accurate solutions (Table 1).
One way to overcome this limitation is by leveraging additional knowledge. In this case, we know from the BVP structure that its solution must be periodic with period $\pi/15$. By incorporating this invariance using (SCL)(S), we can obtain accurate solutions despite our use of fixed collocation points for (SCL)(M) (Fig. 2c). Note that the worst-case loss in (SCL)(S) is fundamental to avoid degenerate solutions.

Figure: Figure 2: Using invariance in convection BVPs (BC and PDE losses in (PI) and (SCL) are evaluated using a *fixed set of collocation points*).
Refer to caption: x3.png

### 5.4 Supervised solution of BVPs

We conclude by exploring supervised methods that rely only on pairs of initial conditions $h_{j}(x,0)$ and corresponding (weak) solutions $u^{\dagger}_{j}$. To showcase the versatility of SCL, we train FNOs rather than MLPs, fixing $\psi_{\alpha}^{\text{OB}}$ to a uniform distribution over a regular grid to accommodate their predictions. We also cast the SCL problem using only the observational constraints (SCL)(O), i.e., without any objective. In contrast to (PII), which minimizes the average error, this formulation enforces a maximum error of $\epsilon_{\text{o}}$ across samples. Though apparently minor, this leads to better prediction quality, as shown in 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 a maximum error gives more weight to challenging data points. This becomes clear in Fig. 3, which compares the magnitude of each IC in the training set with its final dual variables $\lambda^{\text{OB${}_{j}$}}$ 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 to guide the collection of additional data points or improve the NO architecture. Indeed, any model improvement can immediately take advantage of SCL since it is (virtually) independent of the choice of$u_{\theta}$. The benefits of $\lambda^{\text{OB}}$ are exacerbated when predicting which IC is hard to fit is intricate, as is the case for the Navier-Stokes equation (see App. E).

Figure: Table 2: Relative $L_{2}$ error on test set (FNO).
Refer to caption: x4.png

## 6 Conclusion

This paper developed SCL, a technique for solving BVPs based on constrained learning. It demonstrated that finding (weak) solutions of PDEs is equivalent to solving constrained learning problems with worst-case losses. It also showed how additional prior knowledge can be incorporated to the solution of the BVP, such as structural constraints (e.g., invariances), real-world measurements, and previously known solutions. It then developed a practical algorithm to tackle these problems and showcased its performance across a variety of PDEs, NN architectures, and prior knowledge levels. SCL not only yields accurate BVP solutions, but tackles many of the challenges faced by previous methods, such as extensive hyperparameter tuning, degenerate solutions, and fine discretizations of domain and/or coefficients. Future work includes exploring structures beyond invariance, worst-case losses for functional parameters such as $\tau$, and active learning methods based on dual variables.

#### Acknowledgments

This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (EXC 2075-390740016). This work 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. The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS)
for supporting Viggo Moro.

## Appendix A Applications of ( BVP )

We showcase here a variety of phenomena that can be described using (BVP). The PDEs detailed in this section are used throughout our experiments to illustrate the performance of SCL.

### A.1 Convection equation

The one-dimensional convection equation models the transport of a scalar quantity $u(x,t)$, such as temperature or concentration, along the spatial dimension $x$. We consider the convection problem

$$ $\displaystyle\frac{\partial u(x,t)}{\partial t}+\beta\frac{\partial u(x,t)}{ \partial x}$ $\displaystyle=0,$ $\displaystyle(x,t)$ $\displaystyle\in(0,2\pi)\times(0,1]$ (8a) $\displaystyle u(x,0)$ $\displaystyle=\sin(x),$ $\displaystyle x$ $\displaystyle\in[0,2\pi]$ (8b) $\displaystyle u(0,t)$ $\displaystyle=u(2\pi,t),$ $\displaystyle t$ $\displaystyle\in(0,1]$ (8c) $$

where the coefficient $\beta$ denotes the *convection speed*. Despite its use of periodic BCs, namely, (8c), this problem (and its variations) can be cast as a straightforward extension of (BVP).

### A.2 Reaction-diffusion equation

The one-dimensional reaction-diffusion equation describes a variety of phenomena, including chemical 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,

$$ $\displaystyle\frac{\partial u(x,t)}{\partial t}-\nu\frac{\partial^{2}u(x,t)}{ \partial x^{2}}$ $\displaystyle=\rho u(x,t)\big{(}1-u(x,t)\big{)},$ $\displaystyle(x,t)$ $\displaystyle\in(0,2\pi)\times(0,1]$ (9a) $\displaystyle u(x,0)$ $\displaystyle=\exp\bigg{(}-\frac{1}{2}\Big{(}\frac{x-\pi}{\pi/4}\Big{)}^{2} \bigg{)},$ $\displaystyle x$ $\displaystyle\in[0,2\pi]$ (9b) $\displaystyle u(0,t)$ $\displaystyle=u(2\pi,t),$ $\displaystyle t$ $\displaystyle\in(0,1]$ (9c) $$

where $\nu>0$ and $\rho$ are the diffusion and reaction coefficients, respectively.

### A.3 Eikonal equation

The eikonal equation is encountered in many applications involving wave propagation, e.g., electromagnetism. It also describes the (signed) distance between any point $x\in\Omega$ and some fixed boundary $\partial\mathcal{S}$, hence its usage in vision applications. In this case, we consider the BVP

$$ $\displaystyle\left\|\nabla u(x,y)\right\|$ $\displaystyle=1,\quad$ $\displaystyle(x,y)$ $\displaystyle\in\Omega$ (10a) $\displaystyle u(x,y)$ $\displaystyle=0,$ $\displaystyle(x,y)$ $\displaystyle\in\partial\mathcal{S}$ (10b) $$

where $\Omega=(-1,1)^{2}$ and $\partial\mathcal{S}$ is a complex shape, in our case, the gears figure from (Daw et al., ). In order to ensure that negative (positive) distances are assigned to the interior (exterior) of the shape, we add a structural constraint to the solution which enforces that $u$ is non-negative on the boundary of $\Omega$. Explicitly, we enforce that

$$ $u(x,y)\geq 0\text{,}\quad(x,y)\in\partial\Omega\text{.}$ (11) $$

This can be done by using different losses for the IC in (10b) and the BC (11), replacing the objective of (SCL) by

$$ $\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{\alpha}^{\text{IC}}}\!\Big{[}\big{(} u_{\theta}(x,t)-h(x,t)\big{)}^{2}\Big{]}+\operatorname{\mathbb{E}}_{(x,y)\sim \psi_{\alpha}^{\text{BC}}}\!\Big{[}\big{[}-u_{\theta}(x,y)\big{]}_{+}\Big{]}$ (12) $$

for $\psi_{\alpha}^{\text{IC}}$ supported on $\partial\mathcal{S}$ and $\psi_{\alpha}^{\text{BC}}$ supported on $\partial\Omega$.

### A.4 Helmholtz equation

The two-dimensional Helmholtz equation models wave propagation and vibration phenomena in various physical contexts. Here, we consider the problem

$$ $\displaystyle\nabla^{2}u(x,y)+k^{2}u(x,y)$ $\displaystyle=\tau(x,y,\pi),$ $\displaystyle(x,y)$ $\displaystyle\in\text{interior}(\mathcal{S})$ (13a) $\displaystyle u(x,y)$ $\displaystyle=\sin(\pi a_{1}x)\sin(\pi a_{2}y),$ $\displaystyle(x,y)$ $\displaystyle\in\partial\mathcal{S}$ (13b) $$

where $\mathcal{S}=[0,1]^{2}$; $k>0$ is the wave number; coefficients $\pi=(a_{1},a_{2})$ represent the spatial frequencies in the $x$ and $y$ directions, respectively; and the forcing function is $\tau(x,y,\pi)=\left(k^{2}-\pi^{2}a_{1}^{2}-\pi^{2}a_{2}^{2}\right)\left(\sin(
\pi a_{1}x)\sin(\pi a_{2}y)\right)$.

### A.5 Burgers’ equation

The one-dimensional Burgers’ equation models the behavior of a scalar field $u(x,t)$ under the combined effects of nonlinear convection and diffusion. It 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

$$ $\displaystyle\frac{\partial u(x,t)}{\partial t}+\dfrac{1}{2}\frac{\partial u^{ 2}(x,t)}{\partial x}$ $\displaystyle=\nu\frac{\partial^{2}u(x,t)}{\partial x^{2}},$ $\displaystyle(x,t)$ $\displaystyle\in(0,1)\times(0,1]$ (14a) $\displaystyle u(x,0)$ $\displaystyle=h_{0}(x),$ $\displaystyle x$ $\displaystyle\in[0,1]$ (14b) $\displaystyle u(0,t)$ $\displaystyle=u(1,t),$ $\displaystyle t$ $\displaystyle\in(0,1]$ (14c) $$

where $\nu>0$ is the viscosity coefficient, which governs the strength of the diffusion term. We consider ICs $h_{0}$ from the same distribution as Li et al. ().

### A.6 Diffusion-sorption equation

The diffusion-sorption equation models the transport of a scalar field $u(x,t)$ (e.g., a contaminant concentration) in a porous medium. It is commonly used in environmental science and chemical engineering. In terms of this PDE, we consider the following BVP

$$ $\displaystyle\frac{\partial u(x,t)}{\partial t}$ $\displaystyle=\frac{\nu}{R\big{(}u(x,t)\big{)}}\frac{\partial^{2}u(x,t)}{ \partial x^{2}},$ $\displaystyle(x,t)$ $\displaystyle\in(0,1)\times(0,500]$ (15a) $\displaystyle u(x,0)$ $\displaystyle=h_{0}(x),$ $\displaystyle x$ $\displaystyle\in[0,1]$ (15b) $\displaystyle u(0,t)$ $\displaystyle=1\text{,}\quad u(1,t)=\nu\frac{\partial u(1,t)}{\partial x}\text {,}$ $\displaystyle t$ $\displaystyle\in(0,500]$ (15c) $$

where $\nu=5\times 10^{-4}$ is the diffusion coefficient and $R$ is the retardation factor defined as

$$ $R(u)=1+\frac{(1-\phi)}{\phi}\rho_{s}kn_{f}u^{n_{f}-1}\text{,}$ $$

where $\rho_{s}=2880$ is the bulk density, $\phi=0.29$ is the medium porosity, $k=3.5\times 10^{-4}$ is Freundlich’s sorption parameter, and $n_{f}=0.874$ is Freundlich’s exponent. We consider the distribution of ICs $h_{0}$ from Takamoto et al. ().

A notable aspect of (15) is its nonlinearity due to the dependence on $u(x,t)$ of the effective diffusion coefficient $\frac{\nu}{R(u)}$. What is more, it can become singular when $u=0$, making this equation particularly challenging to solve.

### A.7 Navier-Stokes equation

The two-dimensional, incompressible Navier-Stokes equation in vorticity form removes the pressure term and focuses on the dynamics of rotational flow. It is used to describe the local rotation of a fluid, i.e., the vorticity $\omega(x,t)=\nabla\times u(x,t)$, where $u$ is the two-dimensional velocity field and $\nabla\times f$ denotes the *curl* of $f$. Although vorticity is more challenging to model than velocity, it offers a deeper understanding of the flow dynamics. The velocity field can be recovered from the vorticity using Poisson’s equation.

Explicitly, for $x=[x_{1},x_{2}]^{\top}$, we consider the BVP

$$ $\displaystyle\frac{\partial\omega(x,t)}{\partial t}+u(x,t)^{\top}\nabla\omega( x,t)$ $\displaystyle=\nu\nabla^{2}\omega+\tau(x),\quad$ $\displaystyle(x,t)$ $\displaystyle\in(0,1)^{2}\times(0,T]$ (16a) $\displaystyle\omega(x,0)$ $\displaystyle=\omega_{0},\quad$ $\displaystyle x$ $\displaystyle\in[0,1]^{2},$ (16b) $\displaystyle\nabla\cdot u(x,t)$ $\displaystyle=0,\quad$ $\displaystyle(x,t)$ $\displaystyle\in(0,1)^{2}\times(0,T]$ (16c) $$

where $\nu>0$ is the viscosity coefficient and $\nabla\cdot f$ denotes the divergence of $f$. The forcing function is taken to be

$$ $\tau(x)=0.1\big{(}\sin\big{(}2\pi(x_{1}+x_{2})\big{)}+\cos\big{(}2\pi(x_{1}+x_ {2})\big{)}\big{)}$ $$

and the ICs $\omega_{0}$ are taken from the same distribution as in Li et al. (). We consider three settings, namely, $\nu=10^{-3}$ with $T=50$, $\nu=10^{-4}$ with $T=30$, and $\nu=10^{-5}$ and $T=20$.

## Appendix B Weak solutions and robust learning

The space $\mathcal{T}$ of test functions plays a fundamental role when defining the weak formulations (1). Typically, it is chosen to be a Sobolev space due to its natural compatibility with BVPs and the fact that it leads to less stringent differentiability requirements on $u$. It therefore overcomes the main issues with the strong formulation (BVP). A Sobolev space consists of functions in some Lebesgue space $L^{p}$ whose *weak derivatives* are also in $L^{p}$. Recall that $L^{p}$ is the space of $p$-integrable functions. To define a Sobolev space, we therefore need to start by defining a weak derivative.

A locally integrable function $f$ defined on an open set $\mathcal{Z}\subset\mathbb{R}^{d}$ is weakly differentiable with respect to $z_{i}$ if there exists $Df$ also locally integrable (i.e., $f,Df\in L^{1}_{\text{loc}}(\mathcal{Z})$) such that

$$ $\int_{\mathcal{Z}}Df(z)\xi(z)dz=-\int_{\mathcal{D}}\varphi(z)\dfrac{\partial \xi(z)}{\partial z_{i}}dz,\quad\text{for all }\xi\in C_{c}^{\infty}(\mathcal{Z })\text{,}$ (17) $$

where $C_{c}^{\infty}(\mathcal{Z})$ is the space of infinitely differentiable, compactly supported functions. We say $Df$ is the weak derivative of $f$. To generalize (17) to higher-order derivatives, consider the multi-index $\alpha=(\alpha_{1},\ldots,\alpha_{d})$ to be a $d$-tuple of non-negative integers and let ${\left|\alpha\right|}=\sum_{i=1}^{d}\alpha_{i}$. We then define the $\alpha$-weak derivatives of $f$, denoted $D^{\alpha}f$, as

$$ $\int_{\mathcal{Z}}D^{\alpha}f(z)\xi(z)dz=-\int_{\mathcal{D}}\varphi(z)\frac{ \partial^{{\left|\alpha\right|}}\xi(z)}{\partial z_{1}^{\alpha_{1}}\cdots \partial z_{d}^{\alpha_{d}}}dz,\quad\text{for all }\xi\in C_{c}^{\infty}( \mathcal{Z})\text{.}$ (18) $$

We can now define what we mean by Sobolev space (Evans, ).

###### Definition B.1 (Sobolev space) .

For an integer $k\geq 0$ and $1\leq p\leq\infty$, we define the Sobolev space $W^{k,p}(\mathcal{Z})=\{f\in L^{p}(\mathcal{Z})\mid D^{\alpha}f\in L^{p}(
\mathcal{Z})\text{ for all multi-indices }\alpha\text{ with }{\left|\alpha
\right|}\leq k\}$.

We write $W^{k,p}$ whenever the set $\mathcal{Z}$ is clear from the context. Note that since $L^{p}\subset L^{1}_{\text{loc}}$ for $p\geq 1$, Sobolev spaces impose more restrictions than weak differentiability. Also, while $W^{k,p}$ is in general a Banach space, $W^{k,2}$ is a Hilbert space (Evans, ).

Having set the groundwork, we can now proceed with the proof of Prop. 3.2.

###### Definition B.1 (Sobolev space) .

### B.1 Proof of Prop. 3.2

The proof follows by constructing a measure of the deviation from the weak formulation (1) and showing that it is dominated by the proposed worst-case statistical loss. This immediately implies Prop. (3.2). Explicitly, note that the weak formulation (1) can equivalently be expressed as

$$ $\left[\int_{\mathcal{D}}\Big{(}D[u](x,t)-\tau(x,t)\Big{)}\varphi(x,t)dxdt \right]^{2}=0\text{,}\quad\text{for all }\varphi\in W^{k^{\prime},2}$ (19) $$

where we omitted the dependence on the coefficients $\pi$ for conciseness.
Using Jensen’s inequality, we can upper bound (19) for any $\varphi$ as in

$$ $\left[\int_{\mathcal{D}}\Big{(}D[u](x,t)-\tau(x,t)\Big{)}\varphi(x,t)dxdt \right]^{2}\leq\int_{\mathcal{D}}\Big{(}D[u](x,t)-\tau(x,t)\Big{)}^{2}\varphi^ {2}(x,t)dxdt\text{.}$ (20) $$

Since $\varphi\in W^{k^{\prime},2}\subset L^{2}$, the partition function $Z_{\psi}=\int_{\mathcal{D}}\varphi(x,t)^{2}dxdt=\left\|\varphi\right\|_{L^{2}}
^{2}<\infty$ is well-defined. We can therefore consider the normalized $\psi=\varphi^{2}/Z_{\psi}$ in (20) to get

$$ $\left[\int_{\mathcal{D}}\Big{(}D[u](x,t)-\tau(x,t)\Big{)}\varphi(x,t)dxdt \right]^{2}\leq Z_{\psi}\int_{\mathcal{D}}\Big{(}D[u](x,t)-\tau(x,t)\Big{)}^{2 }\psi(x,t)dxdt\text{.}$ (21) $$

Since $\psi$ is a non-negative, normalized function, it is the density of a probability measure, i.e., $\psi\in\mathcal{P}$. The following technical lemma shows that it is in fact in $\mathcal{P}^{2}$, i.e., it is a square-integrable probability density.

###### Lemma B.2 .

Let $\varphi\in W^{k^{\prime},2}$ and $\psi\propto\varphi^{2}$ be a probability distribution. Then, $\psi\in\mathcal{P}^{2}$.

Before proving Lemma B.2, let us conclude the proof. The hypothesis on $u^{\dagger}$ implies that

$$ $\sup_{\psi\in\mathcal{P}^{2}}\ \int_{\mathcal{D}}\Big{(}D[u^{\dagger}](x,t)- \tau(x,t)\Big{)}^{2}\psi(x,t)dxdt=0$ $$

and since $Z_{\psi}$ is bounded, we obtain from (21) that

$$ $\left[\int_{\mathcal{D}}\Big{(}D[u^{\dagger}](x,t)-\tau(x,t)\Big{)}\varphi(x,t )dxdt\right]^{2}\leq 0\text{.}$ $$

Noticing that this holds for all $\psi\in W^{k^{\prime},2}$, we recover (19), which concludes the proof. $\blacksquare$

###### Proof of Lemma B.2 .

To show $\psi\in\mathcal{P}^{2}$, we must show that $\phi\in L^{4}$. Indeed,

$$ $\int_{\mathcal{D}}\psi(x,t)^{2}dxdt=\frac{1}{Z^{2}}\int_{\mathcal{D}}\varphi(x ,t)^{4}dxdt=\left(\frac{\left\|\varphi\right\|_{L^{4}}}{\left\|\varphi\right\| _{L^{2}}}\right)^{4}\text{.}$ (22) $$

Since $\varphi\in W^{k^{\prime},2}\subset L^{2}$, suffices it to show that the numerator is finite. To do so, we can use the Gagliardo-Nirenberg interpolation inequality (Nirenberg, ; Fiorenza et al., ) to write

$$ $\left\|\varphi\right\|_{L^{4}}\leq C\big{(}\left\|D^{m}\varphi\right\|_{L^{2}} ^{\alpha}\left\|\varphi\right\|_{L^{2}}^{1-\alpha}+\left\|\varphi\right\|_{L^{ 2}}\big{)}\text{,}\quad\text{for }\alpha\in[0,1]\text{,}$ (23) $$

as long as (i) $k^{\prime}\geq m\in\mathbb{N}$ and (ii) $4\alpha m=d+1$. An additional condition must hold in particular cases:

$$ $\text{(iii) if }m-\dfrac{d+1}{2}\text{ is a non-negative integer, then }\alpha <1\text{.}$ $$

Note that $d+1$ is the dimension of the space-time domain $\mathcal{D}$. Since $\varphi\in W^{k^{\prime},2}$, the right-hand of (23) is finite. We therefore only need to show that (i)–(iii) are satisfied under the hypothesis of the theorem. To do so, we consider three cases:

- (a)
$\bm{d=0}$: in this case, $d+1$ is odd so that condition (iii) does not apply. Immediately, (23) holds for $m=1$ and $\alpha=0.25$.
- (b)
$\bm{d=1}$: since (ii) requires $m\geq 1$, we now have that $m-\frac{d+1}{2}$ is a non-negative integer. Hence, condition (iii) applies. Nevertheless, we can still take $m=1$ and $\alpha=0.5$ in (23).
- (c)
$\bm{d\geq 2}$: we can now consider all other cases by taking
$m=\left\lceil\dfrac{d+1}{4}\right\rceil\quad\text{and}\quad\alpha=\frac{d+1}{4
m}\text{.}$
Indeed, (ii) holds by construction. Additionally, from the hypothesis of the theorem, we have $d+1\leq 4k^{\prime}$, which by the monotonicity of the ceiling operation implies that (i) holds. It suffices to show that (iii) never applies. For $2\leq d\leq 3$, we have $m=1$ and $(d+1)/2>1$. For $d>3$, we have
$\left\lceil\dfrac{d+1}{4}\right\rceil-\dfrac{d+1}{2}\leq\left(\dfrac{d+1}{4}+1
\right)-\dfrac{d+1}{2}<0\text{.}$
(24)

Thus, under the hypothesis of the theorem, (23) implies that $\left\|\varphi\right\|_{L^{4}}<\infty$. From (22), this in turn implies that $\psi\in\mathcal{P}^{2}$.
∎

###### Lemma B.2 .

###### Proof of Lemma B.2 .

## Appendix C Sampling with the Metropolis-Hastings algorithm

Algorithm 1 uses samples from one or more $\psi_{0}$ in order to compute the losses in steps 4–7. This is, however, not straightforward unless those distributions are fixed to, e.g., a uniform (as we typically do for the BCs). That is because we only know $\psi_{0}$ up to a normalization factor. We overcome this issue using MCMC techniques, more specifically, the Metropolis-Hastings algorithm (Robert & Casella, ). In Alg. 2, we consider the general case of sampling from $\psi_{0}\propto\ell$, where $\ell$ is a non-negative, scalar-valued loss. We denote by $z_{n}$ the desired samples and $\mathcal{R}$ their support. In Alg. 1, for instance, we would take $z_{n}=(x^{\text{pde}}_{n},t^{\text{pde}}_{n},\pi^{\text{pde}}_{n})$, $\mathcal{R}=\mathcal{D}\times\Pi$, and

$$ $\ell(z_{n})=\Big{(}D_{\pi^{\text{pde}}_{n}}\left[u_{\theta_{k}}(\pi^{\text{pde }}_{n})\right](x^{\text{pde}}_{n},t^{\text{pde}}_{n})-\tau(x^{\text{pde}}_{n}, t^{\text{pde}}_{n})\Big{)}^{2}$ (25) $$

in step 5 and $z_{n}=(x^{\text{o${{}_{j}}$}}_{n},t^{\text{o${{}_{j}}$}}_{n})$, $\mathcal{R}=\mathcal{D}$, and

$$ $\ell(z_{n})=\Big{(}u_{\theta_{k}}(\pi_{j},\tau_{j})(x^{\text{o${}_{j}$}}_{n},t ^{\text{o${}_{j}$}}_{n})-u_{j}^{\dagger}(x^{\text{o${}_{j}$}}_{n},t^{\text{o${ }_{j}$}}_{n})\Big{)}^{2}$ $$

in step 7.

Typically, the covariance $\Sigma$ is taken to be diagonal (independent proposals), e.g., $\sigma^{2}I$. The choice of parameter $\sigma^{2}$ of the proposal (step 3) affects the *mixing rate*, i.e., how fast the samples converge to the desired distribution. Smaller values of $\sigma^{2}$ 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 $\sigma^{2}$ is adapted during a burn-in phase to hit a specific acceptance rate, around $30\%$ as a rule-of-thumb (Robert & Casella, ). In our experiments, we find a reasonable value for $\sigma^{2}$ and keep it fixed throughout training. We also consider a burn-in period by using only the last $N_{0}$ samples generated by Alg. 2. For the single BVP experiments in Sec. 5.1 we use $N_{0}=1000$ and for the parameterized BVPs in Sec. 5.2 we use $N_{0}=2500$.

Note that there are better sampling methods than Algorithm 2, i.e., methods with faster mixing rates. That is the case, for instance, of Langevin Monte Carlo (LMC) (Robert & Casella, ). Yet, the LMC algorithm uses first-order information of $\ell$. For the PDE loss in (25), this means higher-order space-time derivatives of $u_{\theta}$ 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 Algorithm 2.

Figure: Algorithm 2 Metropolis-Hastings algorithm with Gaussian proposal

## Appendix D Experimental details

### D.1 Hyperparameters and implementation details

Throughout our experiments, we use the relative $L_{2}$ error as a performance metric, which we define as

$$ $e_{\text{rel}}(\pi,h)=\sqrt{\dfrac{\sum_{n=1}^{N}\big{[}u_{\theta}(\pi,h)(x_{n },t_{n})-u^{\dagger}(\pi,h)(x_{n},t_{n})\big{]}^{2}}{\sum_{n=1}^{N}\big{[}u^{ \dagger}(\pi,h)(x_{n},t_{n})\big{]}^{2}}}\text{,}$ (26) $$

where $u^{\dagger}$ is the solution of (BVP) obtained either analytically or by using classical numerical methods. For MLPs, the collocation points $\{x_{n},t_{n}\}$ are taken from a dense regular grid of points (see exact numbers below), and for FNOs, they are determined by the test sets from (Li et al., ; Takamoto et al., ). For parametrized problems, we report the average error

$$ $\bar{e}_{\text{rel}}=\dfrac{1}{J}\sum_{j=1}^{J}e_{\text{rel}}(\pi_{j},h_{j}) \text{,}$ $$

evaluated either on a dense regular grid of points (for coefficients $\pi$, see exact numbers below) or based on the test sets from (Li et al., ; Takamoto et al., ).

#### D.1.1 Solving a specific BVP (Sec. 5.1 )

In this section, we formulated SCL problems of the form (SCL)(M) in order to use the same information that PINNs traditionally rely on. Recall that we do use the worst-case distribution $\psi_{0}$ (or even random points) for the BCs, but instead consider fixed, regularly distributed points. This reduces the overall computational complexity of the problem at essentially no performance cost. Explicitly, we consider the following problem

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\frac{1}{N}\sum_{n=1}^{N}\Big{(}u_{\theta}(x^{\text{bc}}_{n},t^{ \text{bc}}_{n})-h(x^{\text{bc}}_{n},t^{\text{bc}}_{n})\Big{)}^{2}$ $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{0}^{\text{PDE}}}\!\Big{ [}\big{(}D[u_{\theta}](x,t)-\tau(x,t)\big{)}^{2}\Big{]}\leq\epsilon_{\text{pde }}\text{.}$ $$

To compute the objective for the convection and reaction-diffusion PDEs, we use $256$ points $(x^{\text{bc}},0)$, $x^{\text{bc}}\in[0,2\pi]$, for the IC and $100$ points equally spaced in $t\in(0,1]$ to evaluate the period BC. For the eikonal PDE, recall from (12) that the BC and IC use different losses. In this case, we therefore formulate the SCL problem

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\frac{1}{M}\sum_{m=1}^{M}\big{[}u_{\theta}(x_{m},y_{m})\big{]}^{2 }+\frac{1}{N}\sum_{n=1}^{N}\big{[}-u_{\theta}(x_{n},y_{n})\big{]}_{+}$ $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{0}^{\text{PDE}}}\!\Big{ [}\big{(}D[u_{\theta}](x,t)-\tau(x,t)\big{)}^{2}\Big{]}\leq\epsilon_{\text{pde }}\text{.}$ $$

where we once again consider fixed collocation points for all BCs, namely, $M=2234$ points on $\partial\mathcal{S}$ (the gears figure from (Daw et al., )) and $N=40$ points on $\partial\Omega$. We use the exact same points for (PI).

##### Problem hyperparameters.

For SCL, the tolerance $\epsilon_{\text{pde}}$ was selected by starting with a small value (e.g., $10^{-4}$) 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 $\mu$ in (PI) used in (Daw et al., ). Note that we used different weights for the BC and IC to solve the eikonal PDE, since in this case the BCs play a less critical role. All values are displayed in Table 3.

**Table 3: Problem hyperparameters for solving a specific BVP**
|  | $\mu_{\text{D}}$ | $\mu_{\text{BC}}$ | $\mu_{\text{IC}}$ | $\epsilon_{\text{pde}}$ |
| --- | --- | --- | --- | --- |
| Convection: $\beta=30$ | 1 | 100 | 100 | $10^{-3}$ |
| Convection: $\beta=50$ | 1 | 100 | 100 | $5\times 10^{-3}$ |
| Reaction-diffusion: $(\nu,\rho)=(3,3)$ | 1 | 100 | 100 | $10^{-2}$ |
| Reaction-diffusion: $(\nu,\rho)=(3,5)$ | 1 | 100 | 100 | $5\times 10^{-3}$ |
| Eikonal | 1 | 10 | 500 | $5\times 10^{-1}$ |

##### Model.

We used MLPs with 4 hidden layers for $u_{\theta}$ each with 50 neurons for the convection and reaction-diffusion equations or 128 neurons for the eikonal equation and hyperbolic tangent activation function.

##### 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 et al., ) (R3), or using Alg. 2 (SCL). For R3, we use the hyperparameters from (Daw et al., ). For Alg. 2, we use $\Sigma=\operatorname{diag}(0.25,0.01)$ for drawing proposals for $x$ and $t$ respectively for both convection and reaction-diffusion. For the eikonal PDE, we use $\Sigma=0.04\times I$. In both cases, we run the algorithm for $N=5000$ and use only the last $1000$ samples. All methods were trained using Adam with the default parameters from (Kingma & Ba, ) and learning rates described in Table 4. Note that the baselines only use learning rate $\eta_{p}$, since they do not use dual methods.

**Table 4: Training hyperparameters for solving a specific BVP**
|  | $\eta_{p}$ | $\eta_{d}$(only SCL) | Learning rate decay | Iterations |
| --- | --- | --- | --- | --- |
| Convection: $\beta=30$ | $10^{-3}$ | $10^{-4}$ | $0.9\eta$ every $5\,000$ iter. | $175\,000$ |
| Convection: $\beta=50$ | $10^{-3}$ | $10^{-4}$ | $0.9\eta$ every $5\,000$ iter. | $175\,000$ |
| Reaction-diffusion: $(\nu,\rho)=(3,3)$ | $10^{-3}$ | $10^{-4}$ | — | $200\,000$ |
| Reaction-diffusion: $(\nu,\rho)=(3,5)$ | $10^{-3}$ | $10^{-4}$ | — | $200\,000$ |
| Eikonal | $10^{-3}$ | $10^{-5}$ | $0.9\eta$ every $5\,000$ iter. | $50\,000$ |

##### Testing.

The solution of the convection and reaction-diffusion PDEs were tested on a dense regular grid of $256\times 100$ points $(x,t)\in\mathcal{D}$ against their analytical solutions. The solution of the eikonal PDE was tested on a dense regular grid of $384\times 384$ points $(x,y)\in\Omega$ against the ground truth predictions from (Daw et al., ).

#### D.1.2 Solving parametric families of BVPs (Sec. 5.2 )

The SCL problem we formulate here is similar to the previous section, although we used the parameterized version (SCL^′)(M). Once again, we replace the $(x,t)$ marginals of the worst-case distribution $\psi_{0}^{\text{BC}}$ by a fixed, uniform distribution. Note, however, that we keep the worst-case formulation for the coefficients $\pi$. Explicitly, we consider the SCL problem

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\operatorname{\mathbb{E}}_{\pi\sim\psi_{0}^{\text{BC}}}\!\bigg{[} \frac{1}{N_{\text{bc}}}\sum_{n=1}^{N_{\text{bc}}}\Big{(}u_{\theta}(\pi)(x^{ \text{bc}}_{n},t^{\text{bc}}_{n})-h(\pi)(x^{\text{bc}}_{n},t^{\text{bc}}_{n}) \Big{)}^{2}\bigg{]}$ $\displaystyle\operatorname{subject\ to}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t,\pi)\sim\psi_{0}^{\text{PDE}}}\! \Big{[}\big{(}D_{\pi}[u_{\theta}(\pi)](x,t)-\tau(\pi)(x,t)\big{)}^{2}\Big{]} \leq\epsilon_{\text{pde}}$ $$

Once again, we compute the objective for the convection and reaction-diffusion PDEs using $256$ points $(x,0)$, $x\in[0,2\pi]$, for the IC and $100$ points equally spaced in $t\in(0,1]$ to evaluate the period BC. For the Helmholtz PDE, we use $4\times 256$ points equally space around $\partial\Omega$ to evaluate the BC. We use the exact same points for (PI). Note that we include the coefficients $\pi$ in the forcing function to account for the Helmholtz BVP (see Sec. A).

##### Problem hyperparameters.

Once again, the tolerance $\epsilon_{\text{pde}}$ were selected by starting with a small value (e.g., $10^{-4}$) and increasing when the dual variables achieved too large a value during training to accommodate difficult problems. The weights $\mu$ in (PI) for the baselines were taken from (Daw et al., ). Exact values are displayed in Table 5.

**Table 5: Problem hyperparameters for solving a parametric family of BVPs**
|  | $\mu_{\text{D}}$ | $\mu_{\text{BC}}$ | $\epsilon_{\text{pde}}$ |
| --- | --- | --- | --- |
| Convection | 1 | 100 | $10^{-3}$ |
| Reaction-diffusion: $(\nu,\rho)\in[0,5]^{2}$ | 1 | 100 | $5\times 10^{-3}$ |
| Reaction-diffusion: $(\nu,\rho)\in[0,10]^{2}$ | 1 | 100 | $10^{-2}$ |
| Reaction-diffusion: $(\nu,\rho)\in[0,20]^{2}$ | 1 | 100 | $10^{-1}$ |
| Helmholtz: $(a_{1},a_{2})\in[1,2]^{2}$ | 1 | 100 | $5\times 10^{-1}$ |
| Helmholtz: $(a_{1},a_{2})\in[1,3]^{2}$ | 1 | 100 | $5$ |

##### Model.

We used MLPs with 4 hidden layers of 50 neurons.

##### Training.

When training using (PI), we used $1000$ collocation points per coefficient value, sampled uniformly at random at the beginning of each epoch. For (SCL^′)(M), we used Alg. 2. For the PDE loss, we used $\Sigma=\operatorname{diag}(0.25,0.01,\sigma_{\pi}^{2})$ to sample from $(x,t,\pi)$ for both the convection ($\sigma_{p}i^{2}=9$ for coefficient $\beta$) and reaction-diffusion [$\sigma_{p}i^{2}=(1,1)$ for coefficients $(\nu,\rho)$] equations. For the Helmholtz PDE, we used $\Sigma=0.04\times I$ to sample from $(x,t,\pi)$ for coefficients $\pi=(a_{1},a_{2})$. The same variances $\sigma_{\pi}^{2}$ were used to 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 $(a_{1},a_{2})\in[1,3]^{2}$, where we use 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, ) and learning rate of $10^{-3}$ for (PI). For SCL, the dual learning rate was $\eta_{d}=10^{-4}$. In all cases, we decayed the learning rates by a factor of $0.9$ every 5000 epochs, with the exception of the reaction-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\times 100\times 1000$ points $(x,t,\pi)\in\mathcal{D}\times\Pi$ and $256\times 100\times 100\times 100$ points $(x,t,\nu,\rho)\in\mathcal{D}\times\Pi$, respectively, against their analytical solutions. The solution of the Helmholtz PDE was tested on a dense regular grid of $256\times 256\times 100\times 100$ points $(x,y,a_{1},a_{2})\in\Omega\times\Pi$ against its analytical solution.

#### D.1.3 Leveraging invariance when solving BVPs (Sec. 5.3 )

The SCL problems formulated in this section are of the form (SCL)(M+I). To showcase the advantages of integrating additional knowledge, such as the structure of the BVP solution, we consider fixed collocation points for the constraints (SCL)(M). This is in fact not uncommon for PINNs, see, e.g., (Raissi et al., ; Lu et al., ). 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 periodic with period $2\pi/\beta$. We therefore use the problem

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle\frac{1}{N}\sum_{n=1}^{N}\Big{(}u_{\theta}(x_{n},t_{n})-h(x_{n},t _{n})\Big{)}^{2}$ $\displaystyle\operatorname{subject\ to}$ $\displaystyle\frac{1}{M}\sum_{m=1}^{M}\Big{(}D[u_{\theta}](x_{m},t_{m})-\tau(x _{m},t_{m})\Big{)}^{2}\leq\epsilon_{\text{pde}}$ $\displaystyle\operatorname{\mathbb{E}}_{(x,t)\sim\psi_{0}^{\text{ST}}}\!\bigg{ [}\Big{(}u_{\theta}(x,t)-u_{\theta}\Big{[}x,t+\frac{2\pi}{\beta}\Big{]}\Big{)} ^{2}\bigg{]}\leq\epsilon_{\text{s}}$ $$

For both (PI) and (SCL^′)(M) we use a total of $N=456$ collocation points, namely $256$ points $(x,0)$, $x\in[0,2\pi]$, for the IC and $100$ points equally spaced in $t\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 $\epsilon_{\text{pde}}=10^{-3}$ and $\epsilon_{s}=10^{-3}$. For (PI), we use the weights $\mu$ from (Daw et al., ), namely, $\mu_{\text{D}}=1$, $\mu_{\text{BC}}=100$, and $\mu_{\text{IC}}=100$.

##### Model.

We used MLPs with 4 hidden layers of 50 neurons.

##### Training.

For (SCL^′)(I), we used Alg. 2. For the invariance loss, we used $\Sigma=\operatorname{diag}(0.5,0.1)$ to sample from $(x,t)$. All models were trained for $200\,000$ epochs using Adam with the default parameters from (Kingma & Ba, ) and learning rate of $10^{-3}$ for (PI). For SCL, the dual learning rate was $\eta_{d}=10^{-4}$. We decayed the learning rates by a factor of $0.9$ every 5000 epochs.

##### Testing.

The solution was tested on a dense regular grid of $256\times 100$ points $(x,t)\in\mathcal{D}$ against its analytical solution.

#### D.1.4 Supervised solution of BVPs (Sec. 5.4 )

For supervised experiments, we formulate an SCL without objective using only data constraints (observational knowledge). Since we use FNOs, that can only make predictions on uniform grids, we replace $\psi_{0}^{OB}$ in (SCL) with a uniform distribution over a fixed regular grid. The problem the FNOs tackle is that of predicting the solution $u^{\dagger}$ of a BVP given its IC $h(x,0)$. Hence, the training data is composed of pairs $(u_{j}^{\dagger},h_{j})$ describing ICs and their corresponding solution. We therefore pose the SCL problem

$$ $\displaystyle\operatorname*{minimize}_{\theta\in\Theta}$ $\displaystyle 0$ $\displaystyle\operatorname{subject\ to}$ $\displaystyle\frac{1}{N}\sum_{n=1}^{N}\left(u_{\theta}(h_{j})(x_{n},t_{n})-u_{ j}^{\dagger}(x_{n},t_{n})\right)^{2}\leq\epsilon_{\text{o}},\quad j=1,\dots,J.$ $$

##### Problem hyperparameters.

For SCL, the tolerance was chosen as before, using a coarse hyperparameter search. The final values are reported in Table 6.

**Table 6: Problem hyperparameters for supervised solutions**
|  | $\epsilon_{\text{o}}$ | # training samples | # test samples | FNO architecture |
| --- | --- | --- | --- | --- |
| Burgers’ | $10^{-3}$ | 1000 | 200 | 16 modes, 4 layers |
| Diffusion-sorption | $10^{-3}$ | 1000 | 500 | 8 modes, 5 layers |
| Navier-Stokes: $\nu=10^{-3}$ | $10^{-2}$ | 1000 | 200 | 8 modes, 8 layers |
| Navier-Stokes: $\nu=10^{-4}$ | $5\times 10^{-2}$ | 1000 | 500 | 8 modes, 8 layers |
| Navier-Stokes: $\nu=10^{-5}$ | $10^{-2}$ | 1000 | 200 | 8 modes, 8 layers |

##### Model.

We used the FNO architecture from (Li et al., ) with 64 hidden channels, 128 projection channels, and no lifting channels. The number of modes and layers are reported in Table 6.

##### Training.

All models were trained for $500$ epochs using Adam with the default settings from (Kingma & Ba, ) with learning rate $10^{-3}$ and batch size of $20$. For SCL, the dual learning rate was $\eta_{d}=10^{-4}$. All learning rates were decreased by a factor of $0.5$ every $100$ epochs. All test errors are reported for the best model across epochs.

##### Testing.

The test sets from (Li et al., ) were used for Burgers’ and Navier-Stokes, whereas the diffusion-sorption test data was taken from (Takamoto et al., ). The sizes of the test sets are reported in Table 6.

## Appendix E Additional experiments

### E.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,

$$ $\text{Relative complexity}=\frac{\text{\# differential operator evaluations per epoch for \eqref{eq:scl_parametric}}}{\text{\# differential operator evaluations per epoch for \eqref{eq:PINN}}}\times 100\%$ $$

Recall that in order to evaluate the PDE loss, (PI) uses 1000 collocation points per discretized coefficient $\pi_{j}$ whereas (SCL^′) takes $5000$ steps of Alg. 2.

##### Convection equation.

Table 7 considers simultaneously solving all BVPs corresponding to the convection equation with $\beta\in[1,30]$ and compares (SCL^′)(M) with (PI). We see that (SCL^′)(M) outperforms or matches (PI) while being more efficient. In particular, (SCL^′)(M) significantly outperforms (PI) in terms of relative $L_{2}$ 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 a better compromise between error and computational cost. This is even clearer from Fig. 4, particularly when we normalize the $x$-axis in terms of differential operator evaluations.

**Table 7: Relative $L_{2}$ error and computational efficiency for the parametric convection problem.**
| Discretization for (PI) | Average Relative $L_{2}$ Error | Relative complexity<br>(SCL^′) $\div$ (PI) |  |
| --- | --- | --- | --- |
| (PI) | (SCL^′)(M) |  |  |
| $[1.0,10.0,20.0,30.0]$ | 0.365 | 0.0110 | 125% |
| $[1.0,5.0,10.0,15.0,20.0,25.0,30.0]$ | 0.220 | 71% |  |
| $[1.0,2.0,3.0,4.0,5.0,\ldots,30]$ | 0.0476 | 16% |  |

Figure: Figure 4: Relative $L_{2}$ error as a function of training epoch and differential operator evaluations for the parametric convection problem.
Refer to caption: x5.png

##### Reaction-diffusion equation.

Additional results for the parametric reaction-diffusion BVP are shown in Table 8. Same as for the convection PDE, (SCL^′)(M) is more efficient than (PI), striking a better compromise between computational complexity and performance. Indeed, in order to achieve the same error as (SCL^′)(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. 5), especially when the $x$-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 $\psi_{0}$ throughout training to understand where the advantage of SCL comes from (Fig. 6). First, we do not note any interesting behavior over the $x$-marginal (the samples are mostly uniform and the histogram is therefore omitted). Once again, we see that (SCL^′)(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. 1. As for the $(\nu,\rho)$, we see that the distributions shift during training, indicating the change in difficulty of fitting the solution of the reaction-diffusion PDE. In the end, the samples for $\rho$ are quite uniform, while we notice that there remains a strong focus on smaller values of $\nu$. Note that these distributions reflect the error patterns of the final solution (Fig. 7).

**Table 8: Relative $L_{2}$ error and computational efficiency for the parametric reaction-diffusion problem.**
| Coefficients<br>range | Discretization for (PI) | Average relative $L_{2}$ error | Relative complexity<br>(SCL^′) $\div$ (PI) |  |
| --- | --- | --- | --- | --- |
| (PI) | (SCL^′)(M) |  |  |  |
| $\nu\in[0,5]$<br>$\rho\in[0,5]$ | $\{0.0,2.5,5.0\}^{2}$ | 0.0793 | 0.0126 | 55.6% |
| $\{0.0,1.67,3.33,5.0\}^{2}$ | 0.0190 | 31.3% |  |  |
| $\{0.0,1.25,2.5,3.75,5.0\}^{2}$ | 0.0119 | 20% |  |  |
| $\{0.0,1.0,2.0,3.0,4.0,5.0\}^{2}$ | 0.0105 | 13.9% |  |  |
| $\nu\in[0,10]$<br>$\rho\in[0,10]$ | $\{0.0,5.0,10.0\}^{2}$ | 0.636 | 0.0133 | 55.6% |
| $\{0.0,2.5,5.0,7.5,10.0\}^{2}$ | 0.0228 | 20% |  |  |
| $\{0.0,2.0,4.0,6.0,8.0,10.0\}^{2}$ | 0.0131 | 13.9% |  |  |
| $\nu\in[1,20]$<br>$\rho\in[1,20]$ | $\{1.0,10.0,20.0\}^{2}$ | 0.841 | 0.0204 | 55.6% |
| $\{1.0,5.0,10.0,15.0,20.0\}^{2}$ | 0.0128 | 20% |  |  |

Figure: Figure 5: Relative $L_{2}$ error as a function of training epoch and differential operator evaluations for the parametric reaction-diffusion problem. (PI) uses the discretization $(\nu,\rho)\in\{0.0,2.0,4.0,6.0,8.0,10.0\}^{2}$
Refer to caption: x7.png

Figure: Figure 6: Histogram of (marginalized) MH samples of $\psi_{0}$ for the parametric reaction-diffusion equation.
Refer to caption: x9.png

Figure: Figure 7: Relative $L_{2}$ error for reaction-diffusion solutions trained using (SCL^′) and (PI) with discretization (a) $(\nu,\rho)\in\{0.0,2.5,5.0\}^{2}$ and (b) $(\nu,\rho)\in\{0.0,2.0,4.0,6.0,8.0,10.0\}^{2}$.
Refer to caption: x10.png

##### Helmholtz equation.

Once again, we see from Table 8 that (SCL^′)(M) makes more efficient use of computations than (PI). 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 $x$-axis in terms of computational complexity (Fig. 5). Naturally, taking finer discretizations eventually leads to lower errors (Fig. 9), but the computational cost associated also rises considerably. On the other hand, we 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 $\psi_{0}$ throughout training to gain a better understanding of the difficulties perceived by the MLP to fit solutions of this problem (Fig. 10). We display only the $x$ and $a_{1}$ marginals, seen as they display the same behaviors as $y$ and $a_{2}$ respectively due to the symmetry of the Helmholtz equation. Here, we notice that the distribution of $x$ has an alternating pattern initially. This makes sense seen as the solution of the Helmholtz equation is periodic. SCL clearly picks up on this pattern, focusing on 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 MLP clearly has difficulties fitting the solution. With respect to the problem coefficients, we notice that $\psi_{0}$ concentrates on larger values of $a_{1}$, 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. 9).

**Table 9: Relative $L_{2}$ error and computational efficiency for the parametric reaction-diffusion problem.**
| Coefficients<br>range | Discretization for (PI) | Average relative $L_{2}$ error | Relative complexity<br>(SCL^′) $\div$ (PI) |  |
| --- | --- | --- | --- | --- |
| (PI) | (SCL^′)(M) |  |  |  |
| $a_{1}\in[1,2]$<br>$a_{2}\in[1,2]$ | $\{1.0,1.5,2.0\}^{2}$ | 0.0307 | 0.0125 | 55.6% |
| $\{1.0,1.25,1.5,1.75,2.0\}^{2}$ | 0.00593 | 20% |  |  |
| $\{1.0,1.2,1.4,1.6,1.8,2.0\}^{2}$ | 0.00463 | 13.9% |  |  |
| $a_{1}\in[1,3]$<br>$a_{2}\in[1,3]$ | $\{1.0,2.0,3.0\}^{2}$ | 1.34 | 0.0549 | 55.6% |
| $\{1.0,1.5,2.0,2.5,3.0\}^{2}$ | 0.00943 | 20% |  |  |
| $\{1.0,1.4,1.8,2.2,2.6,3.0\}^{2}$ | 0.00953 | 13.9% |  |  |

Figure: Figure 8: Relative $L_{2}$ error as a function of training epoch and differential operator evaluations for the parametric Helmholtz problem. (PI) uses the discretization $(a_{1},a_{2})\in\{1.0,1.2,1.4,1.6,1.8,2.0\}^{2}$
Refer to caption: x12.png

Figure: Figure 9: Relative $L_{2}$ error for Helmholtz solutions trained using (SCL^′) and (PI) with discretization (a) $(a_{1},a_{2})\in\{1.0,1.5,2.0\}^{2}$ and (b) $(a_{1},a_{2})\in\{1.0,1.2,1.4,1.6,1.8,2.0\}^{2}$.
Refer to caption: x14.png

Figure: Figure 10: Histogram of (marginalized) MH samples of $\psi_{0}$ for the parametric Helmholtz equation.
Refer to caption: x16.png

### E.2 Supervised solution of BVPs

##### Burgers’ equation.

Fig 11 shows the box plots for the relative $L_{2}$ error across samples. They show that not only is the average error across samples smaller when using (SCL)(O), but in fact the whole 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 12). 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 in Fig. 3.

Figure: Figure 11: Distribution of train and test errors (across data points) for the Burgers’ equation (orange line indicates the median).
Refer to caption: x18.png

Figure: Figure 12: Dual variables obtained by Alg. 1 for the Burgers’ equation: (a) evolution during training and (b) distribution of the dual variables after training.
Refer to caption: x19.png

##### Diffusion-sorption equation.

We once again display the distribution of the relative $L^{2}$ errors across training and test data points for models trained using (PII) and (SCL)(O) (Fig. 13). Here, we clearly see that by bounding the maximum error rather than minimizing its average leads to a more homogeneous fit across samples. This is once again due to the different weight assigned to each data sample, weights that also evolve throughout training (Fig 14a). By inspecting the ICs with large and small values of $\lambda$, we notice the pattern showcased in Fig 14b, where IC with either large or small magnitude are more challenging to fit than those with moderate ones.

Figure: Figure 13: Distribution of train and test errors (across data points) for the diffusion-sorption equation (orange line indicates the median).
Refer to caption: x21.png

Figure: Figure 14: Dual variables obtained by Alg. 1 for the diffusion-sorption equation: (a) evolution during training and (b) value as a function of the IC magnitude.
Refer to caption: x22.png

##### Navier-Stokes equation.

Fig 15 shows that (SCL^′)(O) not only improves the average relative $L_{2}$ error, but its entire distribution across train and test data points. For the Navier-Stokes equations, however, it is harder to find a relation between IC properties and the difficulty of fitting the solution. That is because, as we show in Fig 16, the dual variables do not have such extremely different values. This is certainly due to the fact that the tolerance $\epsilon_{\text{o}}$ is set very loose ($0.01$) and that in these situations, all ICs are similarly difficult to fit. Still, some ICs have outlier values of $\lambda$ (Fig 16b), which does point to the fact that the FNO does struggle more to fit certain conditions. That being said, it not easy to identify what in those conditions make them hard (Fig 17). Nevertheless, this is not an issue as we need not know beforehand which ICs are challenging: suffices it to run Alg. 1 to solve (SCL)(O).

Figure: Figure 15: Distribution of train and test errors (across data points) for the Navier-Stockes equation with $\nu=10^{-3}$ (orange line indicates the median).
Refer to caption: x24.png

Figure: Figure 16: Dual variables obtained by Alg. 1 for the Navier-Stokes equation with $\nu=10^{-3}$: (a) evolution during training and (b) distribution of the dual variables after training.
Refer to caption: x25.png

Figure: Figure 17: Initial conditions corresponding to the smallest and largest final dual variables for Navier-Stokes equation with $\nu=10^{-3}$.
Refer to caption: x27.png

## Appendix F Additional related work

### F.1 Physics-informed neural networks

Fitting MLPs to the solution of BVPs goes back to (Psichogios & Ungar, ; Dissanayake & Phan-Thien, ; Lagaris et al., ). 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., ; Wight & Zhao, ; Chen et al., ; Lu et al., ; Basir & Senocak, ; Yu et al., ; Xu et al., )). This led to new architectures tailored for PDEs (Raissi et al., ; Fathony et al., ; Gao et al., ; Wang et al., ; Kang et al., ; Moseley et al., ; Cho et al., ; Chalapathi et al., ), leveraging positional embedding (Wang et al., ) and adaptive activation functions (Jagtap et al., ).

##### 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., ; Markidis, ; Wight & Zhao, ; Wang et al., ; ; ; Grossmann et al., ). Based on these observations, adaptive heuristics have been proposed to select the collocation points based on importance sampling (Nabian et al., ; Wu et al., ), adversarial training Wang et al. (), rejection sampling (Daw et al., ), and causality-inspired rules (Wang et al., ). Similarly, empirical rules for determining the loss weights [$\mu$ in (PI)] have been developed using the magnitude of the gradients (Wang et al., ), eigenvalues of the neural tangent kernel (Wang et al., ), inverse-Dirichlet weighting (Maddu et al., ), soft attention mechanisms (McClenny & Braga-Neto, ), and (augmented) Lagrangian formulations (Lu et al., ; Basir & Senocak, ). Other works have addressed these challenges by changing the problem formulation inspired by traditional numerical methods (Kharazmi et al., ; Chiu et al., ; Patel et al., ), proposing different objective functions (Yu et al., ; Son et al., ), and using sequential (Wight & Zhao, ; Krishnapriyan et al., ) and transfer learning (Goswami et al., ; Chakraborty, ; Desai et al., ) techniques.

In contrast to 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 use either adversarial training to estimate the worst-case loss as in (Wang et al., ) or constrained learning to manipulate the loss weights as in (Lu et al., ; Basir & Senocak, ). Both are required simultaneously. Leveraging findings from adversarially robust learning, we also replace the gradient methods used in (Wang et al., ), i.e., the technique from (Mądry et al., ), by the sampling-based approach in (Robey* et al., ).

### F.2 Neural Operators

In contrast to the MLPs and convolutional NNs (CNNs) typically used in PINNs, NOs are NNs capable of handling infinite-dimensional inputs and outputs. They can therefore be trained to find BVP solutions for different IC or forcing functions. Many different architectures have been proposed, such as DeepONets (Lu et al., ), FNOs (Li et al., ), and NO based on U-Nets (Gupta & Brandstetter, ). FNOs in particular have become quite popular and garnered many efforts towards addressing their limitations, such as improving memory efficiency (Rahman et al., ), designing equivariant FNOs (Helwig et al., ), extending FNOs to general geometries (Li et al., ), factorizing the Fourier transform (Tran et al., ), and leveraging multiwavelets (Gupta et al., ).

##### Training NOs.

Regardless of these improvements, the vast majority of NOs are trained in a supervised manner by minimizing their average error across samples as in (PII) (see, e.g., (Lu et al., ; Li et al., ; Kovachki et al., )). 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-supervised techniques involving PDE losses have also been used (Li et al., ), computing the space-time derivatives needed to evaluate $D_{\pi}[u]$ is challenging for NOs.

In this paper, we do not develop new NO architectures, but focus on the problem of training them. Explicitly, rather than targeting the average error, we target the maximum error across samples. This is 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.

### F.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 a technique to train ML systems under requirements, such as fairness (Kearns et al., ; Cotter et al., ; Chamon & Ribeiro, ; Chamon et al., ) and robustness (Chamon & Ribeiro, ; Robey* et al., ; Hounie et al., ; Chamon et al., ), or to handle applications in which we want to attain good performance with respect to more than one metric. As in unconstrained learning, it is formulated as 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 to a practical learning rule with generalization guarantees (Chamon & Ribeiro, ; Chamon et al., ). These duality results have also been exploited to automatically adapt each constraint specification to their underlying difficulty, striking better compromises between objective and requirements Hounie et al. ().

Aside from dealing with the nominal accuracy vs. robustness trade-off typical in ML systems, constrained learning has itself been used to optimize robust losses as in (4). Typically, this is done using some combination of gradient ascent and random initialization, restart, and pruning heuristics (Goodfellow et al., ; Mądry et al., ; Dhillon et al., ; Wu et al., ; Cheng et al., ). Using semi-infinite optimization techniques, however, these deterministic methods can be replaced by a sampling approach that has been successful in a variety of domains (Robey* et al., ). Different MCMC methods (Robert & Casella, ) have been used in this context, including Langevin Monte Carlo (LMC) (Robey* et al., ) and Metropolis-Hastings (MH) (Hounie et al., ).

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 4
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).
soundness: 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.
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.

The 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.
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:

McGreivy, 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.

The 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.

And 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.

So we might instead take the contribution of the paper to be the guarantee, except:

1) As stated in the paper, existing methods already have guarantees.

2) 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.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
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.
soundness: 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.
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.
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. 

One 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. 

The 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.

[1] Kentaro Hoshisashi, Carolyn E Phelan, and Paolo Barucca. No-arbitrage deep calibration for volatility smile and skewness. arXiv preprint arXiv:2310.16703, 2023.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: The specific questions (associated with the items above) are:

+ 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)?
soundness: 3
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
with worst-case losses.

The 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.
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
with worst-case losses.
weaknesses: The reviewer is having a difficult time distinguishing the contribution from the work of Paris Perdikaris and his group.  Examples are:

A. Daw, J. Bu, S. Wang, P. Perdikaris, A. Karpatne, Rethinking the Importance of Sampling in Physics-informed Neural
Networks, arXiv preprint arXiv:2207.02338

 S. Wang, S. Sankaran, P. Perdikaris, Respecting causality is all you need for training physics-informed neural networks,
arXiv preprint arXiv:2203.07404

Paris 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).   

The 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.

At 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:

https://www.sciencedirect.com/science/article/pii/S1385894724012919#b78

Is 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.

In terms of concrete actions, can the authors:
+  Provide a more comprehensive literature review that clearly positions their work relative to existing constrained learning approaches for PINNs?

+  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)?

+ Highlight any novel theoretical results or empirical improvements over state-of-the-art methods

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
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? 

- (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?

- (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.

- (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?

- Under what conditions, does Algorithm 1 converge?  (Chamon et al., 2023; Elenter et al., 2024) could be valuable references, to my understanding.
soundness: 2
strengths: - The authors tackle important problems in physics-informed learning: reducing computational costs and avoiding hyperparameter tuning.

- The experimental results show that the proposed formulation outperforms baselines in some cases.
summary: This paper developed SCL, a technique for solving BVPs based on constrained learning. 
It then developed a practical algorithm to tackle these problems and showcased its performance across a variety of PDEs.
SCL not only yields accurate BVP solutions, but tackles many of the challenges faced by previous methods, such as extensive hyperparameter tuning and computational costs.
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.

- The dependence on random seeds is unknown, potentially causing reproducibility issues. Error bars should be added.

- 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.

- I have several critical questions. Please see Questions below. I am open to discussion.

- (Minor) (Lines 1730, 1733, etc.) Fig -> Fig.

Meta-review: This work shows that obtaining weak solutions to PDEs can be formulated as a constrained learning problem with worst case losses. Furthermore, it develops an efficient algorithm for solving this problem and incorporating further structural information that maybe be known. The method is evaluated on a variety of tasks and the comparisons are extensive.

Decision: Accept (Poster)