# Inverse-Dirichlet Weighting Enables Reliable Training of Physics Informed
Neural Networks

Suryanarayana Maddu1,2,3,7, Dominik Sturm9,10, Christian L. Müller4,5,6, Ivo
F. Sbalzarini1,2,3,7,8<EMAIL_ADDRESS>1 Technische Universität Dresden,
Faculty of Computer Science, 01069 Dresden, Germany  2 Max Planck Institute of
Molecular Cell Biology and Genetics, 01307 Dresden, Germany  3 Center for
Systems Biology Dresden, 01307 Dresden, Germany  4 Center for Computational
Mathematics, Flatiron Institute, New York, NY, USA  5 Department of
Statistics, LMU München, Munich, Germany  6 Institute of Computational
Biology, Helmholtz Zentrum München, Germany  7 Center for Scalable Data
Analytics and Artificial Intelligence ScaDS.AI, Dresden/Leipzig, Germany  8
Cluster of Excellence Physics of Life, TU Dresden, Germany  9 Helmholtz-
Zentrum Dresden-Rossendorf, D-01328 Dresden, Germany  10 Center for Advanced
Systems Understanding (CASUS), D-02826 Görlitz, Germany

###### Abstract

We characterize and remedy a failure mode that may arise from multi-scale
dynamics with scale imbalances during training of deep neural networks, such
as Physics Informed Neural Networks (PINNs). PINNs are popular machine-
learning templates that allow for seamless integration of physical equation
models with data. Their training amounts to solving an optimization problem
over a weighted sum of data-fidelity and equation-fidelity objectives.
Conflicts between objectives can arise from scale imbalances,
heteroscedasticity in the data, stiffness of the physical equation, or from
catastrophic interference during sequential training. We explain the training
pathology arising from this and propose a simple yet effective inverse-
Dirichlet weighting strategy to alleviate the issue. We compare with Sobolev
training of neural networks, providing the baseline of analytically
$\bm{\epsilon}$-optimal training. We demonstrate the effectiveness of inverse-
Dirichlet weighting in various applications, including a multi-scale model of
active turbulence, where we show orders of magnitude improvement in accuracy
and convergence over conventional PINN training. For inverse modeling using
sequential training, we find that inverse-Dirichlet weighting protects a PINN
against catastrophic forgetting.

††preprint: APS/123-QED

##  I Introduction

Data-driven modeling has emerged as a powerful and complementary approach to
first-principles modeling. It was made possible by advances in imaging and
measurement technology, computing power, and innovations in machine learning,
in particular neural networks. The success of neural networks in data-driven
modeling can be attributed to their powerful function approximation properties
for diverse functional forms, from image recognition [1] and natural language
processing [2] to learning policies in complex environments [3].

Recently, there has been a surge in exploiting the versatile approximation
properties of neural networks for surrogate modeling in scientific computing
and for data-driven modeling of physical systems [4, 5]. These applications
hinge on the approximation properties of neural networks for Sobolev-regular
functions [6, 7]. The popular Physics Informed Neural Networks (PINNs) rely on
knowing a differential equation model of the system in order to solve a soft-
constrained optimization problem [8, 5]. Thanks to their mesh-free character,
PINNs can be used for both forward and inverse modeling in domains including
material science [9, 10, 11], fluid dynamics [8, 12, 13, 14] and turbulence
[15, 16], biology [17], medicine [18, 19], earth science [20], mechanics [21],
uncertainty quantification [22], as well as stochastic [23], high-dimensional
[24], and fractional differential equations [25].

The wide application scope of PINNs rests on their parameterization as deep
neural networks, which can be tuned by minimizing a weighted sum of data-
fitting and equation-fitting errors. However, finding network parameters that
impartially optimize for several errors is challenging when errors impose
conflicting requirements. For Multi-Task Learning (MTL) problems in machine
learning, this issue has been addressed by various strategies for relative
weighting of conflicting objectives, e.g., based on uncertainty [26], gradient
normalization [27], or Pareto optimality [28]. The PINN community in parallel
has developed heuristics for loss weighting with demonstrated gains in
accuracy and convergence [29]. However, there is neither consensus on when to
use such strategies in PINNs, nor are there optimal benchmark solutions for
objectives based on differential equations.

Here, we characterize training pathologies of PINNs and provide criteria for
their occurrence. We explain these pathologies by showing a connection between
PINNs and Sobolev training, and we propose a strategy for loss weighting in
PINNs based on the Dirichlet energy of the task-specific gradients. We show
that this strategy reduces optimization bias and protects against catastrophic
forgetting. We evaluate the proposed inverse-Dirichlet weighting by comparing
with Sobolev training in a case where provably optimal weights can be derived,
and by empirical comparison with two conceptually different state-of-the-art
PINN weighting approaches.

##  II A Physics Informed Neural Network is a Multi-Objective Optimization
Problem

An optimization problem involving multiple tasks or objectives can be written
as:

$\textrm{minimize }\mathcal{L}_{k}(\bm{\theta}^{sh},\bm{\theta}^{k}),\qquad
k=1,\ldots,K.$ (1)

The total number of tasks is given by $K$. $\bm{\theta}^{sh}$ are shared
parameters between the tasks, and $\bm{\theta}^{k}$ are the task-specific
parameters. For conflicting objectives, there is no single optimal solution,
and Pareto dominated solutions can be found using heuristics like evolutionary
algorithms [30]. Alternatively, the multi-objective optimization problem can
be converted to a single-objective optimization problem using scalarization,
e.g., as a weighted sum with real-valued relative weights $\lambda_{k}$
associated with each loss objective $\mathcal{L}_{k}$:

$\min_{\bm{\theta}^{sh},\bm{\theta}^{1},\ldots,\bm{\theta}^{K}}\sum_{k=1}^{K}\lambda_{k}\mathcal{L}_{k}(\bm{\theta}^{sh},\bm{\theta}^{k}).$
(2)

Often, the $\lambda_{k}$ are dynamically tuned along training epochs [28]. The
weights $\lambda_{k}$ intuitively quantify the extent to which we aim to
minimize the objective $\mathcal{L}_{k}$. Large $\lambda_{k}$ favor small
$\mathcal{L}_{k}$, and vice-versa [31]. For $K\geq 2$, manual tuning and grid
search become prohibitively expensive for large parameterizations as common in
deep neural networks.

The loss function of a PINN is a scalarized multi-objective loss where each
objective corresponds to either a data-fidelity term or a constraint derived
from the physics of the underlying problem [5, 4]. A typical PINN loss
contains five objectives:

$\mathcal{L}(\bm{u}_{\bm{\theta}},\bm{u})=\frac{\lambda_{\text{1}}}{N_{\text{int}}}\sum_{i,j=1}^{N_{\text{int}}}\Big{|}\partial_{t}\bm{u}_{\bm{\theta}}-\mathcal{N}(\bm{u}_{\bm{\theta}}(t_{j},\bm{x}_{i}^{\Omega}),\Sigma)\Big{|}^{2}+\\\
\frac{\lambda_{2}}{N_{\mathcal{B}}}\sum_{i,j=1}^{N_{\mathcal{B}}}\Big{|}\mathcal{B}(t_{j},\bm{x}_{i}^{\partial\Omega})\Big{|}^{2}+\frac{\lambda_{3}}{N_{\text{I}}}\sum_{i=1}^{N_{\text{I}}}\Big{|}\mathcal{I}(t=0,\bm{x}_{i}^{\Omega})\Big{|}^{2}+\\\
\frac{\lambda_{4}}{N_{\text{int}}}\sum_{i,j=1}^{N_{\text{int}}}\Big{|}\mathcal{F}(\bm{u}_{\bm{\theta}}(t_{j},\bm{x}_{i}^{\Omega}))\Big{|}^{2}+\frac{\lambda_{5}}{N_{\text{int}}}\sum_{i,j=1}^{N_{\text{int}}}\Big{|}\bm{u}_{\bm{\theta}}(t_{j},\bm{x}_{i}^{\Omega})-\bm{u}_{i})\Big{|}^{2}.$
(3)

We distinguish the neural network approximation
$\bm{u}_{\bm{\theta}}(t_{j},\bm{x}_{i}^{\Omega})$ from the training data
$\bm{u}_{i}=\bm{u}(t_{j},\bm{x}_{i}^{\Omega})$ sampled at the $N_{\text{int}}$
space-time points $(t_{j},\bm{x}_{i}^{\Omega})$ in the interior of the
solution domain $\Omega$, $N_{\mathcal{B}}$ points
$(t_{j},\bm{x}_{i}^{\partial\Omega})$ on the boundary $\partial\Omega$, and
$N_{\text{I}}$ samples of the initial state at time $t=0$. The $\lambda_{1}$,
$\lambda_{2}$, and $\lambda_{3}$ weight the residuals of the physical equation
$\partial_{t}\bm{u}=\mathcal{N}(\bm{u},\Sigma)$ with right-hand side
$\mathcal{N}$ and coefficients $\Sigma$, the boundary condition $\mathcal{B}$,
and the initial condition $\mathcal{I}$, respectively. The objective with
weight $\lambda_{4}$ is an additional constraint $\mathcal{F}$ on the state
variable $\bm{u}$, such as divergence-freeness ($\nabla\cdot\bm{u}=0$) of the
velocity field of an incompressible flow. The weight $\lambda_{5}$ penalizes
deviations of the neural network approximation from the training data. From
now on, we simplify the notation by designating space-time points as
$\bm{x}_{i}:=(t,\bm{x})_{i}=(t_{j},\bm{x}_{k})$ in time-dependent problems.

PINNs can be used to solve both forward and inverse modeling problems of
ordinary differential equations (ODEs) or partial differential equations
(PDEs). In the forward problem, a numerical approximation to the solution of a
ODE or PDE is to be learned. For this, $\bm{\theta}^{\text{sh}}$ are the
networks parameters, and $\bm{\theta}^{k}$ is an empty set. In the inverse
problem, the coefficients of an ODE or PDE are to be learned from data such
that the equation describes the dynamics in the data. Then, the task-specific
parameters $\bm{\theta}^{k}$ correspond to the coefficients
($\bm{\xi}\subseteq\Sigma$) of the physical equation that shall be inferred.
During training of a PINN, the parameters
$(\bm{\theta}^{\text{sh}},\bm{\theta}^{k})$ are determined by minimizing the
total loss $\mathcal{L}$ for given data $\bm{u}(\bm{x}_{i})$.

###  II.1 Sobolev training is a special case of PINN training

Like most deep neural networks, PINNs are usually trained by first-order
parameter update based on (approximate) loss gradients. The update rule can be
interpreted as a discretization of the gradient flow, i.e.,

$\displaystyle\bm{\theta}(\tau+1)$
$\displaystyle=\bm{\theta}(\tau)-\eta(\tau)\sum_{k=1}^{K}\lambda_{k}\nabla_{\bm{\theta}}\mathcal{L}_{k}.$
(4)

The learning rate $\eta(\tau)$ depends on the training epoch $\tau$,
reflecting the adaptive step size used in optimizers like Adam [32]. The
gradients are batch gradients
$\nabla_{\bm{\theta}}\mathcal{L}_{k}=\frac{1}{|\textrm{B}|}\sum_{i\in\textrm{B}}\nabla\ell_{k}(\bm{u}_{\bm{\theta}}(\bm{x}_{i}),\bm{u}_{i})$,
where B is a mini-batch created by random splitting of the $N$ training data
points into $N/|\textrm{B}|$ distinct partitions and $\ell_{k}$ is one summand
(for one data point $\bm{x}_{i}$) of the loss $\mathcal{L}_{k}$ corresponding
to the $k^{\text{th}}$ objective. For batch size $|\textrm{B}|=1$, the mini-
batch gradient descent reduces to stochastic gradient descent.

Adding derivatives $\mathcal{D}_{\textbf{x}}^{m}$, $m\geq 1$, of the state
variable into the training loss has been shown to improve data efficiency and
generalization power of neural networks [7]. Minimizing the resulting loss

$\sum_{i=1}^{N}\Bigg{[}\lambda_{0}|\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\bm{u}(\bm{x}_{i})|^{2}+\sum_{k=1}^{K}\lambda_{k}|\mathcal{D}_{\textbf{x}}^{k}\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\mathcal{D}_{\textbf{x}}^{k}\bm{u}(\bm{x}_{i})|^{2}\Bigg{]}$
(5)

using the update rule in Eq. 4 can be interpreted as a finite approximation to
the Sobolev gradient flow of the loss functional. The neural network is then
trained with access to the ($K+2$)-tuples
$\\{(\bm{x}_{i},\bm{u}(\bm{x}_{i}),\mathcal{D}_{\textbf{x}}^{1},\ldots,\mathcal{D}_{\textbf{x}}^{K})\\}_{i=1}^{N}$
instead of the usual training set
$\\{(\bm{x}_{i},\bm{u}(\bm{x}_{i}))\\}_{i=1}^{N}$. Such Sobolev training is an
instance of scalarized multi-objective optimization with each objective based
on the approximation of a specific derivative. The weights are usually chosen
as $\lambda_{0}=\ldots=\lambda_{\text{K}}=1$ [7]. The additional information
from the derivatives introduces an inductive bias, which has been found to
improve data efficiency and network generalization [7]. This is akin to PINNs,
which use derivatives from the PDE model to achieve inductive bias. Indeed,
the Sobolev loss in Eq. 5 is a special case of the PINN loss from Eq. 3 with
each objective based on a PDE with right-hand side
$\mathcal{N}=\mathcal{D}_{\textbf{x}}^{k}\bm{u}$, $k=1,\ldots,K$.

###  II.2 Vanishing task-specific gradients are a failure mode of PINNs

The update rule in Eq. 4 is known to lead to stiff learning dynamics for loss
functionals whose Hessian matrix has large positive eigenvalues [29]. This is
the case for Sobolev training, even with $K=1$, when the function
$\bm{u}(\bm{x})$ is highly oscillatory, leading us to the following
proposition, proven in the Supplementary Material:

###### Proposition II.2.1

[Stiff learning dynamics in first-order parameter update] For Sobolev training
on the tuples
$\\{\bm{x}_{i},\bm{u}(\bm{x}_{i}),\partial_{\bm{x}}^{m}\bm{u}(\bm{x}_{i})\\}_{i=1}^{N},m\geq
1$, with total loss
$\mathcal{L}=\lambda_{1}\sum_{i=1}^{N}|\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\bm{u}(\bm{x}_{i})|^{2}+\lambda_{2}|\mathcal{D}_{\bm{x}}^{m}\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\mathcal{D}_{\bm{x}}^{m}\bm{u}(\bm{x}_{i})|^{2}$,
the rate of change of the residual of the dominant Fourier mode
$\widetilde{\bm{r}}(\bm{k}_{0})=\widetilde{\bm{u}}_{\bm{\theta}}(\bm{k}_{0})-\widetilde{\bm{u}}(\bm{k}_{0})$
using first-order parameter update is:

$\displaystyle\Big{|}\frac{\mathrm{d}\widetilde{\bm{r}}({\bm{k}_{0}})}{\mathrm{d}\tau}\Big{|}$
$\displaystyle=\frac{\eta}{O(\bm{k}_{0})}\Big{[}\lambda_{1}O(1)+\lambda_{2}O(\bm{k}_{0}^{2m})\Big{]}\Big{|}\frac{\partial\widetilde{\bm{r}}({\bm{k}_{0}})}{\partial\bm{\theta}}\Big{|},$

where $O(\cdot)$ is the Bachmann-Landau big-O symbol. For
$\lambda_{1}=\lambda_{2}=1$, this leads to training dynamics with slow time
scale $\mathrm{d}/\mathrm{d}t\in O(\bm{k}_{0}^{-1})$ and fast time scale
$\mathrm{d}/\mathrm{d}t\in O(\bm{k}_{0}^{2m-1})$ for
$\widetilde{\bm{r}}({\bm{k}_{0}})\in O(1)$ ( Proof in supplementary ).

This proposition suggests that losses containing high derivatives $m$,
dominant high frequencies $\bm{k}_{0}$, or a combination of both exhibit
dominant gradient statistics. This can lead to biased optimization of the
dominant objectives at the expense of the other objectives with relatively
smaller gradients. This is akin to the well known vanishing gradients
phenomenon across the layers of a neural network [33], albeit now across
different objectives of a multi-objective optimization problem. We call this
phenomenon _vanishing task-specific gradients_.

Proposition II.2.1 also admits a qualitative comparison to the phenomenon of
numerical stiffness in PDEs of the form

$\frac{\partial{\widetilde{\bm{u}}}}{\partial
t}=\mathbf{L}\widetilde{\bm{u}}+\mathbf{F}(\widetilde{\bm{u}},t)$ (6)

If the real parts of the eigenvalues of $\mathbf{L}$,
$c=\mathrm{Re}(\lambda(\mathbf{L}))\gg 1$, then the solution of the above PDE
has fast modes $\widetilde{\bm{u}}\in O(1)$ with $\mathrm{d}/\mathrm{d}t\in
O(c)$ and slow modes $\widetilde{\bm{u}}\in O(1/c)$ with
$\mathrm{d}/\mathrm{d}t\in O(1)$. In words, high frequency modes evolve on
shorter time scales than low frequency modes. The rate amplitude for spatial
derivatives of order $m$ is $\mathrm{d}/\mathrm{d}t\in O(\bm{k}^{-m})$ for the
wavenumber $\bm{k}$ [34]. Due to this analogy, we call the learning dynamics
of a neural network with vanishing task-specific gradients _stiff_. Stiff
learning dynamics leads to discrepancy in the convergence rates of different
objectives in a loss function.

Taken together, vanishing task-specific gradients constitute a likely failure
mode of PINNs, since their training amounts to an extended version of Sobolev
training.

##  III Strategies for Balanced PINN Training

Since PINN training minimizes a weighted sum of multiple objectives (Eq. 3),
the result depends on appropriately chosen weights. The ratio between any two
weights $\lambda_{i}/\lambda_{j}$, $i\neq j$, defines the relative importance
of the two objectives during training. Suitable weighting is therefore
quintessential to finding a Pareto-optimal solution that balances all
objectives. While ad hoc manual tuning is still commonplace, it has recently
been proposed to determine the weights as moving averages over the inverse
gradient magnitude statistics of a given objective as [29]:

$\displaystyle\hat{\lambda}_{k}(\tau)$
$\displaystyle=\frac{\max\\{|\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{\text{1}}(\tau)|\\}}{\lambda_{k}(\tau)\overline{|\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}(\tau)|}},$
(7) $\displaystyle\lambda_{k}(\tau+1)$
$\displaystyle=\alpha\lambda_{k}(\tau)+(1-\alpha)\hat{\lambda}_{k}(\tau),$

where $\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{\text{1}}$ is the gradient of the
residual, i.e. of the first objective in Eq. 3, and $\alpha$ is a user-defined
learning rate (here always $\alpha=0.5$, see Supplementary Material). All
gradients are taken with respect to the shared parameters across tasks and
$|\cdot|$ is the component-wise absolute value. The overbar signifies the
algebraic mean over the vector components. The maximum in the numerator is
over all components of the gradient vector. While this weighting has been
shown to improve PINN performance [29], it does not prevent vanishing task-
specific gradients.

###  III.1 Inverse-Dirichlet weighting

Instead of determining the relative weights proportional to the inverse
average gradient magnitude, we here propose to use weights based on the
gradient variance. In particular, we propose to use weights for which the
variances over the components of the back-propagated weighted gradients
$\lambda_{k}\nabla_{\bm{\theta}}\mathcal{L}_{k}$ become equal across all
objectives, thus directly preventing vanishing task-specific gradients. This
can be achieved in any (stochastic) gradient descent optimizer by using the
update rule in Eq. 4 with weights

$\displaystyle\hat{\lambda}_{k}(\tau)$ $\displaystyle=\frac{\max_{t=1\ldots
K}\left(\text{Var}[\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{t}(\tau)]\right)}{\text{Var}[\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}(\tau)]},$
(8) $\displaystyle\lambda_{k}(\tau+1)$
$\displaystyle=\alpha\lambda_{k}(\tau)+(1-\alpha)\hat{\lambda}_{k}(\tau).$

This guarantees that all weighted gradients have the same variance over
gradient components, i.e., that
$\text{Var}[\lambda_{k}\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}(\tau)]\approx\max_{t=1\ldots
K}\text{Var}[\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{t}(\tau)]\>\forall
k=1,\ldots,K$. Since the popular Adam optimizer is invariant to diagonal
scaling of the gradients [32], the weights in Eq. 8 can in this case
efficiently be computed as the inverse of averaged squared gradients of the
task-specific objectives, which in turn is proportional to the Dirichlet
energy of the objective, i.e.,
$\mathbb{E}[\left(\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}\right)^{2}]\propto\int_{\Omega_{\bm{\theta}}}|\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}|^{2}\mathrm{d}\bm{\theta}$.
We therefore refer to this weighting strategy as _inverse-Dirichlet
weighting_.

The weighting strategies in Eqs. 7 and 8 are fundamentally different. Eq. 7
weights objectives in inverse proportion to _certainty_ as quantified by the
average gradient magnitude. Eq. 8 uses weights that are inversely proportional
to _uncertainty_ as quantified by the variance of the loss gradients. Equation
8 is also different from previous uncertainty-weighted approaches [26] because
it measures the _training uncertainty_ stemming from noise in the loss
gradients rather than the uncertainty from the model’s _observational noise_.

###  III.2 Gradient-based multi-objective optimization

We compare weighting-based scalarization approaches with gradient-based multi-
objective approaches using Karush-Kuhn-Tucker (KKT) local optimality
conditions to find a descent direction that decreases all objectives.
Specifically, we adapt the multiple gradient descent algorithm (MGDA) [28],
which leverages the KKT conditions

$\exists\bm{\lambda}=\\{\lambda_{k}\\}\in\mathbb{R}_{+}^{K}\>\textrm{
s.t.}\sum_{k=1}^{K}\lambda_{k}\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}(\bm{\theta}^{sh},\bm{\theta}^{k})=0,\>\sum_{k=1}^{K}\lambda_{k}=1$
(9)

to solve an MTL problem. The MGDA is guaranteed to converge to a Pareto-
stationary point [35]. Given the gradients
$\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{k}$ of all objectives $k=1,\ldots,K$
with respect to the shared parameters, we find weights $\lambda_{k}$ that
satisfy the above KKT conditions. The resulting $\bm{\lambda}$ leads to a
descent direction that improves all objectives [28], see also Suppl. Fig. A.1.
We adapt the MGDA to PINNs as described in the Materials and Methods section.
We refer to this adapted algorithm as _pinn-MGDA_.

##  IV Results

As we posit above, PINNs are expected to fail for problems with dominant high
frequencies, i.e, with high $\bm{k}_{0}$ or high derivatives $m$, as is
typical for multi-scale PDE models. We therefore present results of numerical
experiments for such cases and compare the approximation properties of PINNs
with different weighting schemes. This showcases learning failure modes
originating from vanishing task-specific gradients.

We first consider Sobolev training of neural networks, which is a special case
of PINN training, as shown above. Due to the linear nature of the Sobolev loss
(Eq. 5), this test case allows us to compute $\epsilon$-optimal analytical
weights (see Materials and Methods), providing a baseline for our comparison.
Second, we consider a real-world example for a nonlinear PDE, using PINNs to
model active turbulence.

In the first example, we specifically consider the Sobolev training problem
with target function derivatives up to fourth order, i.e. $K=4$, and total
loss

$\sum_{i=1}^{N}\Bigg{[}\big{|}\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\bm{u}(\bm{x}_{i})\big{|}^{2}+\sum_{k=1}^{4}\lambda_{k}\big{|}\hat{\xi}_{k}\mathcal{D}_{\textbf{x}}^{k}\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\xi_{k}\mathcal{D}_{\textbf{x}}^{k}\bm{u}(\bm{x}_{i})\big{|}^{2}\Bigg{]}.$
(10)

The neural network is trained on the data $\bm{u}(\bm{x}_{i})$ and their
derivatives. To mimic inverse-modeling scenarios, in addition to producing an
accurate estimate $\bm{u}_{\bm{\theta}}(\bm{x}_{i})$ of the function, the
neural network is also tasked to infer unknown scalar prefactors
$\bm{\hat{\xi}}=(\hat{\xi}_{1},\hat{\xi}_{2},\hat{\xi}_{3},\hat{\xi}_{4})$
with true values chosen as $\xi_{k}=1$. This mimics scenarios in which unknown
coefficients of the PDE model also needed to be inferred from data. For this
loss, $\epsilon$-optimal weights can be analytically determined (see Materials
and Methods) such that all objectives are minimized in an unbiased fashion
[36]. This provides the baseline against which we benchmark the different
weighting schemes.

###  IV.1 Inverse-Dirichlet weighting avoids vanishing task-specific gradients

We characterize the phenomenon of vanishing task-specific gradients,
demonstrate its impact on the accuracy of a learned function approximation,
and empirically validate Proposition II.2.1. For this, we consider a generic
neural network with 5 layers and 64 neurons per layer, tasked with learning a
2D function that contains a wide spectrum of frequencies and unknown
coefficients
$\bm{\hat{\xi}}=(\hat{\xi}_{1},\hat{\xi}_{2},\hat{\xi}_{3},\hat{\xi}_{4})$
using Sobolev training with the loss from Eq. 10. The details of the test
problem and the training setup are given in the Supplementary Material.

Figure 1: Gradient histograms illustrating vanishing task-specific gradients.
The gradient distributions of the different objectives ($k=0,\ldots,4$, see
main text) are shown in different colors (see inset legend) at training epoch
5000. The panels correspond to different $\lambda$ weighting schemes as named.

We first confirm that in this test case vanishing task-specific gradients
occur when using uniform weights, i.e., when setting $\lambda_{k}=1$ for all
$k=0,\ldots,4$. For this, we plot the gradient histogram of the five loss
terms in Fig. 1A. The histograms clearly show that training is biased in this
case to mainly optimize the objective containing the highest derivative $k=4$.
All other objectives are neglected, as predicted by Proposition II.2.1. This
leads to a uniformly high approximation error on the test data (training vs.
test data split 50:50, see Supplementary Material), as shown in Fig. 2B, as
well as a uniformly high error in the estimated coefficients $\bm{\xi}$ (Fig.
2C). When determining the weights using the Pareto-seeking pinn-MGDA, the
accuracy considerably improves (Fig. 2B,C). However, the phenomenon of
vanishing task-specific gradients is still present, as seen in Fig. 1B, this
time favoring the objective for $k=0$ and neglecting the derivatives. Max/avg
weighting as proposed in Ref. [29] and given in Eq. 7 improves the accuracy
for higher derivatives (Fig. 2B), but still suffers from unbalanced gradient
distributions (Fig. 1C), which leads to a uniformly high error in the
estimated coefficients $\bm{\xi}$ (Fig. 2C). The inverse-Dirichlet weighting
proposed here and detailed in Eq. 8 shows balanced gradient histograms (Fig.
1D) closest to those observed when using $\epsilon$-optimal analytical weights
(Fig. 1E). This leads to orders of magnitude better accuracy in both the
function approximation and the coefficient estimates, as shown in Fig. 2B,C,
reaching the performance achieved when using $\epsilon$-optimal weights.

This is confirmed when looking at the power spectrum of the function
approximation $\bm{u}_{\bm{\theta}}$ learned by the neural network in Fig. 2D.
Inverse-Dirichlet weighting leads to results that are as good as those
achieved by optimal weights in this test case, in particular to approximate
high frequencies, as predicted by Proposition II.2.1. A closer look at the
power spectra learned by the different weighting strategies also reveals that
the max/avg strategy fails to capture high frequencies, whereas the pinn-MDGA
performs surprisingly well (Fig. 2D) despite its unbalanced gradients (Fig.
1B). This may be problem-specific, as the weight trajectories taken by pinn-
MDGA during training (Suppl. Fig. B.1) are fundamentally different from the
behavior of the other methods. This could be because the pinn-MGDA uses a
fundamentally different strategy purely based on satisfying Pareto-
stationarity conditions (Eq. 9). As shown in Suppl. Fig. B.2, this leads to a
rapid attenuation of the Fourier spectrum of the residual
$|\tilde{\bm{r}}(\bm{k})|$ over the first few learning epochs in this test
case, allowing pinn-MGDA to escape the stiffness defined in Proposition
II.2.1.

In summary, we observe stark differences between different weighting
strategies. Of the strategies compared, only the inverse-Dirichlet weighting
proposed here is able to avoid vanishing task-specific gradients, performing
on par with the $\epsilon$-optimal solution. The pinn-MGDA seems to be able to
escape, rather than avoid, the stiffness caused by vanishing task-specific
gradients. We also observe a clear correlation between having balanced
effective gradient distributions and achieving good approximation and
estimation accuracy.

Figure 2: Approximation accuracy and convergence for Sobolev training using
different weighting strategies. A) One instance of the function $\bm{u}$ given
by Eq. 20 visualized on a $128\times 128$ regular Cartesian grid. B) The
relative L2 error of the learned function approximation, i.e.
$\|\bm{u}_{\bm{\theta}}(\bm{x}_{i})-\bm{u}(\bm{x}_{i})\|_{2}/\|\bm{u}(\bm{x}_{i})\|_{2}$,
across all grid nodes $\bm{x}_{i}$ of the test data $\bm{u}(\bm{x}_{i})$ for
different objective weighting methods (symbols, see inset legend). C) The
relative L1 error in the estimated model coefficients, i.e.
$\|\bm{\xi}-\bm{\hat{\xi}}\|_{1}/\|\bm{\xi}\|_{1}$ for the different objective
weighting methods. D) The power spectrum $E(k)$ of the true function $\bm{u}$
(solid black line) compared with that of the learned estimators
$\bm{u}_{\bm{\theta}}$ for the different weighting methods (colored lines with
symbols). All results are shown after training for 20 000 epochs using the
Adam optimizer. The colored bands show to the standard deviation over 5
repetitions for different realizations (i.e., RNG seeds) of the random
function $\bm{u}$ (see Suppl. Eq. 20).

###  IV.2 Inverse-Dirichlet weighting enables multi-scale modeling using PINNs

We investigate the use of PINNs in modeling multi-scale dynamics in space and
time where vanishing task-specific gradients are a common problem. As a real-
world test problem, we consider meso-scale active turbulence as it occurs in
living fluids, such as bacterial suspensions [37, 38] and multi-cellular
tissues [39]. Unlike inertial turbulence, active turbulence has energy
injected at small length scales comparable to the size of the actively moving
entities, like the bacteria or individual motor proteins. The spatiotemporal
dynamics of active turbulence can be modeled using the generalized
incompressible Navier-Stokes equation [37, 38]:

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+\xi_{0}(\bm{u}\cdot\nabla\bm{u})$ $\displaystyle=-\nabla
p+\xi_{1}\nabla|\bm{u}|^{2}-\alpha\bm{u}-\beta|\bm{u}|^{2}\bm{u}+$
$\displaystyle\qquad\Gamma_{0}\Delta\bm{u}-\Gamma_{2}\Delta^{2}\bm{u}$ (11)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,$

where $\bm{u}(t,\bm{x})$ is the mean-field flow velocity and $p(t,\bm{x})$ the
pressure field. The coefficient $\xi_{0}$ scales the contribution from
advection and $\xi_{1}$ the active pressure. The $\alpha$ and $\beta$ terms
correspond to a quartic Landau-type velocity potential. The higher-order terms
with coefficients $\Gamma_{0}$ and $\Gamma_{2}$ model passive and active
stresses arising from hydrodynamic and steric interactions, respectively [38].
Supplementary Figure C.1 illustrates the rich multi-scale dynamics observed
when numerically solving this equation on a square for the parameter regime
$\Gamma_{0}<0$ and $\Gamma_{2}>0$ [37] (see Supplementary Material for details
on the numerical methods used). An important characteristic of Eq. 11 is the
presence of disparate length and time scales in addition to fourth-order
derivatives and nonlinearities in $\bm{u}$, which jointly lead to considerable
numerical stiffness [34]. Given such large dominant wavenumbers $\bm{k}_{0}$,
any data-driven approach, forward or inverse, is expected to encounter
significant learning pathologies as per Proposition II.2.1.

We demonstrate these pathologies by applying PINNs to approximate the forward
solution of Eq. 11 in 2D for the parameter regime explored in Ref. [37] in
both square (Fig. C.1) and annular (Fig. 3E) domains. The details of
simulation and training setup are given in the Supplementary Material. The
loss function includes four terms as given in Eq. 3, corresponding to the PDE
residual, initial and boundary conditions, and the divergence-freeness
constraint.

Figure 3A,B shows the average (over different initializations of the neural
network weights) relative L2 errors of the vorticity field over training
epochs when using different PINN weighting strategies. As predicted by
Proposition II.2.1, the uniformly weighted PINN fails completely both in the
square (Fig. 3A) and in the annular geometry (Fig. 3B). In this case,
optimization is entirely biased towards minimizing the PDE residual,
neglecting the initial, boundary, and incompressibility constraints. This
imbalance is also clearly visible in the loss gradient distributions shown in
Suppl. Fig. C.4A. Comparing the approximation accuracy after 7000 training
epochs confirms the findings from the previous section with inverse-Dirichlet
weighting outperforming the other methods, followed by pinn-MGDA and max/avg
weighting. This is made possible by assigning adequate weights for the initial
and boundary conditions in relative to the equation residual (See Suppl. Fig.
C.2A,D).

In Fig. 3C we compare the convergence of the learned approximation
$\bm{u}_{\bm{\theta}}$ to the ground-truth numerical solution for different
spatial resolutions $h$ of the discretization grid. The errors scale as
$O(h^{r})$ with convergence orders $r$ as identified by the dashed lines.
While the pinn-MGDA and max/avg weighting achieve linear convergence ($r=1$),
the inverse-Dirichlet weighted PINN converges with order $r=4$. The uniformly
weighted PINN does not converge at all and is therefore not shown in the plot.

Inverse-Dirichlet weighting also perfectly recovers the power spectrum $E(k)$
of the active turbulence velocity field, as shown in Fig. 3D. This confirms
that inverse-Dirichlet weighting enables PINNs to capture high frequencies in
the PDE solution, corresponding to fine structures in the fields, by avoiding
the stiffness problem from Proposition II.2.1. Indeed, the PINNs trained with
the other weighting methods fail to predict small structures in the flow
field, as shown by the spatial distributions of the prediction errors in Fig.
3F–I compared to the ground-truth numerical solution in Fig. 3E. The evolution
of the Fourier spectra $|\widetilde{r}(\bm{k})|$ of the residuals of both
velocity components $(u,v)=\bm{u}$ shown in Suppl. Fig. C.3 confirm the rapid
solution convergence achieved by inverse-Dirichlet weighted PINNs.

In the Supplementary Material, we provide results for two additional examples:
approximation the solution of the 2D Poisson equation on a square with a large
spectrum of frequencies, and learning the spatiotemporal dynamics of
curvature-driven level-set flow developing small structures. In both
supplementary multi-scale examples, we notice similar behavior with inverse-
Dirichlet weighting outperforming the other approaches either in terms of
accuracy or convergence rate.

In summary, our results suggest that avoiding vanishing task-specific
gradients can enable the use of PINNs in multi-scale problems previously not
amenable to neural-network modeling. Moreover, we observe that seemingly small
changes in the weighting can lead to significant differences in the
convergence order of the learned approximation.

Figure 3: Forward modeling of active turbulence using PINNs. A,B) Average
relative L2 prediction errors for the vorticity field $\omega$ of the active
turbulence PDE in Eq. 11 in square (A) and annular domains (B) using different
weighting strategies for PINN training (symbols, see inset legend). C)
Convergence of the different PINN predictions $\bm{u}_{\bm{\theta}}$ of the
flow vorticity to the ground-truth numerical solution in a square domain for
different grid resolutions. Dashed lines visualize integer convergence rates.
The PINN solution is extracted at epoch $7000$. The colored bands show the
standard deviations over three different initializations of the neural network
weights. D) Power spectrum $E(k)$ of the turbulent flow velocity field (solid
black line) compared with the power spectra of the approximations learned by
PINNs with different weighting (see inset legend) at epoch 500. E)
Visualization of the ground-truth vorticity ($\omega=\nabla\times\bm{u}$)
field in an annular domain obtained by numerically solving Eq. 11 as detailed
in the Supplementary Material. F–I) Spatial distribution of the relative
errors in the vorticity field when using PINNs with the weighting strategies
named in the panel titles.

###  IV.3 Inverse-Dirichlet weighting protects against catastrophic forgetting

Artificial neural networks tend to “forget” their parameterization when
training for different tasks sequentially [40], leading to a phenomenon
referred to as catastrophic forgetting. Sequential training is often used,
e.g., when learning to de-noise data and estimate parameters at the same time,
like combined image segmentation and denoising [41], or when learning
solutions of fluid mechanics models with additional constraints like
divergence-freeness of the flow velocity field [16].

In PINNs, catastrophic forgetting can occur, e.g., in the frequently practiced
approach of first training on measurement data alone and adding the PDE model
only at later epochs [42]. Such sequential training can be motivated by the
computational cost of evaluating the derivatives in the PDE using automatic
differentiation, so it seems prudent to pre-condition the network for data
fitting and introduce the residual later to regress for the missing network
parameters.

Figure 4: Catastrophic forgetting in sequential PINN training of active
turbulence. A): Test error in the velocity field $\bm{u}$ during the three-
step sequential training described in the main text with different weighting
strategies (colors). The solid vertical lines at epochs 1000 and 2000 mark the
introduction of the divergence constraint and PDE residual, respectively. The
dotted vertical lines mark epochs at which the learning rate $\eta(\tau)$ is
adjusted. B) Average relative L1 errors in the inferred values for the
coefficients $\xi_{0}$ ($\xi_{1}$ is absorbed into the effective pressure, see
Supplementary Material), $\alpha$, $\beta$, $\Gamma_{0}$, and $\Gamma_{2}$ of
Eq. 11 over training epochs. C) Average relative L2 error for inferring the
latent effective pressure gradient $\nabla p^{*}$ across training epochs.

Using the active turbulence test case described above, we investigate how the
different weighting schemes compared here influence the phenomenon of
catastrophic forgetting. For this, we consider the inverse problem of
inferring unknown model parameters and latent fields like the effective
pressure $p$ from flow velocity data $\bm{u}$. We therefore train a PINN in
three sequential steps: (1) first fitting flow velocity data, (2) later
including the additional objective for the incompressibility constraint
$\nabla\cdot\bm{u}=0$, (3) finally including also the objective for the PDE
model residual. This creates a setup in which catastrophic forgetting occurs
when using uniform PINN weights, as confirmed in Fig. 4A (arrows). As soon as
additional objectives (divergence-freeness at epoch 1000 and equation residual
at epoch 2000) are introduced, the PINN loses the parameterization learned on
the previous objectives, leading to large and increasing test error for the
field $\bm{u}$.

Dynamic weighting strategies can protect PINNs against catastrophic
forgetting, by adequately adjusting the weights whenever an additional
objective is introduced. This is confirmed in Fig. 4A with inverse-Dirichlet
and pinn-MGDA only minimally affected by catastrophic forgetting. The max/avg
weighting strategy [29], however, does not protect the network quite as well.
When looking at the test errors for the learned model coefficients (Fig. 4B)
and the accuracy of estimating the latent pressure field $p(\bm{x},t)$ (Fig.
4C), for which no training data is given, the inverse-Dirichlet weighting
proposed here outperforms the other approaches.

In summary, inverse-Dirichlet weighting can protect a PINN against
catastrophic forgetting, making sequential training a viable option. In our
tests, it offered the best protection of all compared weighting schemes.

##  V Conclusion and Discussion

Physics Informed Neural Networks (PINNs) have rapidly been adopted by the
scientific community for diverse applications in forward and inverse modeling.
It is straightforward to show that the training of a PINN amounts to a Multi-
Task Learning (MTL) problem, which is sensitive to the choice of
regularization weights. Trivial uniform weights are known to lead to failure
of PINNs, empirically found to exacerbate for physical models with disparate
scales or high derivatives.

As we have shown here, this failure can be explained by analogy to numerical
stiffness of Partial Differential Equations (PDEs) when understanding Sobolev
training of neural networks [7] as a special case of PINN training. This
connection enabled us to prove a proposition stating that PINN learning
pathologies are caused by dominant high-frequency components or high
derivatives in the physics prior. The asymptotic arguments in our proposition
inspired a novel dynamic adaptation strategy for the regularization weights of
a PINN during training, which we termed _inverse-Dirichlet weighting_.

We empirically compared inverse-Dirichlet weighting with the recently proposed
max/avg weighting for PINNs [29]. Moreover, we looked at PINN training from
the perspective of multi-objective optimization, which avoids a-priori choices
of regularization weights. We therefore adapted to PINNs the state-of-the-art
multiple gradient descent algorithm (MGDA) for training multi-objective
artificial neural networks, which has previously been proven to converge to a
Pareto-stationary solution [28]. This provided another baseline to compare
with. Finally, we compared against analytically derived $\epsilon$-optimal
static weights in a simple linear test problem where those can be derived,
providing an absolute baseline. In all comparisons, and across a wide range of
problems from 2D Poisson to multi-scale active turbulence and curvature-driven
level-set flow, inverse-Dirichlet weighting empirically performed best.

The inverse-Dirichlet weighting proposed here can be interpreted as an
uncertainty weighting with training uncertainty quantified by the batch
variance of the loss gradients. This is different from previous approaches
[26] that used weights quantifying the model uncertainty, rather than the
training uncertainty. The proposition proven here provides an explanation for
why the loss gradient variance may be a good choice, and our results confirmed
that weighting based on the Dirichlet energy of the loss objectives
outperforms other weighting heuristics, as well as Pareto-front seeking
algorithms like pinn-MGDA. The proposed inverse-Dirichlet weighting also
performed best at protecting a PINN against catastrophic forgetting in
sequential training settings where loss objectives are introduced one-by-one,
making sequential training of PINNs a viable option in practice.

While we have focussed on PINNs here, our analysis equally applies to Sobolev
training of other neural networks, and to other multi-objective losses, as
long as training is done using an explicit first-order stochastic gradient
descent (SGD) optimizer, of which Adam [32] is most commonly used. It would
also be interesting to see how the weighting strategies presented in this
study generalize to other first-order optimizers [43].

Taken together, we have exploited a connection between PINNs and Sobolev
training to explain a likely failure mode of PINNs, and we have proposed a
simple yet effective strategy to avoid training failure. The proposed inverse-
Dirichlet weighting strategy is easily included into any SDG-type optimizer at
almost no additional computational cost (see Suppl. Fig. C.8). It has the
potential to improve the accuracy and convergence of PINNs by orders of
magnitude, to enable sequential training with less risk of catastrophic
forgetting, and to extend the application realm of PINNs to multi-scale
problems that were previously not amenable to data-driven modeling.

###### Acknowledgements.

This work was supported by the German Research Foundation (DFG) – EXC-2068,
Cluster of Excellence “Physics of Life”, and by the Center for Scalable Data
Analytics and Artificial Intelligence (ScaDS.AI) Dresden/Leipzig, funded by
the Federal Ministry of Education and Research (BMBF).

## References

  * He  _et al._ [2016] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2016) pp. 770–778.
  * van den Oord  _et al._ [2016] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, Wavenet: A generative model for raw audio, in  _9th ISCA Speech Synthesis Workshop_ (2016) pp. 125–125.
  * Silver  _et al._ [2016] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, _et al._ , Mastering the game of go with deep neural networks and tree search, nature 529, 484 (2016).
  * Sirignano and Spiliopoulos [2018] J. Sirignano and K. Spiliopoulos, Dgm: A deep learning algorithm for solving partial differential equations, Journal of computational physics 375, 1339 (2018).
  * Raissi  _et al._ [2019a] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378, 686 (2019a).
  * Gühring  _et al._ [2020] I. Gühring, G. Kutyniok, and P. Petersen, Error bounds for approximations with deep relu neural networks in w s, p norms, Analysis and Applications 18, 803 (2020).
  * Czarnecki  _et al._ [2017] W. M. Czarnecki, S. Osindero, M. Jaderberg, G. Swirszcz, and R. Pascanu, Sobolev training for neural networks, in  _Advances in Neural Information Processing Systems_ (2017) pp. 4278–4287.
  * Raissi  _et al._ [2019b] M. Raissi, Z. Wang, M. S. Triantafyllou, and G. E. Karniadakis, Deep learning of vortex-induced vibrations, Journal of Fluid Mechanics 861, 119 (2019b).
  * Fang and Zhan [2019] Z. Fang and J. Zhan, Deep physical informed neural networks for metamaterial design, IEEE Access 8, 24506 (2019).
  * Liu and Wang [2019] D. Liu and Y. Wang, Multi-fidelity physics-constrained neural network and its application in materials modeling, Journal of Mechanical Design 141 (2019).
  * Chen  _et al._ [2020] Y. Chen, L. Lu, G. E. Karniadakis, and L. Dal Negro, Physics-informed neural networks for inverse problems in nano-optics and metamaterials, Optics express 28, 11618 (2020).
  * Raissi  _et al._ [2020] M. Raissi, A. Yazdani, and G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367, 1026 (2020).
  * Sun  _et al._ [2020] L. Sun, H. Gao, S. Pan, and J.-X. Wang, Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data, Computer Methods in Applied Mechanics and Engineering 361, 112732 (2020).
  * Mao  _et al._ [2020] Z. Mao, A. D. Jagtap, and G. E. Karniadakis, Physics-informed neural networks for high-speed flows, Computer Methods in Applied Mechanics and Engineering 360, 112789 (2020).
  * Raissi  _et al._ [2019c] M. Raissi, H. Babaee, and P. Givi, Deep learning of turbulent scalar mixing, Physical Review Fluids 4, 124501 (2019c).
  * Jin  _et al._ [2020] X. Jin, S. Cai, H. Li, and G. E. Karniadakis, Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations, Journal of Computational Physics 426, 109951 (2020).
  * Yazdani  _et al._ [2020] A. Yazdani, L. Lu, M. Raissi, and G. E. Karniadakis, Systems biology informed deep learning for inferring parameters and hidden dynamics, PLoS computational biology 16, e1007575 (2020).
  * Sahli Costabal  _et al._ [2020] F. Sahli Costabal, Y. Yang, P. Perdikaris, D. E. Hurtado, and E. Kuhl, Physics-informed neural networks for cardiac activation mapping, Frontiers in Physics 8, 42 (2020).
  * Kissas  _et al._ [2020] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, and P. Perdikaris, Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 358, 112623 (2020).
  * Zheng  _et al._ [2020] Q. Zheng, L. Zeng, and G. E. Karniadakis, Physics-informed semantic inpainting: Application to geostatistical modeling, Journal of Computational Physics 419, 109676 (2020).
  * Rao  _et al._ [2021] C. Rao, H. Sun, and Y. Liu, Physics-informed deep learning for computational elastodynamics without labeled data, Journal of Engineering Mechanics 147, 04021043 (2021).
  * Yang  _et al._ [2020] L. Yang, X. Meng, and G. E. Karniadakis, B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data, Journal of Computational Physics 425, 109913 (2020).
  * Zhang  _et al._ [2020] D. Zhang, L. Guo, and G. E. Karniadakis, Learning in modal space: Solving time-dependent stochastic pdes using physics-informed neural networks, SIAM Journal on Scientific Computing 42, A639 (2020).
  * Han  _et al._ [2018] J. Han, A. Jentzen, and E. Weinan, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences 115, 8505 (2018).
  * Pang  _et al._ [2019] G. Pang, L. Lu, and G. E. Karniadakis, fpinns: Fractional physics-informed neural networks, SIAM Journal on Scientific Computing 41, A2603 (2019).
  * Kendall  _et al._ [2018] A. Kendall, Y. Gal, and R. Cipolla, Multi-task learning using uncertainty to weigh losses for scene geometry and semantics, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2018) pp. 7482–7491.
  * Chen  _et al._ [2018] Z. Chen, V. Badrinarayanan, C.-Y. Lee, and A. Rabinovich, Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks, in  _International Conference on Machine Learning_ (PMLR, 2018) pp. 794–803.
  * Sener and Koltun [2018] O. Sener and V. Koltun, Multi-task learning as multi-objective optimization, Advances in Neural Information Processing Systems 31, 527 (2018).
  * Wang  _et al._ [2020] S. Wang, Y. Teng, and P. Perdikaris, Understanding and mitigating gradient pathologies in physics-informed neural networks, arXiv preprint arXiv:2001.04536 (2020).
  * Zitzler and Thiele [1999] E. Zitzler and L. Thiele, Multiobjective evolutionary algorithms: a comparative case study and the strength pareto approach, IEEE Transactions on Evolutionary Computation 3, 257 (1999).
  * Boyd  _et al._ [2004] S. Boyd, S. P. Boyd, and L. Vandenberghe,  _Convex optimization_ (Cambridge university press, 2004).
  * Kingma and Ba [2015] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, in  _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_, edited by Y. Bengio and Y. LeCun (2015).
  * Glorot and Bengio [2010a] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in  _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ (JMLR Workshop and Conference Proceedings, 2010) pp. 249–256.
  * Cox and Matthews [2002] S. M. Cox and P. C. Matthews, Exponential time differencing for stiff systems, Journal of Computational Physics 176, 430 (2002).
  * Désidéri [2012] J.-A. Désidéri, Multiple-gradient descent algorithm (mgda) for multiobjective optimization, Comptes Rendus Mathematique 350, 313 (2012).
  * van der Meer  _et al._ [2020] R. van der Meer, C. Oosterlee, and A. Borovykh, Optimally weighted loss functions for solving pdes with neural networks, arXiv preprint arXiv:2002.06269 (2020).
  * Wensink  _et al._ [2012] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulence in living fluids, Proceedings of the National Academy of Sciences 109, 14308 (2012).
  * Dunkel  _et al._ [2013] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Fluid dynamics of bacterial turbulence, Physical review letters 110, 228102 (2013).
  * Ramaswamy and Jülicher [2016] R. Ramaswamy and F. Jülicher, Activity induces traveling waves, vortices and spatiotemporal chaos in a model actomyosin layer, Scientific reports 6, 1 (2016).
  * Kirkpatrick  _et al._ [2017] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska, _et al._ , Overcoming catastrophic forgetting in neural networks, Proceedings of the national academy of sciences 114, 3521 (2017).
  * Buchholz  _et al._ [2020] T.-O. Buchholz, M. Prakash, D. Schmidt, A. Krull, and F. Jug, Denoiseg: joint denoising and segmentation, in  _European Conference on Computer Vision_ (Springer, 2020) pp. 324–337.
  * Both  _et al._ [2021] G.-J. Both, S. Choudhury, P. Sens, and R. Kusters, Deepmod: Deep learning for model discovery in noisy data, Journal of Computational Physics 428, 109985 (2021).
  * Schmidt  _et al._ [2020] R. M. Schmidt, F. Schneider, and P. Hennig, Descending through a crowded valley–benchmarking deep learning optimizers, arXiv preprint arXiv:2007.01547 (2020).
  * Hardy and Titchmarsh [1931] G. Hardy and E. Titchmarsh, A note on parseval’s theorem for fourier transforms, Journal of the London Mathematical Society 1, 44 (1931).
  * Rahaman  _et al._ [2019] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville, On the spectral bias of neural networks, in  _International Conference on Machine Learning_ (PMLR, 2019) pp. 5301–5310.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Kassam and Trefethen [2005] A.-K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff pdes, SIAM Journal on Scientific Computing 26, 1214 (2005).
  * Osher and Fedkiw [2006] S. Osher and R. Fedkiw,  _Level set methods and dynamic implicit surfaces_ , Vol. 153 (Springer Science & Business Media, 2006).
  * Malladi and Sethian [1995] R. Malladi and J. A. Sethian, Image processing via level set curvature flow, proceedings of the National Academy of sciences 92, 7046 (1995).
  * Zhao  _et al._ [1996] H.-K. Zhao, T. Chan, B. Merriman, and S. Osher, A variational level set approach to multiphase motion, Journal of computational physics 127, 179 (1996).
  * Bergdorf  _et al._ [2010] M. Bergdorf, I. F. Sbalzarini, and P. Koumoutsakos, A lagrangian particle method for reaction–diffusion systems on deforming surfaces, Journal of mathematical biology 61, 649 (2010).
  * Alvarez  _et al._ [1993] L. Alvarez, F. Guichard, P.-L. Lions, and J.-M. Morel, Axioms and fundamental equations of image processing, Archive for rational mechanics and analysis 123, 199 (1993).
  * Alamé  _et al._ [2020] K. Alamé, S. Anantharamu, and K. Mahesh, A variational level set methodology without reinitialization for the prediction of equilibrium interfaces over arbitrary solid surfaces, Journal of Computational Physics 406, 109184 (2020).
  * Glorot and Bengio [2010b] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in  _In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS10). Society for Artificial Intelligence and Statistics_ (2010).

##  Appendix A Material and Methods

###  A.1 Determining analytical optimal weights

For Sobolev loss functions based on elliptic PDEs, optimal weights $\lambda$
that minimize all objectives in an unbiased way can be determined by
$\epsilon$-closeness analysis [36].

######  Definition 1 ($\epsilon$-closeness)

A candidate solution $\bm{u}$ is called $\epsilon$-close to the true solution
$\hat{\bm{u}}$ if it satisfies

$\frac{|\mathcal{D}^{\bm{\beta}}\bm{u}(\mathbf{x})-\mathcal{D}^{\bm{\beta}}\hat{\bm{u}}(\mathbf{x})|}{|\mathcal{D}^{\bm{\beta}}\hat{\bm{u}}(\mathbf{x})|}\leq\epsilon$
(12)

for any multi-index $\bm{\beta}\in\mathbb{N}$ and
$\mathbf{x}\in\Omega\subset\mathbb{R}^{d}$, thus
$\mathcal{D}^{\bm{\beta}}\bm{u}(\mathbf{x})=\frac{\partial^{|\bm{\beta}|}\bm{u}(\textbf{x})}{\partial
x_{1}^{\beta_{1}}\partial x_{2}^{\beta_{2}}\ldots\partial x_{d}^{\beta_{d}}}$
with $|\bm{\beta}|=\sum_{i=1}^{d}\beta_{i}$.

For the loss function typically used in Sobolev training of neural networks,
i.e.,

$\mathcal{L}_{k}(\cdot)=\int_{\Omega}\left(\mathcal{D}_{\mathbf{x}}^{k}\hat{\bm{u}}-\mathcal{D}_{\mathbf{x}}^{k}\bm{u}\right)^{2}\mathrm{d}\bm{x}$

we can use $\epsilon$-closeness analysis to show that

$\mathcal{L}_{k}\leq\epsilon^{2}\int_{\Omega}\left(|\xi_{k}\mathcal{D}_{\bm{x}}^{k}(\bm{\hat{u}})|\right)^{2}\mathrm{d}\bm{x}\,.$
(13)

Using this inequality, we can bound the total loss as

$\mathcal{L}(\bm{u})\leq\epsilon^{2}\sum_{k=1}^{K}\lambda_{k}\int_{\Omega}\left(|\xi_{k}\mathcal{D}_{\bm{x}}^{k}(\bm{\hat{u}})|\right)^{2}\mathrm{d}\bm{x}=\epsilon^{2}\sum_{k=1}^{K}\lambda_{k}\mathcal{I}_{k}$
(14)

with
$\mathcal{I}_{k}=\int_{\Omega}\left(|\xi_{k}\mathcal{D}_{\bm{x}}^{k}(\bm{\hat{u}})|\right)^{2}\\!\mathrm{d}\Omega$.
This bound can be used to determine the weights $\\{\lambda_{k}\\}_{k=1}^{K}$
that result in the smallest $\epsilon$:

$\displaystyle\mathcal{L}_{k}(\bm{u})$
$\displaystyle\leq\frac{1}{\lambda_{k}}\mathcal{L}(\bm{u})\leq\epsilon^{2}\lambda_{k}\mathcal{I}_{k}\>\>\forall
k={1,\ldots,K},$ $\displaystyle\mathcal{L}(u)$
$\displaystyle\leq\epsilon^{2}\min\\{\lambda_{1}\mathcal{I}_{1},\ldots,\lambda_{K}\mathcal{I}_{K}\\}.$
(15)

The weights $\lambda_{k},\forall k\in{1,...,K}$ that lead to the tightest
upper bound on the total loss in Eq. A.1 can be found by solving the
maximization-minimization problem described in following corollary:

###### Corollary A.1.1

A maximization-minimization problem can be transformed into a piecewise linear
convex optimization problem by introducing an auxiliary variable $z$, such
that

$\displaystyle\textrm{maximize}\>\>z\qquad\quad$ $\displaystyle\textrm{subject
to }\mathbf{1}^{\top}\bm{\lambda}^{*}=1;$ $\displaystyle
z\leq\lambda_{k}\mathcal{I}_{k},\quad k={1,\ldots,K}.$

The analytical solution to this problem exists and is given by
$\lambda_{k}^{*}=\frac{\Pi_{j\neq k}\mathcal{I}_{j}}{\sum_{k=1}^{K}\Pi_{j\neq
k}\mathcal{I}_{j}},\quad k=1,\ldots,K$.

For a Sobolev training problem with loss given by Eq. 5, weights that lead to
$\epsilon$-optimal solutions are computed using

$\mathcal{I}_{k}=\int_{\Omega}\left(|\xi_{k}\mathcal{D}_{\bm{x}}^{k}(\bm{\hat{u}})|\right)^{2}\\!\mathrm{d}\bm{x}.$

###  A.2 PINN Multiple Gradient Descent Algorithm (pinn-MGDA)

We adapt the MGDA to PINNs by formulating the conditions for Pareto stationary
as stated in Eq. 9 as a quadratic program:

$\underset{\bm{\lambda}\in\mathbb{R}^{K}}{\text{minimize
}}\frac{1}{2}\bm{\lambda}^{\top}\bm{Q}\bm{\lambda},\textrm{ s.t.
}\lambda_{k}\geq 0\>\forall k\in[K]\textrm{ and }\sum_{k=1}^{K}\lambda_{k}=1,$
(16)

where $\bm{Q}=\bm{U}^{\top}\bm{U}$. The matrix
$\bm{U}\in\mathbb{R}^{|\textrm{B}|\times K}$ is given by

$\bm{U}=\begin{bmatrix}\vdots\>\>\qquad\qquad\vdots\qquad\qquad\>\>\vdots\>\>\\\
\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{1}\>\>\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{2}\>\>\ldots\>\>\nabla_{\bm{\theta}^{sh}}\mathcal{L}_{K}\\\
\vdots\>\>\qquad\qquad\vdots\qquad\qquad\>\>\vdots\>\>\\\ \end{bmatrix}.$

We solve the quadratic program in Eq. 16 using the Frank-Wolfe algorithm.

Figure A.1:  Multiple gradient descent algorithm  The schematic of multi-
objective optimization for minimizing two objectives $f_{1}(x)$ and
$f_{2}(x)$. The red curve corresponds to the Pareto-optimal front. The MGDA
algorithm searches for the descent direction (arrow) that minimizes all the
objectives and is closer to the Pareto-optimal front. The grey region
corresponds to the feasible objective space.

##  Appendix B Sobolev Training

### Proof of Proposition 1

Let us assume the Fourier spectrum of the signal is a delta function, i.e.
$\widetilde{u}(k)=\Gamma_{0}\delta(k-k_{0})$, with $k_{0}$ being the dominant
wavenumber. The residual at training time $(t)$ is given as
$r(\cdot,t)=u_{\theta}-u$, and its corresponding Fourier coefficients for the
wavenumber $k$ is given by
$\widetilde{r}(k,t)=\widetilde{u}_{\theta}(k,t)-\widetilde{u}(k)$. The
learning rate of the residual for the dominant frequency is given as,

$\displaystyle\Big{|}\frac{\partial\widetilde{r}(k_{0},\tau)}{\partial\tau}\Big{|}$
$\displaystyle=\Big{|}\frac{\partial\widetilde{u}_{\theta}(k_{0},\tau)}{\partial\theta}\Big{|}\Big{|}\frac{\partial\theta}{\partial\tau}\Big{|}$
$\displaystyle=\Big{|}\eta\frac{\partial\widetilde{u}_{\theta}(k_{0})}{\partial\theta}\Big{|}\Big{|}\frac{\partial\mathcal{L}}{\partial\theta}\Big{|}$
$\displaystyle=\Big{|}\eta\frac{\partial\widetilde{u}_{\theta}(k_{0})}{\partial\theta}\Big{|}\Big{|}\lambda_{1}\frac{\partial\mathcal{L}_{1}}{\partial\theta}+\lambda_{2}\frac{\partial\mathcal{L}_{2}}{\partial\theta}\Big{|}$
(17)

From Parseval’s theorem [44], we can compute the magnitude of the gradients
for both objectives $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ with respect to
the network parameters $\theta$.

$\displaystyle\frac{\partial\mathcal{L}_{1}}{\partial\theta}$
$\displaystyle=2\sum_{k=-N/2}^{N/2-1}\Big{|}\textbf{Re}\Big{[}\widetilde{u}_{\theta}(k)-\widetilde{u}(k)\Big{]}\frac{\partial\widetilde{u_{\theta}}(k)}{\partial\theta}\Big{|}$
$\displaystyle\leq
2\Big{|}\Gamma_{0}\frac{\partial\widetilde{u}_{\theta}(k_{0})}{\partial\theta}\Big{|}+2\sum_{k=-N/2}^{N/2-1}\Big{|}\widetilde{u}_{\theta}(k)\frac{\partial\widetilde{u}_{\theta}(k)}{\partial\theta}\Big{|}$
$\displaystyle\approx\mathcal{O}(1)\Big{|}\frac{\partial\widetilde{r}(k_{0})}{\partial\theta}\Big{|}$
(18)

The second step in the above calculation involves triangular inequality and
the assumption that the Fourier spectrum of the signal is a delta function. In
a similar fashion, we can repeat the same analysis for the Sobolev part of the
objective function as follows,

$\displaystyle\frac{\partial\mathcal{L}_{2}}{\partial\theta}$
$\displaystyle=2\sum_{k=-N/2}^{N/2-1}\Big{|}\textbf{Re}\Big{[}(ik)^{m}\big{(}\widetilde{u}_{\theta}(k)-\widetilde{u}(k)\big{)}\Big{]}\frac{\partial(ik)^{m}\widetilde{u}_{\theta}(k)}{\partial\theta}\Big{|}$
$\displaystyle=2\sum_{k=-N/2}^{N/2-1}\Big{|}\textbf{Re}\Big{[}(ik)^{2m}\big{(}\widetilde{u}_{\theta}(k)-\widetilde{u}(k)\big{)}\Big{]}\frac{\partial\widetilde{u}_{\theta}(k)}{\partial\theta}\Big{|}$
$\displaystyle\leq
2\Big{|}\Gamma_{0}\textbf{Re}\big{[}(ik_{0})^{2m}\big{]}\frac{\partial\widetilde{u}_{\theta}(k_{0})}{\partial\theta}\Big{|}+2\sum_{k=-N/2}^{N/2-1}\Big{|}\textbf{Re}\big{[}(ik)^{2m}\big{]}\widetilde{u}_{\theta}(k)\frac{\partial\widetilde{u}_{\theta}(k)}{\partial\theta}\Big{|}$
$\displaystyle\approx\mathcal{O}(k_{0}^{2m})\Big{|}\frac{\partial\widetilde{r}(k_{0})}{\partial\theta}\Big{|}$
(19)

By substituting equations (B) and (B) into equation (B), we get,

$\displaystyle\Big{|}\frac{\partial\widetilde{r}(k_{0},\tau)}{\partial\tau}\Big{|}$
$\displaystyle=\Big{|}\eta\frac{\partial\widetilde{u}_{\theta}(k_{0})}{\partial\theta}\Big{|}\left(\lambda_{1}\mathcal{O}(1)\Big{|}\frac{\partial\widetilde{r}(k_{0})}{\partial\theta}\Big{|}+\lambda_{2}\mathcal{O}(k_{0}^{2m})\Big{|}\frac{\partial\widetilde{r}(k_{0})}{\partial\theta}\Big{|}\right)$
$\displaystyle=\frac{\eta}{\mathcal{O}(k_{0})}\Big{[}\lambda_{1}\mathcal{O}(1)+\lambda_{2}\mathcal{O}(k_{0}^{2m})\Big{]}\Big{|}\frac{\partial\widetilde{r}(k_{0})}{\partial\theta}\Big{|}$

We use the fact that the spectral rate of the network function residual is
inversely proportional to the wavenumber, i.e.
$\partial\widetilde{u}_{\theta}(k)/\partial\theta=\mathcal{O}\left(k^{-1}\right)$
[45].

### Sobolev training set-up

We use a neural network $f_{\theta}:(x,y)\rightarrow(u)$ with 5 layers and 64
neurons per layer and $\sin(\cdot)$ activation function. The target function
is a sinusoidal forcing term of the form,

$u(x,y)=\sum_{i=1}^{M}A_{x}^{i}\cos(2\pi
l_{x}^{i}x/L+\phi_{x}^{i})A_{y}^{i}\sin(2\pi l_{y}^{i}y/L+\phi_{y}^{i}),$ (20)

where $M=20$ and $L=2\pi$ is the size of the domain. The parameters
$A_{x},A_{y}$ are drawn independently and uniformly from the interval
$[-5,5]$, and similarly $\phi_{x},\phi_{y}$ are drawn from the interval
$[0,2\pi]$. The parameters $l_{x},l_{y}$ which set the local length scales are
sampled uniformly from the set $\\{1,2,3,4,5\\}$. The target function $u(x,y)$
is given on a 2D spatial grid of resolution $128\times 128$ covering the
domain $[0,2\pi]\times[0,2\pi]$. The grid data is then divided into a $50/50$
train and test split which corresponds to $8192$ training points used. We
train the neural network for 20000 epochs using the Adam optimizer [46]. Each
epoch is composed of 2 mini-batch updates of batch size $|\textrm{B}|=4096$.
The initial leaning rate is chosen as $\eta=10^{-3}$ and decreased by a factor
of $10$ after $10000$ and $15000$ epochs.

Figure B.1: $\lambda$ trajectories for different weighing strategies  The
dynamic weights $\lambda$ over the training epochs for a single instance of
Sobolev training of the neural networks for different methods A) pinn-MGDA B)
$\epsilon$-optimal C) inverse-Dirichlet D) max/avg. Different trajectories in
the same plots correspond to the weights $\lambda_{k},k\in{1,..,K}$,
associated with different objectives. On closer inspective we can notice that
the relative weights or relative importance measure
$(\lambda_{i}/\lambda_{j}):\forall{i,j}\in{1,..,K}$ between objectives for
inverse-Dirichlet and the $\epsilon$-optimal strategy are approximately of the
same order: $\lambda_{0}/\lambda_{1}\approx 10,\lambda_{0}/\lambda_{2}\approx
10^{2}-10^{3},\lambda_{0}/\lambda_{3}\approx
10^{3}-10^{4},\lambda_{0}/\lambda_{4}\approx 10^{4}-10^{5}$ (shown in red
double-headed arrows).

Figure B.2: Evolution of the Fourier spectrum of the function residual
$\widetilde{r}(\bm{k})$ The Fourier spectrum of the residual tracked over the
training epochs for different methods (shown in inset). The color intensity
corresponds to the magnitude of the residual, i.e. $|\widetilde{r}(\bm{k})|$.

##  Appendix C Forward and inverse modeling in active turbulence

The ground truth solution of the 2D incompressible active fluid is computed in
the vorticity formulation which is given as follows,

$\frac{\partial\omega}{\partial
t}+\xi_{0}\bm{u}\cdot\nabla\omega=-(\alpha+\beta|\bm{u}|^{2})\omega+\beta\nabla|\bm{u}|^{2}\times\bm{u}+\Gamma_{0}\Delta\omega-\Gamma_{2}\Delta^{2}\omega.$
(21)

The simulation is conducted with a pseudo-spectral solver, where the
derivatives are computed in the Fourier space and the time integration is done
using the Integration factor method [34, 47]. The equations are solved in the
domain $[0,2\pi]\times[0,2\pi]$ with a spatial resolution of $256\times 256$
and a time-step $dt=0.0001$. The incompressiblity is enforced by projecting
the intermediate velocities on to the divergence-free manifold. The
simulations in the square domain are performed with periodic boundary
conditions, where the values of the parameters are chosen as follows:
$\xi_{0}=3.5,\alpha=-1,\beta=0.5,\Gamma_{0}=-0.045,\Gamma_{2}=|\Gamma|^{3}$
[37].

### Forward problem set-up

We study the forward solution of the active turbulence problem using PINNs in
the $(\omega-\psi)$ formulation. We employ a neural network
$f_{\theta}:(x,y,t)\rightarrow(u,v)$ with 7 layers and 100 neurons per layer
and $\sin$-activation functions. The vorticity $\omega$ is directly computed
from the network outputs using automatic differentiation. For the solution in
the square domain, we select a sub-domain
$[\frac{\pi}{4},\frac{3\pi}{4}]\times[\frac{\pi}{4},\frac{3\pi}{4}]$ with 50
time steps ($dt=0.01$) accounting to a total time of $0.5$ units in simulation
time. We choose $N_{I}=4096$ points for the initial conditions, $N_{B}=12800$
boundary points and a total of $N_{\text{int}}=192200$ residual points in the
interior of the square domain. The network is trained for 7000 epochs using
the Adam optimizer with each epoch consisting of $49$ mini-batch updates with
a batch size $|\textrm{B}|=4000$. The velocities are provided at the
boundaries of the training domain.

For the forward simulation in the annular geometry, the domain is chosen with
the inner radius of $r_{\text{inner}}\approx 0.245$ units and the outer radius
$r_{\text{outer}}\approx 0.85$ units. We choose $N_{I}=3536$ points for the
initial conditions, $N_{B}=14400$ boundary points and a total of
$N_{\text{int}}=162400$ residual points in the interior of the domain. The
network training is again performed over 7000 epochs using Adam with each
epoch consisting of $41$ mini-batch updates with a batch size
$|\textrm{B}|=4000$. Dirichlet boundary conditions are enforced at the inner
and outer rims of the annular geometry by providing the velocities. The
network is optimized subject to the loss function

$\mathcal{L}(\bm{u},\theta)=\lambda_{1}\mathcal{L}_{\text{boundary}}+\lambda_{2}\mathcal{L}_{\text{res}}+\lambda_{3}\mathcal{L}_{\text{div}}+\lambda_{4}\mathcal{L}_{\text{init}}.$

Here, $\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}$ correspond to the
weights of the objectives for boundary condition, equation residual,
divergence, and initial condition, respectively.

Figure C.1:  Training domain for active turbulence problem  Left) Training
domain (inside the dotted square) used for forward solution in the square
domain. Right) Training domain (interior of the dotted annular outline) used
for forward solution in the annular geometry.

### Inverse problem set-up

For studying the inverse problem we resort to the primitive formulation
$(u,v,p^{*})$ of the active turbulence model, with $p^{*}=-\nabla
p+\xi_{1}\nabla|\bm{u}|^{2}$ as the effective pressure. We used a neural
network $f_{\theta}:(x,y,t)\rightarrow(u,v,p^{*})$ with 5 layers and 100
neurons per layer and $\sin$-activation functions. Given the measurement data
of flow velocities $\\{u_{i},v_{i}\\}_{i=1}^{N}$ at discrete points, we infer
the effective pressure $p^{*}$ along with the parameters
$\Sigma=\\{\xi_{0},\alpha,\beta,\Gamma_{0},\Gamma_{2}\\}$. The total of
$204800$ points, sampled from the same domain as for the forward solution in
the square domain, are utilized with a $70/30$ train and test split, which
corresponds to $143360$ training points used. The network is trained for a
total of $7000$ epochs using the Adam optimizer with each epoch consisting of
$35$ mini-batch updates with a batch-size $|\textrm{B}|=4096$. We optimize the
model subject to the loss function

$\mathcal{L}(\bm{u},\theta)=\lambda_{1}\mathcal{L}_{\text{fit}}+\lambda_{2}\mathcal{L}_{\text{res}}+\lambda_{3}\mathcal{L}_{\text{div}},$

with $\lambda_{1},\lambda_{2},\lambda_{3}$ corresponding to the data-fidelity,
equation residual and the objective that penalises non-zero divergence in
velocities, respectively. The training is done in a three-step approach, with
the first $1000$ epochs used for fitting the measurement data (pre-training,
$\lambda_{1}=1,\lambda_{2}=\lambda_{3}=0$), the next $1000$ epochs training
with an additional constraint on the divergence ($\lambda_{1}\neq
0,\lambda_{2}\neq 0,\lambda_{3}=0$), and finally followed by the introduction
of the equation residual ($\lambda_{1}\neq 0,\lambda_{2}\neq 0,\lambda_{3}\neq
0$). At every checkpoint introducing a new task, we reset all $\lambda_{k}$
for the computation of $\hat{\lambda}_{k}$ (see Eq. 7 main text) of the
max/avg variant to $\lambda_{k}=1$. For both, forward and inverse problems, we
choose an initial learning rate $\eta=10^{-3}$ and decrease it by a factor of
$10$ after $3000$ and $6000$ epochs (See Fig. C.2).

Figure C.2: $\lambda$ trajectories for different weighting strategies for the
active turbulence problem. A, B, C) correspond to the $\lambda$ trajectories
for the forward solution of the active turbulence problem in the square domain
(See left Fig. C.1). D,E,F) correspond to the $\lambda$ trajectories for the
forward solution in the annular geometry (See right Fig. C.1). The columns
correspond to the different weighting strategies as named on the top left
corner of each column.  Figure C.3: Evolution of Fourier spectrum of the
function residual $\widetilde{r}(\bm{k})$.  Top row) The Fourier spectrum of
the residual for the horizontal velocity component $u$. Bottom row) The
Fourier spectrum of the residual for the vertical velocity component $v$. The
columns correspond to different methods (shown in the inset). The color
intensity corresponds to the magnitude of the residual, i.e.
$|\widetilde{r}(\bm{k})|$.  Figure C.4:  Gradient histograms illustrating
vanishing taks-specific gradients for the active turbulence problem. The
gradient distributions of the different objectives (insets on the right) at
training epoch $5000$ for the forward problem in a squared domain. The panels
correspond to different $\lambda$ weighting schemes as named. Figure C.5:
Catastrophic forgetting in sequential PINN training of active turbulence from
noisy measurement data. A): Test error in the velocity field $\bm{u}$ during
the three-step sequential training described in the main text with different
weighting strategies (colors). The solid vertical lines at epochs 1000 and
2000 mark the introduction of the divergence constraint and PDE residual,
respectively. The dotted vertical lines mark epochs at which the learning rate
$\eta(\tau)$ is adjusted. B) Average relative L1 errors in the inferred values
for the coefficients $\xi_{0}$ ($\xi_{1}$ is absorbed into the effective
pressure, see Supplementary Material), $\alpha$, $\beta$, $\Gamma_{0}$, and
$\Gamma_{2}$ of Eq. 11 over training epochs. C) Average relative L2 error for
inferring the latent pressure gradient $\nabla p$ across training epochs. The
simulated data $\bm{u}$ is corrupted with additive Gaussian noise
$\hat{\bm{u}}=\bm{u}+0.25\bm{\epsilon}$, where
$\bm{\epsilon}\sim\mathcal{G}(0,\psi^{2})$ is a vector of elementwise
independent and identically distributed Gaussian random numbers with mean zero
and empirical variance $\psi^{2}=\textrm{Var}\\{u_{1},...,u_{N}\\}$.  Figure
C.6: Inferred effective pressure. The top row corresponds to the prediction of
horizontal component of the effective pressure gradient $\partial_{x}p^{*}$
given the velocity field $\bm{u}$. The bottom row shows the corresponding
point-wise relative error for the prediction. Different columns correspond to
different weighting strategies as per the titles.  Figure C.7: Inferred
effective pressure from noisy measurement data. The top row corresponds to the
prediction of horizontal component of the effective pressure gradient
$\partial_{x}p^{*}$ given the noise corrupted velocity field
$\hat{\bm{u}}=\bm{u}+0.25\bm{\epsilon}$, where
$\bm{\epsilon}\sim\mathcal{G}(0,\psi^{2})$ is a vector of elementwise
independent and identically distributed Gaussian random numbers with mean zero
and empirical variance $\psi^{2}=\textrm{Var}\\{u_{1},...,u_{N}\\}$. The
bottom row shows the corresponding point-wise relative error for the
prediction. Different columns correspond to different weighting strategies
(titles).

Figure C.8:  Computational aspects of dynamic weighing for the active
turbulence problem  Left) Training times per epoch for the inverse problem in
the primitive formulation with $35$ mini-batch updates per epoch of batch size
$|\textrm{B}|=4096$. Right) Training times per epoch for the forward problem
in the vorticity-streamline formulation with $49$ mini-batch updates per epoch
of batch size $|\textrm{B}|=4000$. The extra computational cost of the forward
problem comes from using more training points and from the additional
derivative that needs to be computed for evaluating the vorticity. The
profiling was done on a NVIDIA GTX 1080 graphics card equipped with 8GB VRAM.

##  Appendix D 2D Poisson problem

We consider the problem of solving the 2D Poisson equation given as,

$\Delta u=-2\omega^{2}\cos(\omega x)\sin(\omega y),$ (22)

with the parameter $\omega$ controlling the frequency of the oscillatory
solution $u(x,y)$. The analytical solution $u(x,y)=\cos(\omega x)\sin(\omega
y)$ is computed on the domain $[0,1]\times[0,1]$ with a grid resolution of
$100\times 100$. The boundary condition is given as
$\mathcal{B}(x,y)=\cos{(\omega x)}\sin{(\omega y)}$, $x,y\in\partial\Omega$.
We set up the PINN with two objective functions one for handling boundary
conditions and the other to compute the residual with the loss

$\mathcal{L}(u,\theta)=\lambda_{1}\mathcal{L}_{\text{res}}+\lambda_{2}\mathcal{L}_{\text{boundary}}$
(23)

We set $N_{\mathcal{B}}=400$ sampled uniformly on the four boundaries, and
randomly sample $N_{\text{int}}=2500$ points from the interior of the domain.
We choose a network $f_{\theta}:(x,y)\rightarrow(u)$ with $5$ layers and $50$
nodes per layer and the $\tanh$ activation function is used. The neural
network is trained for $30000$ epochs using the Adam optimizer. We choose the
initial learning rate $\eta=10^{-3}$ and decrease it by a factor of $10$ after
$10000$, $20000$ and $30000$ epochs. In Fig. D.1 we show the relative L2 error
for different weighting strategies across changing wavenumber $\omega$. We
notice a very similar trend as observed in the previous problems, with inverse
Dirichlet scheme performing in par with the analytical $\epsilon$ optimal
weighting strategy, followed by max/avg and pinn-MGDA.

Figure D.1: Comparison between methods for solution of Poisson Equation  Left:
Solution of the Poisson equation (Eq. 22) for the wave-number $\omega=6$.
Comparison between relative L2 error for different weighing strategies for
increasing wave-numbers $(\omega)$.

##  Appendix E Curvature-driven level-set

The level-set method [48] is a popular method used for implicit representation
of surfaces. In contrast to explicit methods for interface capturing, level-
set methods offer a natural means to handle complex changes in the interface
topology. For this reason, level-set methods have found extensive applications
in computer vision [49], multiphase and computational fluid dynamics [50], and
for solving PDEs on surfaces [51]. For this example, we consider the curvature
driven flow problem that evolves the topology of a complex interface based on
its local curvature, and can be tailored for image processing tasks like image
smoothing and enhancement [49, 52]. Specifically, we look at the problem of a
wound spiral which relaxes under its own curvature [48]. The parametrization
of the initial spiral is given by $\theta=2\pi D\sqrt{s}$ with
$s=(k+a)/(np+a)$. The point locations of the zero contour can be computed as,

$\displaystyle x_{s}$
$\displaystyle=x_{c}+m\left(D/(1+D)\sqrt{s}\cos(\theta)\right)$ $\displaystyle
y_{s}$ $\displaystyle=y_{c}+m\left(D/(1+D)\sqrt{s}\sin(\theta)\right)$

The zero level set function can then be computed using the distance function
$d(\bm{x})=\|\bm{x}-\bm{x}_{p}\|$ as $\phi(\bm{x},t=0)=d(\bm{x})-w$ with $w$
being the width of the spirals. The distance function $d(\bm{x})$ measures the
shortest distance between two points $\bm{x}$ and $\bm{x}_{p}$. Given the
level set function $\phi(\bm{x})$, we can compute the 2D curvature $\kappa$
as,

$\kappa=\frac{\phi_{xx}\phi_{y}\phi_{y}-2\phi_{xy}\phi_{x}\phi_{y}+\phi_{yy}\phi_{x}\phi_{x}}{\left(\phi_{x}^{2}+\phi_{y}^{2}\right)^{3/2}}$
(24)

The evolution of the level set function driven by the local curvature $\kappa$
and is computed using the PDE,

$\frac{\partial\phi}{\partial t}=\kappa|\nabla\phi|$

The simulation was conducted on an $128\times 128$ resolution spatial grid
covering the domain $[-2,2]\times[-2,2]$ using the finite-difference method
with necessary reinitialization to preserve the sign distance property of the
level-set function $\phi$. The spiral parameters are chosen to be
$np=400,D=2.5,a=3,m=2$ [53]. Time integration is performed using a second-
order Runge-Kutta scheme with a step size of $dt=0.001$.

### Training set-up

The training data is generated from uniformly sampling in the domain
$[-1.5,1.5]\times[-1.5,1.5]$ with a resolution of $100\times 100$ in space and
$80$ time snapshots separated with time-step $dt=0.001$. This lead to
$N_{I}=10000,N_{B}=32000$ and $N_{\text{int}}=768000$ points in the interior
of the space and time domain. The network is trained for $2500$ epochs using
the Adam optimizer with each epoch consisting of $187$ mini-batch updates with
batch size $|\textrm{B}|=4096$. The neural network
$f_{\theta}:(x,y,t)\rightarrow(\phi)$ with $5$ layers and $64$ nodes is chosen
for this problem. We found that using the ELU activation function resulted in
the best performance. The total loss is given as,

$\mathcal{L}(\bm{u})=\lambda_{1}\mathcal{L}_{\text{boundary}}+\lambda_{2}\mathcal{L}_{\text{res}}+\lambda_{3}\mathcal{L}_{\text{init}}$
(25)

We choose the initial learning rate as $\eta=10^{-3}$ and decrease it after
$1500$ and $2000$ epochs by a factor of $10$.

Figure E.1: Curvature driven level-set problem A) Relative L2 error for for
different methods for the curvature driven level-set problem. B) The Power
spectrum of the PINN level-set function approximation after 2000 epochs
compared to the Power spectrum of the ground-truth (shown as solid black
line).  Figure E.2: Solution snapshot of curvature-driven level-set problem
Top row) Zero contour of the level-set function visualized at simulation time
$t=0.01s$ and after $2000$ epochs of training. Different columns correspond to
different weighing strategies A) uniform weighing B) pinn-MGDA C) inverse-
Dirichlet D) max/avg. Bottom row) The corresponding point-wise error of the
PINN approximation of the level-set function in comparison to the ground truth
solution.

##  Appendix F Network initialization and parameterization

Let $\mathcal{X}\in\mathbb{R}^{N\times d_{inp}}$ be the set of all input
coordinates to train the neural network, where $N$ denotes the total number of
data points and $d_{inp}$ the dimension of the input data, i.e. $d_{inp}=3$
corresponding to $(x,y,t)$. Moreover,
$\bm{X}_{\textrm{B}}\in\mathbb{R}^{|\textrm{B}|\times d_{inp}}$ will denote a
single batch, by randomly subsampling $\mathcal{X}$ into $N/|\textrm{B}|$
subsets. We then define an $n$-layer deep neural network with parameters
$\bm{\theta}$ as a composite function
$f_{\bm{\theta}}:\mathbb{R}^{|\textrm{B}|\times
d_{inp}}\rightarrow\mathbb{R}^{|\textrm{B}|\times d_{out}}$, where $d_{out}$
corresponds to the dimension of the output variable. The output is usually
computed as

$f_{\bm{\theta}}(\bm{X}_{\textrm{B}})=T_{n}(\sigma(T_{n-1}(\sigma(\dots\sigma(T_{1}(\bm{X}_{\textrm{B}})))))),$
(26)

with the affine linear maps
$T_{i}(\bm{z}_{i-1})=\bm{z}_{i-1}\bm{W}^{T}_{i}+\bm{b}_{i}$, $i=1,\dots,n$,
$\bm{z}_{0}=\bm{X}_{\textrm{B}}$, and output
$\bm{Y}_{\textrm{B}}=\bm{z}_{n-1}\bm{W}^{T}_{n}+\bm{b}_{n}$. Here,
$\bm{W}_{i}\in\mathbb{R}^{q_{i}\times q_{i-1}}$ describes the weights of the
incoming edges of layer $i$ with $q_{i}$ neurons,
$\bm{b}_{i}\in\mathbb{R}^{q_{i}}$ contains a scalar bias for each neuron and
$\sigma(\cdot)$ is a non-linear activation function. The network parameters
$\bm{W}_{i},\bm{b}_{i}\subset\bm{\theta}$ are initialized using normal Xavier
initialization [54] with $\bm{W}_{i}\sim\mathcal{G}(0,\Psi^{2})$, where
$\mathcal{G}$ is a Gaussian distribution with zero mean and standard deviation

$\Psi=g\cdot\sqrt{\frac{2}{q_{i}+q_{i-1}}},$ (27)

and $\bm{b}_{i}=\bm{0}$. The gain $g$ was chosen to be $g=1$ for all
activation functions except for $\tanh$ with a recommended gain of
$g=\frac{5}{3}$. Before feeding $\bm{X}_{\textrm{B}}$ into the first layer of
the neural network, we normalize it with respect to the training data
$\mathcal{X}$ as

$\hat{\bm{X}}_{\textrm{B}}=\frac{\bm{X}_{\textrm{B}}-\bm{\mu}}{\bm{\psi}},$
(28)

where
$\bm{\mu}=((\overline{\mathcal{X}_{1}},\dots,\overline{\mathcal{X}_{d_{inp}}}))$
and
$\bm{\psi}=((\sqrt{\textrm{Var}\\{\mathcal{X}_{1}\\}},\dots,\sqrt{\textrm{Var}\\{\mathcal{X}_{d_{inp}}\\}}))$
are the column-wise mean and standard deviations of the training data
$\mathcal{X}$.

### Weights computation

We initalize all weights to $\lambda_{k}=1,k=1,\dots,K$. For all experiments
we update the dynamic weights at the first batch of every fifth epoch. The
moving average is then performed with $\alpha=0.5$. The mini-batches
$\bm{X}_{\textrm{B}}$ are randomized after each epoch, such that for every
update of the weights, $\lambda_{k}$ will be computed with respect to a
different subset $\bm{X}_{\textrm{B}}$ of all training points $\mathcal{X}$.
Furthermore, batching is only performed over the interior points
$N_{\textrm{int}}$. As $N_{I},N_{B}\ll N_{\textrm{int}}$, the points for
initial and boundary conditions don’t require batching and remain fixed
throughout the training.