h-analysis and data-parallel physics-informed neural networks
We explore the data-parallel acceleration of physics-informed machine
learning (PIML) schemes, with a focus on physics-informed neural networks
(PINNs) for multiple graphics processing units (GPUs) architectures. In order
to develop scale-robust and high-throughput PIML models for sophisticated
applications which may require a large number of training points (e.g.,
involving complex and high-dimensional domains, non-linear operators or
multi-physics), we detail a novel protocol based on $h$-analysis and
data-parallel acceleration through the Horovod training framework. The protocol
is backed by new convergence bounds for the generalization error and the
train-test gap. We show that the acceleration is straightforward to implement,
does not compromise training, and proves to be highly efficient and
controllable, paving the way towards generic scale-robust PIML. Extensive
numerical experiments with increasing complexity illustrate its robustness and
consistency, offering a wide range of possibilities for real-world simulations.
Paul Escapil-Inchauspé, Gonzalo A. Ruz
# \(h\)-analysis and data-parallel physics-informed neural networks

###### Abstract

We explore the data-parallel acceleration of physics-informed machine learning (PIML) schemes, with a focus on physics-informed neural networks (PINNs) for multiple graphics processing units (GPUs) architectures. In order to develop scale-robust PIML models for sophisticated applications (e.g., involving complex and high-dimensional domains, non-linear operators or multi-physics), which may require a large number of training points, we detail a protocol based on the Horovod training framework. This protocol is backed by \(h\)-analysis, including a new convergence bound for the generalization error. We show that the acceleration is straightforward to implement, does not compromise training, and proves to be highly efficient, paving the way towards generic scale-robust PIML. Extensive numerical experiments with increasing complexity illustrate its robustness and consistency, offering a wide range of possibilities for real-world simulations.

## Introduction

Simulating physics throughout accurate surrogates is a hard task for engineers and computer scientists. Numerical methods such as finite element methods, finite difference methods and spectral methods can be used to approximate the solution of partial differential equations (PDEs) by representing them as a finite-dimensional function space, delivering an approximation to the desired solution or mapping [1].

Real-world applications often incorporate partial information of physics and observations, which can be noisy. This hints at using data-driven solutions throughout machine learning (ML) techniques. In particular, deep learning (DL) [2, 3] principles have been praised for their good performance, granted by the capability of deep neural networks (DNNs) to approximate high-dimensional and non-linear mappings, and offer great generalization with large datasets. Furthermore, the exponential growth of GPUs capabilities has made it possible to implement even larger DL models.

Recently, a novel paradigm called physics-informed machine learning (PIML) [4] was introduced to bridge the gap between data-driven [5] and physics-based [6] frameworks. PIML enhances the capability and generalization power of ML by adding prior information on physical laws to the scheme by restricting the output space (e.g., via additional constraints or a regularization term). This simple yet general approach was applied successfully to a wide range complex real-world applications, including structural mechanics [7, 8] and biological, biomedical and behavioral sciences [9].

In particular, physics-informed neural networks (PINNs) [10] consist in applying PIML by means of DNNs. They encode the physics in the loss function and rely on automatic differentiation (AD) [11]. PINNs have been used to solve inverse problems [12], stochastic PDEs [13, 14], complex applications such as the Boltzmann transport equation [15] and large-eddy simulations [16], and to perform uncertainty quantification [17, 18].

Concerning the challenges faced by the PINNs community, efficient training [19], proper hyper-parameters setting [20], and scaling PINNs [21] are of particular interest. Regarding the latter, two research areas are gaining attention.

First, it is important to understand how PINNs behave for an increasing number of training points \(N\) (or equivalently, a decreasing maximum distance between points \(h\)). Throughout this work, we refer to this study as \(h\)-analysis. In their pioneer works [22, 23], Mishra and Molinaro provided a bound for the generalization error with respect to to \(N\) for data-free and unique continuation problems, respectively. More precise bounds have been obtained using characterizations of the DNN [24].

Second, PINNs are typically trained over graphics processing units (GPUs), which have limited memory capabilities. To ensure models scale well with increasingly complex settings, two paradigms emerge: data-parallel and model-parallel acceleration. The former splits the training data over different workers, while the latter distributes the model weights. However, general DL backends do not readily support multiple GPU acceleration. To address this issue, Horovod [25] is a distributed framework specifically designed for DL, featuring a ring-allreduce algorithm [26] and implementations for TensorFlow, Keras and PyTorch.
As model size becomes prohibitive, domain decomposition-based approaches allow for distributing the computational domain. Examples of such approaches include conservative PINNs (cPINNs) [27], extended PINNs (XPINNs) [28, 29], and distributed PINNs (DPINNs) [30]. cPINNs and XPINNs were compared in [31]. These approaches are compatible with data-parallel acceleration within each subdomain. Additionally, a recent review concerning distributed PIML [21] is also available. Regarding existing data-parallel implementations, TensorFlow MirroredStrategy in TensorDiffEq [32] and NVIDIA Modulus [33], which supports Horovod acceleration, should be mentioned. However, to the authors knowledge, there is no systematic study of the background of data-parallel PINNs and their implementation.

In this work, we present a procedure to attain data-parallel efficient PINNs. It relies on \(h\)-analysis and consists of a Horovod-based acceleration. Concerning \(h\)-analysis, we observe PINNs exhibiting three phases of behavior as a function of the number of training points \(N\):

1. A pre-asymptotic regime, where the model does not learn the solution due to missing information;
2. A transition regime, where the error decreases with \(N\);
3. A permanent regime, where the error remains stable.

To illustrate this, Fig. 1 presents the relative \(L^{2}\) error distribution with respect to \(N_{f}\) (number of domain collocation points) for the forward "1D Laplace" case. The experiment was conducted over 8 independent runs with a learning rate of \(10^{-4}\) and 20000 iterations of ADAM [34] algorithm. The transition regime--where variability in the results is high and some models converge while others do not--is between \(N_{f}=32\) and \(N_{f}=256\). For more information on the experimental setting and the definition of precision \(\rho\), please refer to "1D Laplace".

Building on the empirical observations, we use the setting in [22, 23] to supply a rigorous theoretical background to \(h\)-analysis. One of the main contributions of this manuscript is the bound on the "Generalization error for generic PINNs", which allows for a simple analysis of the \(h\)-dependence.

Any practitioner strives to reach the permanent regime for their PIML scheme, and we provide the necessary details for an easy implementation of Horovod-based data acceleration for PINNs, with direct application to any PIML model. Fig. 2 further illustrates the scope of data-parallel PIML. Next, we apply the procedure to increasingly complex problems and demonstrate that Horovod acceleration is straightforward, using the pioneer PINNs code of Raissi as an example. Our main practical findings concerning data-parallel PINNs for up to 8 GPUs are the following :

* They do not require to modify their hyper-parameters;
* They show similar training convergence to the 1 GPU-case;
* They lead to high efficiency for both weak and strong scaling (e.g \(E_{\text{ff}}\geq 81\%\) for Navier-Stokes problem with 8 GPUs).

This work is organized as follows: In "Problem formulation", we introduce the PDEs under consideration, PINNs and convergence estimates for the generalization error. We then move to "Data-parallel PINNs" and present "Numerical experiments". Finally, we close this manuscript in "Conclusion".

Figure 1: Error v/s the number of domain collocation points \(N_{f}\) for the “1D Laplace” case. A pre-asymptotic regime (pink) is followed by a rapid transition regime (blue), and eventually leading to a permanent regime (green). This transition occurs over a few extra training points.
### Problem formulation

#### General notation

Throughout, vector and matrices are expressed using bold symbols. For a natural number \(k\), we set \(\mathbb{N}_{k}:=\{k,k+1,\cdots\}\). For \(p\in\mathbb{N}_{0}=\{0,1,\cdots,\}\), and an open set \(D\subseteq\mathbb{R}^{d}\) with \(d\in\mathbb{N}_{1}\), let \(L^{p}(D)\) be the standard class of functions with bounded \(L^{p}\)-norm over \(D\). Given \(s\in\mathbb{R}^{+}\), we refer to [1, Section 2] for the definitions of Sobolev function spaces \(H^{s}(D)\). Norms are denoted by \(\|\cdot\|\), with subscripts indicating the associated functional spaces. For a finite set \(\mathcal{T}\), we introduce notation \(|\mathcal{T}|:=\text{card}(\mathcal{T})\), closed subspaces are denoted by a \(\underset{\text{cl}}{\subset}\)-symbol and \(t^{2}=-1\).

### Abstract PDE

In this work, we consider a domain \(D\subset\mathbb{R}^{d}\), \(d\in\mathbb{N}_{1}\), with boundary \(\Gamma=\partial D\). For any \(T>0\), \(\mathbb{D}:=D\times[0,T]\), we solve a general non-linear PDE of the form:

\[\begin{cases}\mathcal{N}[u(\mathbf{x},t);\lambda]&=f(\mathbf{x},t)\quad( \mathbf{x},t)\in D_{f}=:D\times[0,T],\\ \mathcal{B}[u(\mathbf{x},t);\lambda]&=g(\mathbf{x},t),\quad(\mathbf{x},t)\in D _{g}=:\Gamma\times[0,T],\\ u(\mathbf{x},0)&=\hbar(\mathbf{x}),\quad\mathbf{x}\in D_{\hbar}=:D,\end{cases} \tag{1}\]

with \(\mathcal{N}\) a spatio-temporal differential operator, \(\mathcal{B}\) the boundary conditions (BCs) operator, \(\lambda\) the material parameters--the latter being unknown for inverse problems--and \(u(\mathbf{x},t)\in\mathbb{R}^{m}\) for any \(m\in\mathbb{N}_{1}\). Accordingly, for any function \(\tilde{u}\) defined over \(\mathbb{D}\), we introduce \(\Lambda=\{f,g,\hbar,u\}\) and define the residuals \(\xi_{v}\) for each \(v\in\Lambda\) and any observation function \(u_{\text{obs}}\):

\[\begin{cases}\xi_{f}(\mathbf{x},t;\lambda)&:=\mathcal{N}[\tilde{u}(\mathbf{x},t);\lambda]-f(\mathbf{x},t)\quad\text{in}\quad D_{f},\\ \xi_{g}(\mathbf{x},t;\lambda)&:=\mathcal{B}[\tilde{u}(\mathbf{x},t);\lambda]- g(\mathbf{x},t)\quad\text{in}\quad D_{g},\\ \xi_{h}(\mathbf{x},0)&:=\hat{u}(\mathbf{x},0)-\hbar(\mathbf{x})\quad\text{in }\quad D_{h},\\ \xi_{u}(\mathbf{x},t)&:=\hat{u}(\mathbf{x},t)-u_{\text{obs}}(\mathbf{x},t) \quad\text{in}\quad D_{u}:=\mathbb{D}.\end{cases} \tag{2}\]

#### PINNs

Following [20, 35], let \(\sigma\) be a smooth activation function. Given an input \((\mathbf{x},t)\in\mathbb{R}^{d+1}\), we define \(\mathcal{N}_{\mathcal{B}}\) as being a \(L\)-layer neural feed-forward neural network with \(W_{0}=d+1\), \(W_{L}=m\) and \(W_{l}\) neurons in the \(l\)-th layer for \(1\leq l\leq L-1\). For constant width DNNs, we set \(W=W_{1}=\cdots=W_{L-1}\). For \(1\leq l\leq L\), let us denote the weight matrix and bias vector in the \(l\)-th layer by \(\mathbf{W}^{l}\in\mathbb{R}^{d_{l}\times d_{l-1}}\) and \(\mathbf{b}^{l}\in\mathbb{R}^{d_{l}}\), respectively, resulting in:

\[\begin{array}{ll}\text{input layer:}&(\mathbf{x},t)\in\mathbb{R}^{d+1},\\ \text{hidden layers:}&\mathbf{z}^{l}(\mathbf{x})=\sigma(\mathbf{W}^{l}\mathbf{z} ^{l-1}(\mathbf{x})+\mathbf{b}^{l})\in\mathbb{R}^{d_{l}}\qquad\text{for}\quad 1 \leq l\leq L-1,\\ \text{output layer:}&\mathbf{z}^{L}(\mathbf{x})=\mathbf{W}^{l}\mathbf{z}^{L-1} (\mathbf{x})+\mathbf{b}^{L}\in\mathbb{R}^{m}.\end{array} \tag{3}\]

Figure 2: Scope of of data-parallel PIML.
This results in representation \(\mathbf{x}^{L}(\mathbf{x},t)\), with

\[\theta:=\left\{(\mathbf{W}^{1},\mathbf{b}^{1}),\cdots,(\mathbf{W}^{L},\mathbf{b} ^{L})\right\}, \tag{4}\]

the (trainable) parameters--or weights--in the network. We set \(\Theta=\mathbb{R}^{|\Theta|}\). Application of PINNs to Eq. (1) yields the approximate \(u_{\theta}(\mathbf{x},t)=\mathbf{x}^{L}(\mathbf{x},t)\).

We introduce the training dataset \(\mathcal{T}_{v}:=\{\tau^{i}_{v}\}_{i=1}^{N_{v}},\ \tau^{i}_{v}\in D_{v}\), \(N_{v}\in\mathbb{N}\) for \(i=1,\cdots,N_{v},v\in\Lambda\) and observations \(u_{\text{obs}}(\tau^{i}_{u})\), \(i=1,\cdots,N_{u}\). Furthermore, to each training point \(\tau^{i}_{v}\) we associate a quadrature weight \(w^{i}_{v}>0\). All throughout this manuscript, we set:

\[M:=N_{u},\quad\hat{N}:=N_{f}+N_{g}+N_{h}\quad\text{and}\quad N:=\hat{N}+M. \tag{5}\]

Note that \(M\) (resp. \(\hat{N}\)) represents the amount of information for the data-driven (resp. physics) part, by virtue of the PIML paradigm (refer to Fig. 2). The network weights \(\theta\) in Eq. (4) are trained (e.g., via ADAM optimizer [34]) by minimizing the weighted loss:

\[\mathcal{L}_{\theta}:=\sum_{v\in\Lambda}\omega_{v}\mathcal{L}^{v}_{\theta}, \quad\text{wherein}\quad\mathcal{L}^{v}_{\theta}:=\sum_{i=1}^{N_{v}}w^{i}_{v} |\xi_{v}(\tau^{i}_{v})|^{2}\quad\text{and}\quad\omega_{v}>0\quad\text{for} \quad v\in\Lambda. \tag{6}\]

We seek at obtaining:

\[\theta^{\star}:=\operatorname*{argmin}_{\theta\in\Theta}(\mathcal{L}_{\theta}). \tag{7}\]

The formulation for PINNs addresses the cases with no data (i.e. \(M=0\)) or physics (i.e. \(\hat{N}=0\)), thus exemplifying the PIML paradigm. Furthermore, it is able to handle time-independent operators with only minor changes; a schematic representation of a forward time-independent PINN is shown in Fig. 3.

Our setting assumes that the material parameters \(\lambda\) are known. If \(M>0\), one can solve the inverse problem by seeking:

\[(\theta^{\star}_{\text{inverse}},\lambda^{\star}_{\text{inverse}}):= \operatorname*{argmin}_{\theta\in\Theta,\lambda}\mathcal{L}_{\theta}[\lambda]. \tag{8}\]

Similarly, unique continuation problems [23], which assume incomplete information for \(f,g\) and \(\hbar\), are solved throughout PINNs without changes. Indeed, "2D Navier-Stokes" combines unique continuation problem and unknown parameters \(\lambda_{1},\lambda_{2}\).

### Automatic Differentiation

We aim at giving further details about back-propagation algorithms and their dual role in the context of PINNs:

1. Training the DNN by calculating \(\frac{\partial\mathcal{L}_{\theta}}{\partial\theta}\);

Figure 3: Schematic representation of a PINN. A DNN with \(L=3\) (i.e. \(L-1=2\) hidden layers) and \(W=5\) learns the mapping \(\mathbf{x}\mapsto u(\mathbf{x})\). The PDE is taken into account throughout the residual \(\mathcal{L}_{\theta}\), and the trainable weights are optimized, leading to optimal \(\theta^{\star}\).
2. Evaluating the partial derivatives in \(\mathcal{N}[u_{\theta}(\mathbf{x},t);\lambda]\) and \(\mathcal{B}[u_{\theta}(\mathbf{x},t);\lambda]\) so as to compute the loss \(\mathcal{L}_{\theta}\).

They consist in a forward pass to evaluate the output \(u_{\theta}\) (and \(\mathcal{L}_{\theta}\)), and a backward pass to assess the derivatives. To further elucidate back-propagation, we reproduce the informative diagram from [11] in Fig. 4. TensorFlow includes reverse mode AD by default. Its cost is bounded with \(|\Theta|\) for scalar output NNs (i.e. for\(m=1\)). The application of back-propagation (and reverse mode AD in particular) to any training point is independent of other information, such as neighboring points or the volume of training data. This allows for data-parallel PINNs. Before detailing its implementation, we justify the \(h\)-analysis through an abstract theoretical background.

### Convergence estimates

To understand better how PINNs scale with \(N\), we follow the method in [22, 23] under a simple setting, allowing to control the data and physics counterparts in PIML. Set \(s\geq 0\) and define spaces:

\[\hat{Y}\underset{\mathrm{cl}}{\subset}Y^{*}\underset{\mathrm{cl}}{ \subset}Y=L^{2}(\mathbb{D},\mathbb{R}^{m})\quad\text{and}\quad\hat{X} \underset{\mathrm{cl}}{\subset}X^{*}\underset{\mathrm{cl}}{\subset}X=H^{s}( \mathbb{D},\mathbb{R}^{m}) \tag{9}\]

We assume that Eq. (1) can be recast as:

\[\mathsf{A}u =b\quad\text{with}\quad\mathsf{A}:X^{*}\to Y^{*}\quad\text{and} \quad b\in Y^{*}, \tag{10}\] \[u =u_{\text{obs}}\quad\text{in}\quad X^{*}. \tag{11}\]

We suppose that Eq. (10) is well-posed and that for any \(u,v\in\hat{X}\), there holds that:

\[\|u-v\|_{Y}\leq C_{\text{pde}}(\|u\|_{\hat{X}},\|v\|_{\hat{X}})\left(\|\mathsf{ A}u-\mathsf{A}v\|_{Y}\right). \tag{12}\]

Eq. (10) is a stability estimate, allowing to control the total error by means of a bound on PINNs residual. Residuals in Eq. (2) are:

\[\begin{cases}\xi_{D}&:=\mathsf{A}u-b\quad\text{in}\quad Y^{*},\\ \xi_{u}&=\hat{u}-u_{\text{obs}}\quad\text{in}\quad X^{*}.\end{cases} \tag{13}\]

From the expression of residuals, we are interested in approximating integrals:

\[\overline{s}=\int_{\mathbb{D}}g(y)dy\quad\text{and}\quad\overline{l}=\int_{ \mathbb{D}}l(z)dz\quad\text{for}\quad g\in\hat{Y},l\in\hat{X}. \tag{14}\]

We assume that we are provided quadratures:

\[\overline{s}_{\hat{X}}=\sum_{i=1}^{\hat{X}}w_{i}g(\tau_{i}^{D})\quad\text{and }\quad\overline{l}_{M}=\sum_{i=1}^{M}q_{i}l(\tau_{i}^{u}) \tag{15}\]

Figure 4: Overview of back-propagation. A forward pass generates activations \(y_{i}\) and computes the error \(\mathcal{L}_{\theta}(y_{3},u)\). This is followed by a backward pass, through which the error adjoint is propagated to obtain the gradient with respect to weights \(\nabla\mathcal{L}_{\theta}\) where \(\theta=(w_{1},\cdots,w_{6})\). Additionally, the gradient \(\nabla_{x}\mathcal{L}_{\theta}\) can also be computed in the same backward pass.
for weights \(w_{i},q_{i}\) and quadrature points \(\tau_{i}^{D},\tau_{i}^{u}\in\mathbb{D}\) such that for \(\alpha,\beta>0\):

\[|\bar{g}-\bar{g}_{\beta}|\leq C_{\mathrm{quad},Y}\hat{N}^{-\alpha}\quad\text{ and}\quad|\bar{l}-\bar{l}_{M}|\leq C_{\mathrm{quad},X}M^{-\beta}. \tag{16}\]

For any \(\omega_{u}>0\), the loss is defined as follows:

\[\mathcal{L}_{\theta} =\sum_{i=1}^{\hat{N}}w_{i}|\xi_{D,\theta}(\tau_{i}^{D})|^{2}+ \omega_{u}\sum_{i=1}^{M}w_{i}^{u}|\xi_{u,\theta}(\tau_{i}^{u})|^{2}\approx\|\xi _{D,\theta}\|_{Y}^{2}+\omega_{u}\|\xi_{u,\theta}\|_{X}^{2}\] \[=\varepsilon_{T,D}^{2}+\omega_{u}\varepsilon_{T,u}^{2},\]

with \(\varepsilon_{T,D}\) and \(\varepsilon_{T,u}\) the training error for collocation points and observations respectively. Notice that application of Eq. (16) to \(\xi_{D,\theta}\) and \(\xi_{u,\theta}\) yields:

\[\|\xi_{D,\theta}\|_{Y}^{2}-\varepsilon_{T,D}^{2}\|\leq C_{\mathrm{quad},Y}\hat {N}^{-\alpha}\quad\text{and}\quad||\xi_{u,\theta}||_{X}^{2}-\varepsilon_{T,u} ^{2}|\leq C_{\mathrm{quad},X}M^{-\beta}. \tag{17}\]

We seek to quantify the _generalization error_:

\[\varepsilon_{G}=\varepsilon_{G}(\theta^{\star}):=\|u-u^{\star}\|_{X}\quad \text{with }u^{\star}:=u_{\theta^{\star}}\quad\text{and}\quad\theta^{\star}:=\arg\min_{ \theta}\mathcal{L}_{\theta}. \tag{18}\]

We detail a new result concerning the generalization error for PINNs.

**Theorem 1** (Generalization error for generic PINNs): _Under the presented setting, there holds that:_

\[\varepsilon_{G}\leq\frac{C_{\mathrm{pde}}}{1+\omega_{u}}\left( \varepsilon_{T,D}+C_{\mathrm{quad},Y}^{1/2}\hat{N}^{-\alpha/2}\right)+\frac{ \omega_{u}}{1+\omega_{u}}\left(\varepsilon_{T,u}+C_{\mathrm{quad},X}^{1/2}M^ {-\beta/2}+\hat{\mu}\right) \tag{19}\]

_with \(\hat{\mu}:=\|u-u_{\mathrm{obs}}\|_{X}\)._

**Proof 1**: _Consider the setting of Theorem 1. There holds that:_

\[(1+\omega_{u}) \varepsilon_{G} =\|u-u^{\star}\|_{X}+\omega_{u}\|u-u^{\star}\|_{X},\quad\text{by Eq. \eqref{eq:pde}}\] \[\leq C_{\mathrm{pde}}\|\hat{\Delta}u-\Delta u^{\star}\|_{Y}+ \omega_{u}\|u-u^{\star}\|_{X},\quad\text{by Eq. \eqref{eq:pde}}\] \[\leq C_{\mathrm{pde}}\|\xi_{D,\theta^{\star}}\|_{Y}+\omega_{u}\|u -u_{\mathrm{obs}}\|_{X}+\omega_{u}\|u_{\mathrm{obs}}-u^{\star}\|_{X},\quad \text{by Eq. \eqref{eq:pde} and triangular inequality}\] \[=C_{\mathrm{pde}}\|\xi_{D,\theta^{\star}}\|_{Y}+\omega_{u}\|\xi_{u,\theta^{\star}}\|_{Y}+\omega_{u}\hat{\mu},\quad\text{by definition}\] \[\leq C_{\mathrm{pde}}\varepsilon_{T,D}+\omega_{u}\varepsilon_{T,u} +C_{\mathrm{pde}}C_{\mathrm{quad},Y}^{1/2}\hat{N}^{-\alpha/2}+\omega_{u}C_{ \mathrm{quad},X}^{1/2}M^{-\beta/2}+\omega_{u}\hat{\mu},\quad\text{by Eq. \eqref{eq:pde}.}\]

The novelty of Theorem 1 is that it describes the generalization error for a simple case involving collocation points and observations. To make the result more intuitive, we rewrite Eq. (19), with \(\sim\) expressing the terms up to positive constants:

\[\varepsilon_{G}\sim(\hat{N}^{-\alpha/2}+M^{-\beta/2})+\varepsilon_{T,D}+ \varepsilon_{T,u}+\hat{\mu}. \tag{20}\]

The generalization error depends on the training errors (which are tractable during training), parameters \(\hat{N}\) and \(M\) and bias \(\hat{\mu}\).

To return to \(h\)-analysis, we now have a theoretical proof of the three regimes presented in introduction. Let us assume that \(\hat{\mu}=0\). For small values of \(\hat{N}\) or \(M\), the bound in Theorem 1 is too high to yield a meaningful estimate. Subsequently, the convergence is as \(\max(\hat{N}^{-\alpha/2},M^{-\beta/2})\), marking the transition regime. It is paramount for practitioners to reach the permanent regime when training PINNs, giving ground to data-parallel PINNs.

## Data-parallel PINNs

### Data-distribution and Horovod

In this section, we present the data-parallel distribution for PINNs. Let us set size\(\in\mathbb{N}_{1}\) and define ranks (or workers):

\[\text{rank}=0,\cdots,\text{size}-1,\]

each rank corresponding generally to a GPU. Data-parallel distribution requires the appropriate partitioning of the training points across ranks.
We introduce \(\hat{N}_{1},M_{1}\in\mathbb{N}_{1}\) collocation points and observations, respectively, for each rank (e.g., a GPU) yielding:

\[\mathcal{T}_{\nu}=\bigcup_{\text{rank}=0}^{\text{size}-1}\mathcal{T}_{\nu}^{ \text{rank}}\quad\text{for}\quad v\in\{D,u\},\]

with

\[\hat{N}=\texttt{size}\times\hat{N}_{1},\quad M=\texttt{size}\times M_{1}\quad \text{and}\quad\mathcal{T}=\mathcal{T}_{\nu}\cup\mathcal{T}_{u}\quad\text{with} \quad N=\hat{N}+M. \tag{21}\]

Data-parallel approach is as follows: We send the same synchronized copy of the DNN \(\mathcal{N}\mathcal{N}_{\theta}\) defined in Eq. (3) to each rank. Each rank evaluates the loss \(\mathcal{L}_{\theta}^{\text{rank}}\) and the gradient \(\nabla_{\theta}\mathcal{L}_{\theta}^{\text{rank}}\). The gradients are then averaged using an all-reduce operation, such as the ring all-reduce implemented in Horovod [26], which is known to be optimal with respect to the number of ranks. The process is illustrated in Figure 5 for \(\texttt{size}=4\). The ring-allreduce algorithm involves each of the size nodes communicating with two of its peers \(2\times(\texttt{size}-1)\) times [26].

It noteworthy that data generation for data-free PINNs (i.e. with \(M=0\)) requires no modification to existing codes, provided that each rank has a different seed for random or pseudo-random sampling. Horovod allows to apply data-parallel acceleration with minimal changes to existing code. Moreover, our approach and Horovod can easily be extended to multiple computing nodes. As pointed out in "Introduction", supports popular DL backends such as TensorFlow, PyTorch and Keras. In Listing 1, we demonstrate how to integrate data-parallel distribution using Horovod with a generic PINNs implementation in TensorFlow 1.x. The highlighted changes in pink show the steps for incorporating Horovod, which include: (i) initializing Horovod; (ii) pinning available GPUs to specific workers; (iii) wrapping the Horovod distributed optimizer and (iv) broadcasting initial variables to the master rank being rank = 0.

Figure 5: Data-parallel framework. Horovod supports ring-allreduce algorithm.
```
1#InitializeHorovod
2importhorovod.tensorflowashvd
3hvd.init()
4
5#PinGPUtobecusedtoprocesslocalrank(oneGPUperprocess)
6config=tf.ConfigProto()
7config.gpu_options.visible_device_list=str(hvd.local_rank())
8
9#BuildthePINN
10loss-...
11opt=tf.train.AdamOptimizer()
12#addMicroodDistributedOptimizer
13opt=hvd.DistributedOptimizer(opt)
14
15train=opt.minimize(loss)
16
17#Initializevariables
18init=tf.global_variables_initializer()
19self.ess.run(init)
20
21#Broadcastvariablesfromrank0totherworkers
22bcast=hvd.broadcast_global_variables(0)
23self.ess.run(bcast)
24
25#Trainthemodel
26whilen<-maxiter:
27sess.run(train)
28n+-1
```

Listing 1: Horovod for PINNs with TensorFlow 1.x. Data-parallel Horovod PINNs require minor changes to existing code.

### Weak and strong scaling

Two key concepts in data distribution paradigms are weak and strong scaling, which can be explained as follows: Weak scaling involves increasing the problem size proportionally with the number of processors, while strong scaling involves keeping the problem size fixed and increasing the number of processors. To reformulate:

* Weak scaling: Each worker has \((\hat{N}_{1},M_{1})\) training points, and we increase the number of workers size.
* Strong scaling: We set a fixed total number of \((\hat{N}_{1},M_{1})\) training points, and we split the data over increasing size workers.

We portray weak and strong scaling in Fig. 6 for a data-free PINN with \(\hat{N}_{1}=16\). Each box represents a GPU, with the number of collocation points as as color.

On the left of each scaling option, we present the unaccelerated case. We introduce the training time \(t_{\texttt{size}}\) for size workers. This allows to define the efficiency and speed-up as:

\[E_{\texttt{ff}}:=\frac{t_{1}}{t_{\texttt{size}}}\quad\text{and}\quad S_{\text{ up}}:=\texttt{size}\frac{t_{1}}{t_{\texttt{size}}}.\]

Figure 6: Weak and strong scaling for \(\hat{N}_{1}=16\) and size\(=8\).
## Numerical experiments

Throughout, we apply our proceeding to three cases of interest:

* "1D Laplace" equation (forward problem);
* "1D Schrodinger" equation (forward problem);
* "2D Navier-Stokes" equation (inverse problem).

For each case, we perform a \(h\)-analysis followed by Horovod data-parallel acceleration, which is applied to the domain training points (and observations for the Navier-Stokes case). Boundary loss terms are negligible due to the sufficient number of boundary data points.

### Methodology

We perform simulations in single float precision on a AMAX DL-E48A AMD Rome EPYC server with 8 Quadro RTX 8000 Nvidia GPUs--each one with a 48 GB memory. We use a Docker image of Horovod 0.26.1 with Python 3.6.9 and Tensorflow 2.6.2. All throughout, we use tensorflow.compat.vl as a backend without eager execution.

All the results are ready for use in HorovodPINNs GitHub repository and fully reproducible. We run experiments 8 times with seeds defined as:

\[\mathsf{seed}+1000\times\mathsf{rank},\]

in order to obtain rank-varying training points. For domain points, Latin Hypercube Sampling is performed with with pyDOE 0.3.8. Boundary points are defined over uniform grids.

We use Glorot uniform initialization [3, Chapter 8]. "Error" refers to the \(L^{2}\)-relative error taken over \(\mathcal{T}^{\mathrm{test}}\), and "Time" stands for the training time in seconds. For each case, the loss in Eq. (6) is with unit weights \(\omega_{v}=1\) and Monte-Carlo quadrature rule \(w_{v}^{i}=\frac{1}{N_{v}}\) for \(v\in\Lambda\). Also, we set \(\mathrm{vol}(\mathbb{D})\) the volume of domain \(\mathbb{D}\) and

\[\rho:=\frac{N_{f}{}^{1/(d+1)}}{\mathrm{vol}(\mathbb{D})^{1/(d+1)}}\]

We introduce \(t^{k}\) the time per to perform \(k\) iterations. The training points processed by second is as follows:

\[\mathrm{pointsec}:=\frac{kN_{f}}{t^{k}}.\]

For the sake of simplicity, we summarize the parameters and hyper-parameters for each case in Table 1.

## 1D Laplace

We first consider the 1D Laplace equation in \(D=[-1,7]\) as being:

\[-\Delta u=f\quad\text{in}\quad D\quad\text{with}\quad f=\pi^{2}\sin(\pi x) \quad\text{and}\quad u(-1)=u(7)=0. \tag{22}\]

Acknowledge that \(u(x)=\sin(\pi x)\). We solve the problem for:

\[N_{f}=2^{i}\quad i=3,\cdots,14\quad\text{and}\quad N_{f}=80,100,105,110,115,120.\]

We set \(N_{g}=2\) and \(N_{u}=0\). Points in \(\mathcal{T}_{D}\) are generated randomly over \(D\), and \(\mathcal{T}_{b}=\{-1,7\}\). The residual in Eq. (2):

\[\xi^{f}=-\Delta u-f\]

yields the loss:

\[\mathcal{L}_{\theta}=\mathcal{L}_{\theta}^{f}+\mathcal{L}_{\theta}^{b}\quad \text{with}\quad\mathcal{L}_{\theta}^{f}:=\frac{1}{N_{f}}\sum_{\mathcal{J}}| \xi_{f}|^{2}\quad\text{and}\quad\mathcal{L}_{\theta}^{b}:=\frac{1}{2}\left(|u (-1)|^{2}+|u(7)|^{2}\right).\]

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline Case & 
\begin{tabular}{c} learning rate \\ \(l_{r}\) \\ \end{tabular} & width & depth & & & & & \\ \hline \hline
1D Laplace & \(10^{-4}\) & 50 & 4 & 20000 & 7801 & \(N_{f}\) & \(\tanh\) \\ \hline
1D Schrödinger & \(10^{-4}\) & 50 & 4 & 30000 & 30802 & \(N_{f}\) & \(\tanh\) \\ \hline
2D Navier-Stokes & \(10^{-4}\) & 20 & 8 & 30000 & 3604 & \(N_{f}\) & \(\tanh\) \\ \hline \end{tabular}
\end{table}
Table 1: Overview of the parameters and hyper-parameters for each case.
\(h\)-analysis

We perform the \(h\)-analysis for the error as portrayed before in Figure 1. The asymptotic regime occurs between \(N_{f}=32\) and \(N_{f}=256\), with a precision of \(\rho=8\) in accordance with general results for \(h\)-analysis of traditional solvers. The permanent regime shows a slight improvement in accuracy, with Error dropping from \(10^{-3}\) to \(10^{-4}\) approximately. To complete the \(h\)-analysis, Figure 7 shows the convergence results of ADAM optimizer for all the values of \(N_{f}\). This plot reveals that each regime exhibits similar patterns. The high variability in convergence during the transition regime is particularly interesting, with some runs converging and others not. In the permanent regime, the convergence shows almost identical and stable patterns irrespective of \(N_{f}\).

#### Data-parallel implementation

We set \(N_{f,1}\equiv N_{1}=32\) and compare both weak and strong scaling to original implementation, referred to as "no scaling".

We provide a detailed description of Fig. 8, as it will serve as the basis for future cases:

* Left-hand side: Error for \(\texttt{size}^{*}\in\{1,2,4,8\}\) corresponding to \(N_{f}\in\{16,32,64,128\}\).
* Middle: Error for weak scaling with \(N_{1}=16\) and for \(\texttt{size}\in\{1,2,4,8\}\);
* Right-hand side: Error for strong scaling with \(N_{1}=128\) and for \(\texttt{size}\in\{1,2,4,8\}\)

Figure 8: Error v/s \(\texttt{size}^{*}\) for the different scaling options.

Figure 7: 1D Laplace: Convergence for ADAM v/s \(N_{f}\) for \(l_{r}=10^{-4}\) and 20000 iterations.
To reduce ambiguity, we use the \(*\)-superscript for no scaling, as size\({}^{*}\) is performed over 1 rank. The color for each violin box in the figure corresponds to the number of domain training points used for each GPU.

Fig. 8 demonstrates that both weak and strong scaling yield similar convergence results to their unaccelerated counterpart. This result is one of the main findings in this work: PINNs scale properly with respect to accuracy, validating the intuition behind \(h\)-analysis and justifying the data-parallel approach. This allows one to move from pre-asymptotic to permanent regime by using weak scaling, or leverage the cost of a permanent regime application by dispatching the training points over different workers. Furthermore, the hyper-parameters, including the learning rate, remained unchanged.

Next, we summarize the data-parallelization results in Table 2 with respect to size.

In the first column, we present the time required to run 500 iterations for ADAM, referred to as \(t^{500}\). This value is averaged over one run of 30000 iterations with a heat-up of 1500 iteration (i.e. we discard the values corresponding to iteration \(0,500\) and \(1000\)). We present the resulting mean value \(\pm\) standard deviation for the resulting vector. The second column, displays the efficiency of the run, evaluated with respect to \(t^{500}\).

The table reveals that data-based acceleration results in supplementary training times as anticipated. Weak scaling efficiency varies between \(70.8\%\) for \(\texttt{size}=1\) to \(60.1\%\) for \(\texttt{size}=8\), resulting in a speed-up of \(4.80\) when using 8 GPUs. Similarly, strong scaling shows similar behavior. Furthermore, it can be observed that \(\texttt{size}=1\) yields almost equal \(t_{500}\) for \(N_{f}=16\) (\(5.65s\)) and \(N_{f}=128\) (\(5.76\%\)).

To conclude, the Laplace case is not reaching its full potential with Horovod due to the small batch size \(N_{1}\). However, increasing the value of \(N_{1}\) in future cases will lead to a noticeable increase in efficiency.

## 1D Schrodinger

We solve the non-linear Schrodinger equation along with periodic BCs (refer to [10, Section 3.1.1]) given over \(D\times T\) with \(D:=(-5,5)\) and \(T:=\pi/2\)

\[u_{t}+0.5u_{xx}+|u|^{2}u =0\quad\text{in}\quad D\times T,\] \[u(-5,t) =u(5,t),\] \[\partial_{x}u(-5,t) =\partial_{x}u(5,t),\] \[u(x,0) =2\,\text{sech}(x),\]

where \(u(x,t)=u^{0}(x,t)+u^{1}(x,t)\). We apply PINNs to the system with \(m=2\) as \((u_{\theta}^{0},v_{\theta}^{1})\in\mathbb{R}^{m}=\mathbb{R}^{2}\). We set \(N_{g}=N_{h}=200\).

\(h\)_-analysis_

To begin with, we perform the \(h\)-analysis for the parameters in Fig. 9. Again, the transition regime for the total execution time distribution begins at a density of \(\rho=6\) and spans between \(N_{f}=350\) and \(N_{f}=4000\). At higher magnitudes of \(N\), the error remains approximately the same. To illustrate this more complex case, we present the total execution time distribution in Fig. 10. We note that the training times remain stable for \(N_{f}\leq 4000\). This observation is of great importance and should be emphasized, as it adds additional parameters to the analysis. Our work primarily focuses on error in \(h\)-analysis, however, execution time and memory requirements are also important considerations. We see that in this case, weak scaling is not necessary (the optimal option is to use \(\texttt{size}=1\) and \(N_{f}=4000\)). Alternatively, strong scaling can be done with \(N_{f,1}=4000\).

To gain further insight into the variability of the transition regime, we focus on \(N_{f}=1000\). We compare the solution for \(\texttt{seed}=1234\) and \(\texttt{seed}=1236\). The upper figures depict \(|u(t,x)|\) predicted by the PINN. The lower figures show a comparison between the exact and predicted solutions are plotted for \(t\in\{0.59,0.79,-.98\}\). It is evident that the solution for \(\texttt{seed}=1234\) closely resembles the exact solution, whereas the solution for \(\texttt{seed}=1236\) fails to accurately represent the solution near \(x=0\), thereby illustrating the importance of achieving a permanent regime.

\begin{table}
\begin{tabular}{|c||c|c||c|c|} \hline  & \multicolumn{2}{c||}{\(t^{500}\)} & \multicolumn{2}{c|}{\(E_{\text{ff}}\)} \\ \hline size & Weak scaling & Strong scaling & Weak scaling & Strong scaling \\ \hline \hline
1 & \(5.65\pm 0.15\) & \(5.76\pm 0.25\) & – & – \\ \hline
2 & \(7.98\pm 0.29\) & \(7.94\pm 0.27\) & \(70.8\%\) & \(72.5\%\) \\ \hline
4 & \(8.72\pm 0.14\) & \(8.72\pm 0.11\) & \(64.8\%\) & \(66.0\%\) \\ \hline
8 & \(9.40\pm 0.18\) & \(9.40\pm 0.18\) & \(60.1\%\) & \(61.2\%\) \\ \hline \end{tabular}
\end{table}
Table 2: 1D Laplace: \(t^{500}\) in seconds and efficiency \(E_{\text{ff}}\) for the weak and strong scaling.
#### Data-parallel implementation

We compare the error for unaccelerated and data-parallel implementations of simulations for \(N_{f,1}\equiv N_{1}=500\) in Fig. 12, analogous to the analysis in Fig. 8. Again, the error is stable with size. Both no scaling and weak scaling are similar. Strong scaling is unaffected by size. We plot the training time in Fig. 13. We observe that both weak and strong scaling increase linearly and slightly with size. Both scaling show similar behaviors. Fig. 14 demonstrates the efficiency with respect to size, with white bars representing the ideal scaling. Efficiency \(E_{\text{ff}}\) shows a gradual decrease with size, with results surpassing those of the previous section. The efficiency for \(\texttt{size}=8\) reaches 77% and 76% respectively for weak and strong scaling, representing a speed-up of 6.2 and 6.1.

Figure 11: Schrödinger: Solution for \(N_{f}=1000\) for two differents seeds.

Figure 10: Schrödinger: Time (in seconds) v/s \(N_{f}\).

Figure 9: Schrödinger: Error v/s \(N_{f}\).
### Inverse problem: Navier-Stokes equation

We consider the Navier-Stokes problem with \(D:=[-1,8]\times[-2,2]\), \(T=20\) and unknown parameters \(\lambda_{1}\lambda_{2}\in\mathbb{R}\). The resulting divergence-free Navier Stokes is expressed as follows:

\[\begin{cases}u_{t}+\lambda_{1}(uu_{x}+vu_{y})&=-p_{x}+\lambda_{2}(u_{xx}+u_{yy}),\\ v_{t}+\lambda_{1}(uv_{x}+vv_{y})&=-p_{y}+\lambda_{2}(v_{xx}+v_{yy}),\\ u_{x}+u_{t}&=0,\end{cases} \tag{23}\]

Figure 14: Schrödinger: Efficiency v/s size.

Figure 12: Schrödinger: Error v/s size\({}^{*}\) for the different scaling options.

Figure 13: Schrödinger: Time (in seconds) v/s size\({}^{*}\) for the different scaling options.
wherein \(u(x,y,t)\) and \(v(x,y,t)\) are the \(x\) and \(y\) components of the velocity field, and \(p(x,y,t)\) the pressure. We assume that there exists \(\varphi(x,y,t)\) such that:

\[u=\varphi_{x},\quad\text{and}\quad v=-\varphi_{x}. \tag{24}\]

Under Eq. (24), last row in Eq. (23) is satisfied. The latter leads to the definition of residuals:

\[\xi_{1}:=u_{t}+\lambda_{1}(uu_{x}+vu_{y})+p_{x}-\lambda_{2}(u_{xx} +u_{yy})\] \[\xi_{2}:=v_{t}+\lambda_{1}(uv_{x}+vv_{y})+p_{y}-\lambda_{2}(v_{xx} +v_{yy}),\]

hence the loss:

\[\mathcal{L}_{\theta}:=\mathcal{L}_{D,\theta}+\mathcal{L}_{u,\theta}\]

with

\[\mathcal{L}_{D,\theta}=\frac{1}{N}\sum_{i=1}^{N}\omega_{i}\left(|u(\mathbf{x} ^{i})-u^{i}|^{2}+|v(\mathbf{x}^{i})-v^{i}|^{2}\right)\quad\text{and}\quad \mathcal{L}_{u,\theta}=\frac{1}{N}\sum_{i=1}^{N}\omega_{i}\left(\xi_{1}^{2}( \mathbf{x}^{i})+\xi_{2}^{2}(\mathbf{x}^{i})\right).\]

Throughout this case, we have \(N_{f}=M\), and we plot the results with respect to \(N_{f}\). Acknowledge that \(N=2N_{f}\), and that both \(\lambda_{1},\lambda_{2}\) and the BCs are unkwown.

#### \(h\)-analysis

We conduct the \(h\)-analysis and show the error in Fig. 15, showing no differences with previous cases.

Surprisingly, the permanent regime is reached only for \(N_{f}=1000\), despite the problem being a 3-dimensional, non-linear, and inverse one. This corresponds to low values of \(\rho\), indicating that PINNs seem to prevent the curse of dimensionality. In fact, the it was achieved with only 1.2 points per unit per dimension. The total training time is presented in Fig. 16, where it can be seen to remain stable up to \(N_{f}=5000\), and then it increases linearly with \(N_{f}\).

Figure 16: Navier-Stokes: Time (in seconds) v/s \(N_{f}\).

Figure 15: Navier-Stokes: Error v/s \(N_{f}\).
#### Data-parallel implementation

We run the data-parallel simulations, setting \(N_{f,1}\equiv N_{1}=M_{1}=500\). As shown in Fig. 17, the simulations exhibit stable accuracy with size. The execution time increases moderately with size, as illustrated in Fig. 18. The training time decreases with \(N\) for the no scaling case.

However, this behavior is temporary (refer to Fig. 15 before). We conclude our analysis by plotting the efficiency with respect to size in Fig. 19. Efficiency lowers with increasing size, but shows the best results so far, with efficiency overpassing 80%. This encouraging result sets the stage for further exploration of more intricate applications.

Figure 19: Navier-Stokes: Efficiency v/s size.

Figure 17: Navier-Stokes: Error v/s size\({}^{*}\) for the different scaling options.

Figure 18: Navier-Stokes: Time (in seconds) v/s size\({}^{*}\) for the different scaling options.
## Conclusion

In this work, we proposed a novel data-parallelization approach for PIML with a focus on PINNs. We provided a thorough \(h\)-analysis and associated theoretical results to support our approach, as well as implementation considerations to facilitate implementation with Horovod data acceleration. Additionally, we ran reproducible numerical experiments to demonstrate the scalability of our approach. Further work include the implementation of Horovod acceleration to DeepXDE[35] library, coupling of localized PINNs with domain decomposition methods, and application on larger GPU servers (e.g., with more than 100 GPUs).

## Data availability

The code required to reproduce these findings are available to download from [https://github.com/pescap/HorovodPINNs](https://github.com/pescap/HorovodPINNs).
