Learning from Integral Losses in Physics Informed Neural Networks
This work proposes a solution for the problem of training physics-informed
networks under partial integro-differential equations. These equations require
an infinite or a large number of neural evaluations to construct a single
residual for training. As a result, accurate evaluation may be impractical, and
we show that naive approximations at replacing these integrals with unbiased
estimates lead to biased loss functions and solutions. To overcome this bias,
we investigate three types of potential solutions: the deterministic sampling
approaches, the double-sampling trick, and the delayed target method. We
consider three classes of PDEs for benchmarking; one defining Poisson problems
with singular charges and weak solutions of up to 10 dimensions, another
involving weak solutions on electro-magnetic fields and a Maxwell equation, and
a third one defining a Smoluchowski coagulation problem. Our numerical results
confirm the existence of the aforementioned bias in practice and also show that
our proposed delayed target approach can lead to accurate solutions with
comparable quality to ones estimated with a large sample size integral. Our
implementation is open-source and available at
https://github.com/ehsansaleh/btspinn.
Ehsan Saleh, Saba Ghaffari, Timothy Bretl, Luke Olson, Matthew West
# Learning from Integral Losses in Physics Informed Neural Networks

###### Abstract

This work proposes a solution for the problem of training physics informed networks under partial integro-differential equations. These equations require infinite or a large number of neural evaluations to construct a single residual for training. As a result, accurate evaluation may be impractical, and we show that naive approximations at replacing these integrals with unbiased estimates lead to biased loss functions and solutions. To overcome this bias, we investigate three types of solutions: the deterministic sampling approach, the double-sampling trick, and the delayed target method. We consider three classes of PDEs for benchmarking; one defining a Poisson problem with singular charges and weak solutions, another involving weak solutions on electro-magnetic fields and a Maxwell equation, and a third one defining a Smoluchowski coagulation problem. Our numerical results confirm the existence of the aforementioned bias in practice, and also show that our proposed delayed target approach can lead to accurate solutions with comparable quality to ones estimated with a large number of samples. Our implementation is open-source and available at [https://github.com/ehsansaleh/btspin](https://github.com/ehsansaleh/btspin).

## 1 Introduction

Physics Informed Neural Networks (PINNs) [23] can be described as solvers of a particular Partial Differential Equation (PDE). Typically, these problems consist of three defining elements. A sampling procedure selects a number of points for learning. Automatic differentiation is then used to evaluate the PDE at these points and define a residual. Finally, a loss function, such as the Mean Squared Error (MSE), is applied to these residuals, and the network learns the true solution by minimizing this loss through back-propagation and stochastic approximation. These elements form the basis of many methods capable of learning high-dimensional parameters. A wealth of existing work demonstrate the utility of this approach to solving a wide array of applications and PDE forms [16; 25; 15].

One particular problem in this area, is the prevalent assumption around our ability to accurately evaluate the PDE residuals for learning. In particular, many partial integro-differential forms may include an integral or a large summation within them. Weak solutions using the divergence or the curl theorems, or the Smoluchowski coagulation equation [28] are two examples of such forms. In such instances, an accurate evaluation of the PDE elements, even at a single point, can become impractical. Naive approximations, such as replacing integrals with unbiased estimates, can result in biased solutions, as we will show later. This work is dedicated to the problem of learning PINNs with loss functions containing a parametrized integral or large summation.

To tackle such challenging problems, this paper investigates three learning approaches: the deterministic sampling approach, the double-sampling trick, and the delayed target method. The deterministic
sampling approach takes the integration stochasticity away. However, it optimizes a different objective and deterministic sampling schemes may not scale well to higher sampling dimensions. The double-sampling trick enjoys strong theoretical guarantees, but relies on a particular point sampling process and a squared loss form and may not be suitable for offline learning. Bootstrapping neural networks with delayed targets has already been studied in other contexts, such as the method of learning from Temporal Differences (TD-learning) in approximate dynamic programming [27]. Time and time again, TD-learning has proven preferable over the ordinary MSE loss minimization (known as the Bellman residual minimization) [24; 5; 31; 4]. Another example is the recent trend of semi-supervised learning, where teacher-student frameworks result in accuracy improvements of classification models by pseudo-labelling unlabelled examples for training [8; 22; 1; 14].

The main contributions of this work are: (1) we formulate the integral learning problem under a general framework and show the biased nature of standard approximated loss functions; (2) we present three techniques to solve such problems, namely the deterministic sampling approach, the double-sampling trick, and the delayed target method; (3) we detail an effective way of implementation for the delayed target method compared to a naive one; and (4) we compare the efficacy of the proposed solutions using numerical examples on a Poisson problem with singular charges, a Maxwell problem with magnetic fields, and a Smoluchowski coagulation problem.

## 2 Problem Formulation

Consider a typical partial integro-differential equation

\[f_{\theta}(x):=\mathbb{E}_{P(x^{\prime}|x)}[g_{\theta}(x^{\prime})]+y. \tag{1}\]

The \(f_{\theta}(x)\) and \(g_{\theta}(x^{\prime})\) are parametrized, and \(y\) includes all the non-parametrized terms in the PDE. The right side of the equation serves as the target value for \(f_{\theta}(x)\). For notation simplicity, we assume a fixed value for \(x\) in the remainder of the manuscript. However, all analyses are applicable to randomized \(x\) without any loss of generality. Equation (1) is a general, yet concise, form for expressing partial integro-differential equations. To motivate, we will express two examples.

Example 1The Poisson problem is to solve the system \(\nabla^{2}U=\rho\) for \(U\) given a charge function \(\rho\). This is equivalent to finding a solution for the system

\[E=\nabla U,\qquad\nabla\cdot E=\rho. \tag{2}\]

A weak solution can be obtained through enforcing the divergence theorem over many volumes:

\[\int_{\partial\Omega}E\cdot\hat{n}\quad\mathrm{d}S=\iint_{\Omega}\nabla\cdot E \quad\mathrm{d}V, \tag{3}\]

where \(\hat{n}\) is the normal vector perpendicular to \(\mathrm{d}S\). The weak solutions can be preferable over the strong ones when dealing with singular or sparse \(\rho\) charges.

Figure 1: Training with the MSE loss under different sample sizes per-surface (\(N\)). The heatmaps show the analytical solution (left), the low-variance training with \(N=100\) (middle), and the high-variance training with \(N=1\) (right). The smaller the \(N\), the more biased the training objective becomes towards finding smoother solutions. The right panel shows the training curves; the training loss and the integration variance represent \(\hat{\mathcal{L}}_{\theta}(x)\) and \(\mathbb{V}_{P(x^{\prime}|x)}[g_{\theta}(x^{\prime})]\) in Equation (15), respectively. For \(N=1\), the training loss seems to be floored at the same value as the integration variance (i.e., approximately \(0.3\)). However, with \(N=100\), the model produces better solutions, lower training losses, and higher integration variances.
To solve this system, we parametrize \(E(x)\) as the gradient of a neural network predicting the \(U\) potentials. To convert this into the form of Equation (1), we replace the left integral in Equation (3) with an arbitrarily large Riemann sum and write

\[\int_{\partial\Omega}E\cdot\hat{n}\quad\mathrm{d}S=\frac{A}{M}\sum_{i=1}^{M}E_{ \theta}(x_{i})\cdot\hat{n}_{i}, \tag{4}\]

where \(A=\int_{\partial\Omega}1\;\mathrm{d}S\) is the surface area and the \(x_{i}\) samples are uniform on the surface. To convert this system into the form of Equation (1), we then define

\[x:=x_{1},\qquad f_{\theta}(x):=\frac{A}{M}E_{\theta}(x)\cdot\hat{ n}_{1},\qquad x^{\prime}\sim\text{Unif}(\{x_{i}\}_{i=2}^{M}),\] \[g_{\theta}(x_{i}):=-\frac{A(M-1)}{M}E_{\theta}(x_{i})\cdot\hat{ n}_{i},\qquad y:=\iint_{\Omega}\rho\mathrm{d}V. \tag{5}\]

Example 2In static electromagnetic conditions, one of the Maxwell Equations, the Ampere circuital law, is to solve the system \(\nabla\times A=B\), \(\nabla\times B=J\) for \(A\) given the current density \(J\) in the 3D space (we assumed a unit physical coefficient for simplicity). Here, \(B\) represents the magnetic field and \(A\) denotes the magnetic potential vector. A weak solution for this system can be obtained through enforcing the Stokes theorem over many volumes:

\[\int_{\partial\Omega}\nabla\times A\cdot\mathrm{d}\mathbf{l}=\iint_{\Omega}J \cdot\mathrm{d}S, \tag{6}\]

where \(\mathrm{d}S\) is an infinitesimal surface normal vector. Just like the Poisson problem, the weak solutions can be preferable over the strong ones when dealing with singular inputs, and this equation can be converted into the form of Equation (1) similarly.

Example 3The Smoluchowski coagulation equation simulates the evolution of particles into larger ones, and is described as

\[\frac{\partial\rho(x,t)}{\partial t}= \int_{0}^{x}K(x-x^{\prime},x^{\prime})\rho(x-x^{\prime},t)\rho(x ^{\prime},t)\mathrm{d}x^{\prime}-\int_{0}^{\infty}K(x,x^{\prime})\rho(x,t) \rho(x^{\prime},t)\mathrm{d}x^{\prime}, \tag{7}\]

where \(K(x,x^{\prime})\) is the coagulation kernel between two particles of size \(x\) and \(x^{\prime}\). The particle sizes \(x\) and \(x^{\prime}\) can be generalized into vectors, inducing a higher-dimensional PDE to solve. To solve this problem, we parametrize \(\rho(x,t)\) as the output of a neural network with parameters \(\theta\) and write

\[f_{\theta}(x):=\frac{\partial\rho_{\theta}(x,t)}{\partial t}, \qquad g_{\theta}^{(1)}(x^{\prime}):=A_{1}K(x-x^{\prime},x^{\prime})\rho_{ \theta}(x-x^{\prime},t)\rho_{\theta}(x^{\prime},t),\] \[g_{\theta}^{(2)}(x^{\prime}):=A_{2}K(x,x^{\prime})\rho_{\theta}(x,t)\rho_{\theta}(x^{\prime},t). \tag{8}\]
The \(x^{\prime}\) values in both \(g_{\theta}^{(1)}\) and \(g_{\theta}^{(2)}\) are sampled from their respective uniform distributions, and \(A_{1}\) and \(A_{2}\) are used to normalize the uniform integrals into expectations. Finally, \(y=0\) and we can define \(g_{\theta}(x^{\prime})\) in a way such that

\[\mathbb{E}_{x^{\prime}}[g_{\theta}(x^{\prime})]:=\mathbb{E}_{x^{\prime}}[g_{ \theta}^{(1)}(x^{\prime})]+\mathbb{E}_{x^{\prime}}[g_{\theta}^{(2)}(x^{\prime })]. \tag{9}\]

The standard way to solve systems such as Examples 1, 2, and 3 with PINNs, is to minimize the following mean squared error (MSE) loss [23; 9]:

\[\mathcal{L}_{\theta}(x):=\big{(}f_{\theta}(x)-\mathbb{E}_{P(x^{\prime}|x)}[g_ {\theta}(x^{\prime})]-y\big{)}^{2}. \tag{10}\]

Since computing exact integrals may be impractical, one may contemplate replacing the expectation in Equation (10) with an unbiased estimate, as implemented in NVIDIA's Modulus package [21]. This prompts the following approximate objective:

\[\hat{\mathcal{L}}_{\theta}(x):=\mathbb{E}_{\{x^{\prime}_{i}\}_{i=1}^{N}}\bigg{[} \big{(}f_{\theta}(x)-\frac{1}{N}\sum_{i=1}^{N}g_{\theta}(x^{\prime}_{i})-y \big{)}^{2}\bigg{]}. \tag{11}\]

We therefore analyze the approximation error:

\[\hat{\mathcal{L}}_{\theta}(x)= \mathbb{E}_{\{x^{\prime}_{i}\}}\bigg{[}\bigg{(}\big{(}f_{\theta}( x)-\mathbb{E}_{x^{\prime\prime}}[g_{\theta}(x^{\prime\prime})]-y\big{)}+ \big{(}\mathbb{E}_{x^{\prime\prime}}[g_{\theta}(x^{\prime\prime})]-\frac{1}{N }\sum_{i=1}^{N}g_{\theta}(x^{\prime}_{i})\big{)}\bigg{)}^{2}\bigg{]}. \tag{12}\]

By decomposing the squared sum, we get

\[\hat{\mathcal{L}}_{\theta}(x)= \mathbb{E}_{\{x^{\prime}_{i}\}}\big{[}\big{(}\big{(}f_{\theta}( x)-\mathbb{E}_{x^{\prime\prime}}[g_{\theta}(x^{\prime\prime})]-y\big{)}^{2} \big{]}+\mathbb{E}_{\{x^{\prime}_{i}\}}\big{[}(\mathbb{E}_{x^{\prime\prime}}[ g_{\theta}(x^{\prime\prime})]-\frac{1}{N}\sum_{i=1}^{N}g_{\theta}(x^{\prime}_{i}))^{2} \big{]}\] \[-2\mathbb{E}_{\{x^{\prime}_{i}\}}\big{[}\big{(}f_{\theta}(x)- \mathbb{E}_{x^{\prime\prime}}[g_{\theta}(x^{\prime\prime})]-y\big{)}\big{(} \mathbb{E}_{x^{\prime\prime}}[g_{\theta}(x^{\prime\prime})]-\frac{1}{N}\sum_{i =1}^{N}g_{\theta}(x^{\prime}_{i})\big{)}\big{]}. \tag{13}\]

Since \(\mathbb{E}_{x^{\prime\prime}}[g_{\theta}(x^{\prime\prime})]=\mathbb{E}_{\{x^{ \prime}_{i}\}}[\frac{1}{N}\sum_{i=1}^{N}g_{\theta}(x^{\prime}_{i})]\), the last term in Equation (13) is zero, and we have

\[\hat{\mathcal{L}}_{\theta}(x)=\mathcal{L}_{\theta}(x)+\mathbb{V}_{P(\{x^{ \prime}_{i}\}|x)}[\frac{1}{N}\sum_{i=1}^{N}g_{\theta}(x^{\prime}_{i})]. \tag{14}\]

Figure 2: Training the same problem as in Figure 1 with delayed targets and \(N=1\). _The top left panel_ shows a diverged training with \(M=100\) in Equation (24). _The lower left panel_ corresponds to \(M=10\), which has a converging training curve even though producing an overly smooth solution. In _the lower right panel_, we set \(\lambda=1\) which allowed setting the simulated \(M\) as \(1000\) while maintaining a stable training loss. In each panel, the left and right heatmaps show the main and the target model predictions, respectively, and the right plots show the training curves. The green curves show the training loss for the delayed target method, and the standard training curves with \(N=1\) and \(100\) are also shown using dotted red and blue lines for comparison, respectively. _The top right panel_ shows an example of deterministic vs. i.i.d. sampling of the surface points \(\{x^{\prime}_{i}\}_{i=1}^{N}\) in the Poisson problem. For each sampled sphere, the surface points and their normal vectors are shown with \(N=100\) samples. With deterministic sampling, the points are evenly spaced to cover the sampling domain.
If the \(x^{\prime}_{i}\) are sampled in an i.i.d. manner, Equation (14) simplifies further to

\[\hat{\mathcal{L}}_{\theta}(x)=\mathcal{L}_{\theta}(x)+\frac{1}{N}\mathbb{V}_{P(x ^{\prime}|x)}[g_{\theta}(x^{\prime})]. \tag{15}\]

The induced excess variance in Equation (15) can bias the optimal solution. As a result, optimizing the approximated loss will prefer smoother solutions over all \(\{x^{\prime}_{i}\}_{i=1}^{N}\) samples. It is worth noting that this bias is mostly harmful due to its parametrized nature; the only link through which this bias can offset the optimal solution is its dependency on \(\theta\). This is in contrast to any non-parametrized stochasticity in the \(y\) term of Equation (10). Non-parameterized terms cannot offset the optimal solutions, since stochastic gradient descent methods are indifferent to them.

## 3 Proposed Solutions

Based on Equation (15), the induced bias in the solution has a direct relationship with the stochasticity of the conditional distribution \(P(x^{\prime}|x)\). If we were to sample the \((x,x^{\prime})\) pairs deterministically, the excess variance in Equation (15) would disappear. However, this condition is satisfied only by modifying the problem conditions. Next, we introduce three preliminary solutions to this problem: the _deterministic sampling strategy_, the _double-sampling trick_, and the _delayed target method_.

**The deterministic sampling strategy**: One approach to eliminate the excess variance term in Equation (14), is to sample the \(\{x^{\prime}_{i}\}_{i=1}^{N}\) set in a way that \(P(\{x^{\prime}_{i}\}|x)\) would be a point mass distribution at a fixed set \(\mathbb{A}_{x}\). This way, \(P(\{x^{\prime}_{i}\}|x)\) yields a zero excess variance: \(\mathbb{V}[\frac{1}{N}\sum_{i=1}^{N}g_{\theta}(\mathbb{A}_{x}^{(i)})]=0\). This induces the following deterministic loss.

\[\hat{\mathcal{L}}_{\theta}^{\text{DET}}(x):=\bigg{(}f_{\theta}(x)-\frac{1}{N} \sum_{i=1}^{N}g_{\theta}(\mathbb{A}_{x}^{(i)})-y\bigg{)}^{2}. \tag{16}\]

Although this approach removes the excess variance term in Equation (14) thanks to its deterministic nature, it biases the optimization loss by re-defining it: \(\mathcal{L}_{\theta}(x)\neq\hat{\mathcal{L}}_{\theta}^{\text{DET}}(x)\). The choice of the \(\mathbb{A}_{x}\) set can influence the extent of this discrepancy. One reasonable choice is to evenly space the \(N\) samples to cover the entire sampling domain as uniformly as possible. For a demonstration, Figure 2 shows a number of example sets used for applying the divergence theorem to the Poisson problem. Of course, this sampling strategy can be impractical in high-dimensional spaces, as the number of samples needed to cover the entire sampling domain grows exponentially with the sampling space dimension.

**The double-sampling trick**: If we have two independent \(x^{\prime}\) samples, namely \(x^{\prime}_{1}\) and \(x^{\prime}_{2}\), we can replace the objective in Equation (11) with

\[\hat{\mathcal{L}}_{\theta}^{\text{DBL}}(x)=\mathbb{E}_{x^{\prime}_{1},x^{ \prime}_{2}\sim P(x^{\prime}|x)}\bigg{[}\big{(}f_{\theta}(x)-g_{\theta}(x^{ \prime}_{1})-y\big{)}\cdot\big{(}f_{\theta}(x)-g_{\theta}(x^{\prime}_{2})-y \big{)}\bigg{]}. \tag{17}\]

It is straightforward to show that \(\hat{\mathcal{L}}_{\theta}^{\text{DBL}}(x)=\mathcal{L}_{\theta}(x)\); the uncorrelation between \(g_{\theta}(x^{\prime}_{1})\) and \(g_{\theta}(x^{\prime}_{2})\) will remove the induced bias on average. However, this approach requires the access to two i.i.d. samples, which may not be plausible in many sampling schemes. In particular, Monte-Carlo samplings used in reinforcement learning do not usually afford the learning method with freedom to choose

Figure 3: The results of the deterministic and double sampling techniques on the Poisson problem. The left plots demonstrate the solutions with \(N=1\), while the right plots show the solutions with \(N=100\). The training curves represent the mean squared error to the analytical solution vs. the training epochs. With \(N=1\), the double sampling trick exhibits divergence in training, and the deterministic sampling process yields overly-smooth functions similar to the standard solution in Figure 1. However, with \(N=100\), both the deterministic and double-sampling approaches exhibit improvements. According to the training curves, the delayed target method with \(N=1\) yields the best solutions in this problem.
multiple next samples or being able to reset to a previous state. Besides reinforcement learning, offline learning using a given collection of samples may make this approach impractical.

**The delayed target method**: This approach replaces the objective in Equation (11) with

\[\mathcal{L}_{\theta}^{\text{DT}}(x)=\mathbb{E}_{P(x^{\prime}|x)}\bigg{[}\big{(}f _{\theta}(x)-g_{\theta^{*}}(x^{\prime})-y\big{)}^{2}\bigg{]}, \tag{18}\]

where we have defined \(\theta^{*}:=\operatorname*{arg\,min}_{\tilde{\theta}}\mathcal{L}_{\tilde{ \theta}}(x)\). Assuming a complete function approximation set \(\Theta\) (where \(\theta\in\Theta\)), we know that \(\theta^{*}\) satisfies Equation (1) at all points. Therefore, we have

\[\nabla_{\theta}\mathcal{L}_{\theta}(x)\big{|}_{\theta=\theta^{*}}=\nabla_{ \theta}\mathcal{L}_{\theta}^{\text{DT}}(x)\big{|}_{\theta=\theta^{*}}=0. \tag{19}\]

Therefore, we can claim

\[\theta^{*}=\operatorname*{arg\,max}_{\theta}\mathbb{E}_{x}[\mathcal{L}_{ \theta}^{\text{DT}}(x)]=\operatorname*{arg\,min}_{\theta}\mathbb{E}_{x}[ \mathcal{L}_{\theta}(x)]. \tag{20}\]

In other words, optimizing Equation (18) should yield the same solution as optimizing the true objective \(\mathcal{L}_{\theta}(x)\) in Equation (10). Of course, finding \(\theta^{*}\) is as difficult as solving the original problem. The simplest heuristic replaces \(\theta^{*}\) with a supposedly independent, yet identically valued, version of the latest \(\theta\) named \(\theta^{\text{Target}}\), hence the delayed, detached, and bootstrapped target naming conventions:

\[\hat{\mathcal{L}}_{\theta}^{\text{DT}}(x)=\mathbb{E}_{P(\{x^{\prime}_{i}\}|x) }\bigg{[}\big{(}f_{\theta}(x)-\frac{1}{N}\sum_{i=1}^{N}g_{\theta^{\text{target }}}(x^{\prime}_{i})-y\big{)}^{2}\bigg{]}. \tag{21}\]

Our hope would be for this approximation to improve as well as \(\theta\) over training. The only practical difference between implementing this approach and minimizing the loss in Equation (10) is to use an incomplete gradient for updating \(\theta\) by detaching the \(g(x^{\prime})\) node from the computational graph in the automatic differentiation software. This naive implementation of the delayed target method can lead to divergence in optimization, as we will show in Section 4 with numerical examples (i.e., Figure 2). Here, we introduce two mitigation factors contributing to the stabilization of such a technique.

Moving target stabilizationOne disadvantage of the aforementioned technique is that it does not define a global optimization objective; even the average target for \(f_{\theta}(x)\) (i.e., \(\mathbb{E}\big{[}g_{\theta^{\text{target}}}(x^{\prime})\big{]}+y\)) changes throughout the training. Therefore, a naive implementation can risk training instability or even divergence thanks to the moving targets.

Figure 4: The solution and performance curves in higher-dimensional Poisson problems. _The left panel_ shows the solution curves for the delayed target (\(N=1\)), the standard (\(N=100\)), and the double-sampling (\(N=100\)) methods. The top and the bottom rows show 2- and 10-dimensional problems, respectively. In these problems, a single charge is located at the origin, so that the analytical solution is a function of the evaluation point radii \(\|x\|\). The horizontal axis shows the evaluation point radii, and covers 99% of points within the training volumes. _The right chart_ is a performance curve against the problem dimension (lower is better). The normalized MSE values were shown to be comparable. These results suggest that (1) higher dimensions make the problem challenging, and (2) delayed targeting with \(N=1\) is comparable to standard trainings with \(N=100\).
To alleviate the fast moving targets issue, prior work suggested fixing the target network for many time-steps [19]. This causes the training trajectory to be divided into a number of episodes, where the target is locally constant and the training is therefore locally stable in each episode. Alternatively, this stabilization can be implemented continuously using Polyak averaging; instead of fixing the target network for a window of \(T\) steps, the target parameters \(\theta^{\text{Target}}\) can be updated slowly with the rule

\[\theta^{\text{Target}}\leftarrow\gamma\theta^{\text{Target}}+(1-\gamma)\theta. \tag{22}\]

This exponential moving average defines a corresponding stability window of \(T=O(1/(1-\gamma))\).

Prior imposition for highly stochastic targetsIn certain instances, the target \(\frac{1}{N}\sum_{i=1}^{N}g_{\theta^{\text{Inpg}}}(x^{\prime}_{i})+y\) in Equation (18) can be excessively stochastic, leading to divergence in the training of the delayed target model. In particular, based on the settings defined in Equation (5) for the Poisson problem, we can write \(g_{\theta}(x_{i})=(M-1)f_{\theta}(x_{i})\). Therefore, we can analyze the target variance as

\[\mathbb{V}_{\{x^{\prime}_{i}\}_{i=1}^{N}}\big{[}\frac{1}{N}\sum_{i=1}^{N}g_{ \theta^{\text{Inpg}}}(x^{\prime}_{i})+y\big{]}=\frac{(M-1)^{2}}{N}\mathbb{V}_ {\{x^{\prime}_{i}\}_{i=1}^{N}}[f_{\theta^{\text{Inpg}}}(x^{\prime}_{i})]+ \mathbb{V}[y]. \tag{23}\]

Ideally, \(M\rightarrow\infty\) in order for Equation (4) to hold. Setting arbitrarily large \(M\) will lead to unbounded target variances in Equation (23), which can slow down the convergence of the training or result in divergence. In particular, such unbounded variances can cause the main and the target models to drift away from each other leading to incorrect solutions.

One technique to prevent this drift, is to impose a Bayesian prior on the main and the target models. Therefore, to discourage this divergence phenomenon, we regularize the delayed target objective in Equation (24) and replace it with

\[\hat{\mathcal{L}}_{\theta}^{\text{DTR}}:=\hat{\mathcal{L}}_{\theta}^{\text{DT }}(x)+\lambda\cdot(f_{\theta}(x)-f_{\theta^{\text{Inpg}}}(x))^{2}. \tag{24}\]

A formal description of the regularized delayed targeting process is given in Algorithm 1, which covers both of the moving target stabilization and the Bayesian prior imposition.

```
1:Initialize \(\theta^{\text{DTR}}\leftarrow\theta^{\text{DTR}}\)
2:for\(i=1,\dots,N\)do

[MISSING_PAGE_POST]
three unit Dirac-delta charges at \([0,0]\), \([-0.5,-0.5]\), and \([0.5,0.5]\). We also study higher-dimensional Poisson problems with a unit charge at the origin in Figure 4. Our second example looks at finding the magnetic potentials and fields around a current circuit, which defines a singular current density profile \(J\). Finally, to simulate particle evolution dynamics, we consider a Smoluchowski coagulation problem where particles evolve from an initial density. We designed the coagulation kernel \(K\) to induce non-trivial solutions in our solution intervals. We employ a multi-layer perceptron as our deep neural network, using 64 hidden neural units in each layer, and the \(\tanh\) activation function. We trained our networks using the Adam [12] variant of the stochastic gradient descent algorithm under a learning rate of \(0.001\). For a fair comparison, we afforded each method 1000 point evaluations for each epoch. Due to space limitations, a wealth of ablation studies with more datasets and other experimental details were left to the supplementary material.

The Poisson problem with singular chargesTo show the solution bias, we first train two models: one with \(N=100\) samples per sphere, and another one with only \(N=1\) sample per sphere. These models represent a baseline for later comparisons. Based on Equation (15), the induced solution bias should be lower in the former scenario. Figure 1 shows the solution defined by these models along with the analytical solution and their respective training curves. The model trained with high estimation variance derives an overly smooth solution. We hypothesize that this is due to the excess variance in the loss. This hypothesis is confirmed by matching the training loss and the excess variance curves; the training loss of the model with \(N=1\) is lower bounded by its excess variance, although it successfully finds a solution with smaller excess variance than the \(N=100\) model. An alternative capable of producing similar quality solutions with only \(N=1\) samples would be ideal.

To investigate the effect of highly stochastic targets on delayed target models, Figure 2 shows the training results with both \(M=100\) and \(M=10\). The former is unstable, while the latter is stable; this confirms the influence of \(M\) in the convergence behavior of the delayed target trainings. Furthermore, when this divergence happens, a clear drift between the main and the target models can be observed. Figure (2) shows that imposing the Bayesian prior of Equation (24) can lead to training convergence even with a larger \(M=1000\), which demonstrates the utility of our proposed solution.

We also investigated the performance of the deterministic and double-sampling techniques in this problem. Figure 3 shows these results when \(N=1\) and \(N=100\) samples are used for integral estimation. With \(N=1\), the training with the deterministic sampling approach is stable and yields similar results to those seen in Figure (1). The double-sampling trick, on the other hand, exhibits unstable trainings and sub-optimal solutions. We suspect that (a) the singular nature of the analytical

Figure 6: Training results on the Smoluchowski coagulation problem. The top left panel shows the ground truth solution, along with the standard \(N=100\) and \(N=1\) solution heatmaps minimizing the \(\hat{\mathcal{L}}_{\theta}(x)\) in Equation (15). The training loss and the integration variance represent the \(\hat{\mathcal{L}}_{\theta}(x)\) and \(\mathbb{V}_{P(x^{\prime}|x)}[g_{\theta}(x^{\prime})]\) quantities in Equation (15). The top right figure shows the training curve for both of the standard trainings. The bottom left panel shows the delayed target solution heatmaps using \(N=1\) sample with its training curve next to it.
solution, and (b) the stochasticity profile of the training loss function \(\hat{\mathcal{L}}_{\theta}^{\text{DBL}}(x)\) in Equation (17) are two of the major factors contributing to this outcome. With \(N=100\), both the deterministic and double-sampling trainings yield stable training curves and better solutions. This suggests that both methods can still be considered as viable options for training integro-differential PINNs, conditioned on that the specified \(N\) is large enough for these methods to train stably and well.

The regularized delayed target training with \(N=1\) sample is also shown in the training curves of Figure 3 for easier comparison. The delayed target method yields better performance than the deterministic or double-sampling in this problem. This may seemingly contradict the fact that the double-sampling method enjoys better theoretical guarantees than the delayed target method, since it optimizes a complete gradient. However, our results are consistent with recent findings in off-policy reinforcement learning; even in deterministic environments where the application of the double-sampling method can be facilitated with a single sample, incomplete gradient methods (e.g., TD-learning) may still be preferable over the full gradient methods (e.g., double-sampling) [24, 5, 31, 4]. Intuitively, incomplete gradient methods detach parts of the gradient, depriving the optimizer from exercising full control over the decent direction and make it avoid over-fitting. In other words, incomplete gradient methods can be viewed as a middle ground between zero-order and first-order optimization, and may be preferable over both of them.

Figure 4 also studies the effect of problem dimensionality on our methods. The results confirm that the problem becomes significantly more difficult with higher dimensions. However, the delayed target models maintain comparable quality to standard trainings with large \(N\).

The Maxwell problem with a rectangular current circuitFigure 5 shows the training results for the Maxwell problem. The results suggest that the standard and the deterministic trainings with small \(N\) produce overly smooth solutions. The double-sampling method with small \(N\) improves the solution quality at first, but has difficulty maintaining a stable improvement. However, delayed targeting with small \(N\) seems to produce comparable solutions to the standard training with large \(N\).

The Smoluchowski coagulation problemFigure 6 shows the training results for the Smoluchowski coagulation problem. Similar to the results in Figure 1, the standard training using \(N=1\) sample for computing the residual summations leads to biased and sub-optimal solution. However, the standard training with \(N=100\) samples suffers less from the effect of bias. The delayed target solution using only \(N=1\) sample produces comparable solution quality to the standard evaluation with \(N=100\) and is not bottlenecked by the integration variance. Figure 7 compares the solution quality for each of the standard and delayed target methods under different problem dimensions. The results suggest that the delayed target solution maintains its quality even in higher dimensional problems, where the excess variance issue leading to biased solutions may be more pronounced.

Figure 7: The solution mean squared error to the ground truth in the 1, 2, and 3-dimensional Smoluchowski coagulation problem. The vertical axis shows the solution error, and the horizontal axis shows the training epochs. The standard solutions were trained by the ordinary MSE loss \(\mathcal{L}_{\theta}(x)\) in Equation (10) with \(N=1\) and \(N=100\) samples. The delayed target solution used \(N=1\) sample, yet produced slightly better results than the standard method with \(N=100\).
Conclusion

In this work, we investigated the problem of learning PINNs in partial integro-differential equations. We presented a general framework for the problem of learning from integral losses, and theoretically showed that naive approximations of the parametrized integrals lead to biased loss functions due to the induced excess variance term in the optimization objective. We confirmed the existence of this issue in numerical simulations. Then, we proposed three solutions to account for this issue. In particular, we detailed a delayed targeting recipe for this class of problems, and verified that it can solve these problems effectively without access to too many samples for integral evaluations. Our numerical results support the utility of our proposed method on three challenging problems, (1) a Poisson problem with singular charges, (2) an electromagnetic problem under a Maxwell equation, and (3) a Smoluchowski coagulation problem. The limitations of our work include its narrow scope to learning PINNs; this work could have broader applications in other areas of machine learning. Also, future work should consider the applications of the proposed method to more problem classes both in scientific and traditional machine learning areas. Developing adaptive processes for setting each method's hyper-paramters is another worthwhile future endeavor.
