# DNN Filter for Bias Reduction in Distribution-to-Distribution Scan Matching

Matt McDermott
Tufts University
Medford MA, USA
<EMAIL_ADDRESS>Jason Rife
Tufts University
Medford MA, USA
<EMAIL_ADDRESS>
###### Abstract

Distribution-to-distribution (D2D) point cloud registration techniques such as
the Normal Distributions Transform (NDT) can align point clouds sampled from
unstructured scenes and provide accurate bounds of their own solution error
covariance– an important feature for safety-of-life navigation tasks. D2D
methods rely on the assumption of a static scene and are therefore susceptible
to bias from range-shadowing, self-occlusion, moving objects, and distortion
artifacts as the recording device moves between frames. Deep Learning-based
approaches can achieve higher accuracy in dynamic scenes by relaxing these
constraints, however, DNNs produce uninterpretable solutions which can be
problematic from a safety perspective. In this paper, we propose a method of
down-sampling LIDAR point clouds to exclude voxels that violate the assumption
of a static scene and introduce error to the D2D scan matching process. Our
approach uses a solution consistency filter– identifying and suppressing
voxels where D2D contributions disagree with local estimates from a PointNet-
based registration network. Our results show that this technique provides
significant benefits in registration accuracy, and is particularly useful in
scenes containing dense foliage.

Figure 1: Comparison of registration of points from two scans captured at
locations 1 and 2, contained within voxels a and b. The DNN registration
estimate disagrees with the D2D estimate for the registration of the head, so
the points in voxel a will be ignored.

##  1 Introduction

Our goal in this paper is to to demonstrate a technique for identifying and
discarding voxels that introduce error in Distribution-to-Distribution point
cloud registration. In order to discuss our proposed method we must first
introduce the current state of the art in point cloud registration, provide
insight into some of the challenges associated with this task, and discuss how
they relate to our proposed approach.

Distribution-to-Distribution methods represent clusters of points as
probability density functions and use Maximum Likelihood Estimation to find
the transformation that best aligns the two PDFs [3, 22]. The accuracy of a
LIDAR point cloud registration is a function of both the quality of the sensor
and the geometry of the scene in which the data is recorded. Analytical
techniques permit the direct calculation of error bounds for solution
covariance by considering the cost function at convergence [4, 22] or directly
via a least-squares analysis of the PDF occupying each voxel [14]. While these
registration techniques boast impressive accuracy prediction on static scenes,
motion of the ego-vehicle can introduce systemic bias through range-shadowing
[15], self-occlusion [27], and ”rolling shutter” motion distortion [9]. For
scan registration algorithms to be useful in real-world applications, they
must maintain robustness in the presence of these forms of bias.

Deep Learning methods can achieve superior accuracy to D2D techniques by
relaxing some of their assumptions and forfeiting solution interpratability.
Deep Neural Network (DNN) based point cloud registration can be broadly
divided between correspondence-free methods, in which a DNN directly estimates
a transformation in an end-to-end manner, and correspondence-based methods,
which alternate between identifying corresponding features in each scan and
minimizing distances between correspondences [28] similar to the early
Iterative Closest Point (ICP) algorithm [2]. Correspondence-based methods can
use a DNN to identify keypoints in each frame before estimating
correspondences analytically [10] or can identify correspondences using an
additional network [11]. Generally, correspondence-based methods achieve
higher accuracy registering larger point clouds (such as full LIDAR scans)
while correspondence-free methods excel on clouds sampled from individual
objects or small portions of a LIDAR scan [28]. One common approach is to
subdivide large LIDAR scans into voxels and apply a single end-to-end network
to each cell simultaneously as a batch operation [1, 29]. Some batch methods
also aggregate local estimates to gain additional contextual information about
the scene [24].

Partial overlap between point clouds, which occurs in regions subject to
range-shadowing and self-occlusion, presents a problem for registration.
Range-shadowing tends to inject systemic bias across a large swath of voxels
in a constant direction normal to the edge of each shadow [15], while
perspective shift bias is highly localized on the surfaces of objects that
appear different from various angles [27]. Geometric strategies exist to
reduce the bias of range-shadowing in D2D algorithms, either through ground
plane removal [19, 7] or through a spherical coordinate representation [9, 15,
5], however, error introduced by self-occlusion is significantly more complex.
One approach to mitigate self-occlusion is to train a network to fill in
occluded regions of point clouds before registration is attempted. In MaskNet,
the authors directly estimate inlier points from one scan using candidate
points from the other [17]. Self-supervised methods overcome partial overlap
by jointly learning geometric representations, keypoint detection, and
correspondence [25, 30]. BTCdet augments preliminary estimates from a backbone
network by filling in points on occluded surfaces using a network trained a
priori on commonly encountered object shapes [27].

Although there is a significant body of literature on both Machine Learning-
based and Distribution-to-Distribution point cloud registration methods, there
is limited work in the context of leveraging ML to improve distribution-based
approaches. Tang et al. use Gaussian Process Regression to reduce bias in NDT
by directly estimating a correction factor to add to their preliminary
registration [23]. While their strategy is capable of reducing registration
error, allowing an uninterpretable model to directly adjust odometry estimates
negates much of the safety benefits of NDT. Our proposed solution differs from
this form of mixed method in that the uninterpretable component solely exists
to down-sample data that is likely to introduce bias to the geometric
registration routine. Crucially, the reduction of total information due to the
down-sampled point cloud can be accounted for in the optimization routine,
automatically adjusting the solution error covariance prediction component of
the registration algorithm as problematic voxels are removed from a scan.

Random Sample Consensus (RANSAC) is a general statistical tool for outlier
removal [21], and has been used for ground plane segmentation in NDT [8].
Kanhere and Gao also use measurement consensus to improve computational
efficiency in NDT, and to devise a metric for quantifying localization
reliability [6]. Other approaches assume the magnitude of residuals at a
converged solution will be normally distributed and simply discard any voxels
that produce residuals more than a few standard deviations away from the mean
[13]. Unfortunately, classical outlier routines rely on a large plurality of
voxels in a scene being well-conditioned. This can lead to these approaches
excluding voxels with accurate solution contributions if they are outnumbered
by voxels with consistently biased estimates. The density of the points within
a voxel is inversely proportional to the distance of the voxel from the
origin, so the regions of a scan that are most susceptible to perspective
shift bias (close to the origin) are also on average the most populated and
heavily weighted.

Moving objects within the frame also present additional challenges for scan
registration. A fine voxel gird will automatically discount small objects
moving at high speed if they do not occupy the same cell at the converged
solution. While our strategy of considering per-voxel registration to identify
bias from occlusion does not directly suppress moving objects, it allows the
rejection threshold of our moving object detection routine to be raised to to
the point that its only function is to discard obviously moving objects.

In section 2 we provide more detail on the detection theory underpinning our
filtering process. In section 3 we describe our network structure and training
procedure. In sections 4 and 5 we introduce and share the results of two
experiments to validate our method; a detection analysis where we demonstrate
how our approach is used to identify and exclude biased voxels, and a
demonstration on various KITTI datasets where we show how our filter preforms
on real LIDAR data. In section 6 we provide further discussion into our
results. Lastly, in section 7 we summarize our findings and provide insights
for future work.

##  2 Filtering

Our approach is to test for bias in an estimate without having access to the
true solution, introducing a secondary estimation technique and considering
the difference between the two. If there is agreement between the two
uncorrelated estimates, it is unlikely that either is significantly flawed. If
the two estimates differ by a sufficient magnitude, however, one of the two
registration techniques is likely biased and the points inside the offending
voxel should be ignored. As long as the errors between D2D and DNN remain
uncorrelated, pruning voxels where the two methods disagree provides the
flexibility of a learned registration scheme while retaining statistical
robustness.

Because the translation errors from DNN and D2D solutions in the $x$, $y$, and
$z$ components are all independent and normally distributed (as we will
demonstrate in the following section), the overall magnitude of error in the
translation vector for each method relative to the ground truth forms a chi-
square distribution with three degrees of freedom. This relationship is
described by equation (1), where $x$, $y$, and $z$ represent the true
translation solutions about each axis and $\hat{x}$, $\hat{y}$, and $\hat{z}$
represent the associated estimates output by each technique.

$Q_{3}\sim(x-\hat{x})^{2}+(y-\hat{y})^{2}+(z-\hat{z})^{2}$ (1)

When constructing our monitor statistic, however, we require a single metric
to summarize the spread of residuals between the DNN and D2D solutions. It is
convenient to difference the two measurements in vector space before unifying
solution components into a single chi-square distribution as the sum of two
normally distributed random variables is also normal, while calculating the
difference of noncentral chi-square distributions is nontrivial. Here, the
vector difference between D2D and DNN estimates (and consequently, the
difference between their errors) is represented as $\Delta$. The distribution
of these errors for a given voxel registration is then defined as
$Q_{3,\Delta}$:

$\begin{split}\mathbf{E}_{DNN}=\mathbf{X}_{DNN}-\mathbf{X}_{true}\\\
\mathbf{E}_{D2D}=\mathbf{X}_{D2D}-\mathbf{X}_{true}\\\
\Delta=\mathbf{E}_{D2D}-\mathbf{E}_{D2D}\\\ Q_{3,\Delta}=|\Delta|\end{split}$
(2)

Figure (2) displays the distribution of $Q_{3,\Delta}$ for two cases. The top
plot demonstrates the spread of $Q_{3,\Delta}$ when neither the D2D or DNN
approach is systemically biased for a particular test sample (such as the wall
inside voxel b in Fig. (1)). The bottom plot in Fig. (2) illustrates an
alternate case of $Q_{3,\Delta}$ on a test sample that produces significant
systemic bias in the D2D solution, as would be the case in the presence of
heavy shadowing or self-occlusion such as the partially occluded human model
inside voxel a in Fig. (1). In Fig. (2), the shaded region within each plot
represents the areas in which the monitor statistic fails to correctly
characterize the behavior of the system. This demonstrates the trade-off
between false alarms and missed detection when selecting a threshold for the
alert limit. The alert limit in figure (2) is set to produce a 5% chance of
false alarm (shown in the top plot) where the DNN incorrectly excludes an
unbiased D2D solution estimate. At this threshold, the lower plot represents
the associated risk of missed detection for test sample where D2D registration
is biased by 5cm. In this example, the selected threshold value produces a
significantly higher level of false negatives than false positives. Framing
outlier rejection as a regression problem, rather than one of binary
classification, allows rejection criteria to be tuned after training depending
on the desired characteristics of the system.

Figure 2: Trade-off between false alarms and missed detections

For both cases in Fig. (2), it is apparent that relying on differences between
two solutions misses a large portion of the useful information. If the
geometic approach was prone to uniformly distributed errors for any scene,
there would be no reason to employ an alternate registration technique (and
thus lose useful info). As is discussed in [15, 27], however, certain
geometric properties in a scene violate the assumption of a static scene
significantly enough to introduce systemic bias. This highlights the inherent
tradeoff between solution accuracy and solution integrity. It may be
acceptable for a noisy filter to reject a large portion of the voxels in a
well-conditioned scene, thereby slightly reducing average registration
accuracy, if the filter is also likely to reject significant bias in other
scenes that would otherwise prevent convergence altogether.

##  3 DNN Structure

Our network is constructed as a variant of Iterative PCR-Net [18], which
leverages PointNet [16] encoding to directly estimate the transformation that
best registers two point clouds. Here, input point clouds are concatenated and
passed through feed-forward layers with shared weights before a Max-Pooling
symmetric function is applied to produce a feature vector, which is then
passed through a series of simple feedforward layers.

Figure 3: DNN structure

The purpose of our network is to either validate or reject the per-voxel D2D
solutions, so we limit the network output to a similar solution space and only
estimate translation within each voxel. Removing rotation from the solution
vector prevents the network from having to contend with any forms of
rotational ambiguity (i.e. when attempting to register an isolated cylinder)
that may detract from the training process. Linearizing each voxel to a
translation-only estimate is consistent with other registration algorithms
such as the Lucas & Kanade (LK) algorithm [12], and its application to the
point cloud registration task with a global PointNet in PointNetLK [1]. A side
benefit of applying our network on a per-voxel basis is that any out of
distribution behavior is compartmentalized to the offending region of the
frame.

Like other local point cloud registration techniques, an initial estimate
within the basin of attraction is required to begin the scan matching process.
This can be achieved using an estimate from another sensor, such as wheel
odometry, a dynamic vehicle model, or simply by zeroing out initial
transformation and attempting a coarse-to-fine approach with large voxels. In
our implementation, a preliminary solution estimate is first achieved by
allowing D2D registration to converge before sampling points to feed to the
DNN.

##  4 Experiments

We validate our method of bias mitigation in distribution-based scan matching
through two experiments

###  4.1 Per-voxel accuracy

Figure 4: 100 point keyframe clouds for uniformly sampled ModelNet40 object
(top), simulated LIDAR scan recorded from from ModelNet40 object (middle), and
voxelized KITTI (bottom)

Our first experiment seeks to quantify the differences in registration
accuracy between the two approaches. Here, we compare solutions from our
PointNet against the distance between the means of each sub-cloud, which is
the basis for the D2D registration contributions for each voxel. The first two
trials of this experiment test simple point clouds generated from objects in
the popular ModelNet40 dataset [26], and represent cases where each voxel
completely encapsulates a single object with no ground plane or external-
occlusion present. The first is an unrealistic case where clouds are formed by
uniformly sampling points on the surface of each object. The second repeats
this process, but generates more realistic point clouds from the ModelNet40
objects by randomly transforming each model and simulating LIDAR scans from
the visible surfaces. This test provides a middle ground between the
unrealistic point clouds in the first trial and real world data in the
following trial.

In the third trial we consider the per-voxel registration accuracy of sub-
clouds sampled from voxels in the KITTI dataset. Here, we gather training data
from each scan using spherical voxels as described by [15] to mitigate the
effects of range-shadowing. Eliminating range-shadows in the training data
greatly reduces the need to share information between cells. For all cases we
train the network on point clouds containing 100 points from each scan per
voxel, a value tuned to align with the optimal minimum number of points
required for D2D registration [15, 7]. Training on small sections of larger
point clouds requires additional considerations, which are discussed in detail
in the following section.

###  4.2 Improving Total Error with Filter

Our second experiment makes use of real LIDAR data, again taken from the KITTI
benchmark dataset. Here, LIDAR data contains shadowing, occlusion, and non-
returns. Surfaces in this test also extend beyond the limits of each cell.
Estimating local registration solutions with no consideration of the larger
context can hinder training, as the information present within a cell is often
insufficient to properly estimate translation [14, 20]. Generally, the spread
of points in the direction normal to a surface is randomly distributed due to
sensor noise, while the spread of points tangent to the surface tends to be
more uniform and structural, following the underlying shape of the object,
rather than the properties of the sensor. Since distributions are cut
arbitrarily short in surface-aligned directions by the voxel grid, attempting
to align centers of these ”overly extended” directions provides no useful
information, and in fact, can negatively impact solution accuracy [14]. For
that reason, we integrate a process of ambiguous axis suppression into our
loss function when training the registration network using equation (4). As
shown in Fig. (5), the our filter is also constructed to only count components
of residuals normal to compact directions of each scan.

Figure 5: Overhead view of a voxel containing partially partially overlapping
sections of the same wall captured in two LIDAR scans. The true translation to
align the clouds is drawn in yellow, however, there is insufficient
information to constrain solutions in the direction aligned with the surface
of the wall. The reduced-dimension ground truth obtained using using Eq. (4)
is shown in purple.

${}^{(j)}\textbf{x}_{D2D}={}^{(j)}\textbf{L}\>{}^{(j)}\textbf{U}^{T}({}^{(j)}\bar{\textbf{x}}_{r}-{}^{(j)}\bar{\textbf{x}}_{n})$
(3)

Equation (3) describes the process of projecting the vector representing the
reference and new scan means $\bar{\mathbf{x}}_{r}$ and
$\bar{\mathbf{x}}_{n}$, respectively for each voxel $j$ into a reduced
dimension residual vector. Here, $\mathbf{U}$ is the rotation matrix (obtained
via eigendecomposition of the covariance of the points in voxel $j$) that
aligns the principal axis of distribution $J$ relative to the world frame, and
$\mathbf{L}$ is an truncation matrix.

${}^{(j)}\textbf{x}_{DNN}={}^{(j)}\textbf{L}\>{}^{(j)}\textbf{U}^{T}{}^{(j)}\hat{\textbf{x}}_{DNN}$
(4)

Equation (4) represents the process through which extended components of the
residual vector estimated by the DNN are excluded. Here,
${}^{(j)}\hat{\textbf{x}}_{DNN}$ represents the raw translation estimate
output by the network for voxel $j$. It is important to note how this process
also makes use of the same $\mathbf{U}$ and $\mathbf{L}$ matrices as equation
(3). This means that the shape of the keyframe distribution (that has already
been calculated in the classical distribution-to-distribution routine) is used
to inform the pruning matrix which components of the DNN residual are
important to consider, because they likely represent a true offset between
surfaces, and which components are merely due to a structural spread of
points. The total offset, ${}^{(j)}\Delta\mathbf{x}$, between useful
components of the D2D and DNN solution is then calculated as follows

${}^{(j)}\Delta\mathbf{x}=\lvert{}^{(j)}\mathbf{x}_{DNN}-{}^{(j)}\mathbf{x}_{D2D}\rvert$
(5)

If ${}^{(j)}\Delta\mathbf{x}$ is larger than a threshold $T$, all points
located inside voxel $j$ are excluded and the associated contributions to the
overall solution vector from voxel $j$ are ignored. When constructing our DNN,
we achieved best results by applying the projection and dimension reduction
matrices inside the loss function, rather than as distinct input parameters
passed to the network. This allows us to train using data from both KITTI, as
well as simulated LIDAR scans of objects in the ModelNet40 dataset.

In Experiment II we validate our strategy on three KITTI datasets. The first
contains 160 sequential LIDAR point clouds recorded from a vehicle moving
between 3 and 7 $\frac{m}{s}$ in a busy urban setting. The second dataset
contains 194 sequential frames where the vehicle moves between 11 and 22
$\frac{m}{s}$ on a highway through a dense forest. The third dataset contains
800 sequential LIDAR frames recorded while moving through a forested
residential neighborhood at 8 $\frac{m}{s}$.

##  5 Results

Table 1: Experiment I: Per Voxel RMS Translation Error (cm) |  DNN (ours) |  D2D
---|---|---
ModelNet40 (uniformly sampled) |  5.48 |  3.61
ModelNet40 (shadowed) |  5.82 |  15.29
KITTI |  9.79 |  14.04

Figure 6: D2D and DNN results produce similar levels of error on uniformly
sampled ModelNet40 data (left). DNN achieves significantly higher accuracy
than D2D on shadowed ModelNet40 point clouds (center). DNN acheives higher
accuracy than D2D on KITTI data (right)

Figure 7: Correlation between Magnitude Difference between D2D/DNN estimates and unobservable error of D2D estimate for uniformly sampled ModelNet40 data (left), simulated LIDAR sampled ModelNet40 data (center), and real LIDAR data from KITTI dataset (right) Table 2: Experiment II: RMS Forward Translation Error (cm) |  no voxel rejection |  2$\sigma$ rejection only |  2$\sigma$ & PointNet rejection (ours)
---|---|---|---
Scene 1 (city) |  1.47 |  1.32 |  1.34
Scene 2 (forest) |  5.09 |  4.63 |  4.45
Scene 3 (residential) |  3.70 |  2.75 |  2.60

The first trial of Experiment I compared the two methods on uniformly sampled
point clouds generated from ModelNet40 objects. In this trial, the D2D method
outperformed the DNN solution. Interestingly however, the DNN achieved almost
the same accuracy in the simulated LIDAR data as in the uniform sampling,
while the D2D estimate performed significantly worse when attempting to
register partial point clouds. In the third trial, our DNN outperformed the
D2D again, though by a slightly smaller margin.

In Experiment II, the city, forest, and residential trials all used the same
network and NDT parameters. For each scene, we compared root-mean-square error
of the global NDT registration against a GNSS/INS baseline. Our first test
case consists of the raw NDT registration with no rejected voxels. We then
consider the commonly employed strategy of removing any voxel with a local D2D
residual greater than two standard deviations from the mean residual distance
after a preliminary convergence, as prescribed by [13]. Finally, we consider
the case where our filter is used in conjunction with the 2$\sigma$ rejection
criteria. For the city test scene, our proposed method achieved slightly worse
accuracy than the 2$\sigma$-only rejection criteria by roughly 1%. The raw
baseline does worse than either, presumably because moving objects are
preventing the registration algorithm from converging on the correct solution.
In the forest scene, however, our PointNet + 2$\sigma$ combined approach does
better than the 2$\sigma$-only method by a margin of roughly 4%. Similarly,
our DNN based approach outperforms the 2$\sigma$-only rejection method in the
residential dataset by a margin of over 5.5%. Furthermore, because the errors
reported in Table (2) are relative to an imperfect baseline the true odometry
error is likely lower in all cases, making the true percentage reduction in
error from our proposed method even more significant.

##  6 Discussion

The results in Table (1) and Fig. (6) demonstrate our PointNet-based
registration network achieves significantly better accuracy than a simple D2D
registration on realistic point cloud data, but slightly worse than the D2D
strategy on the unrealistic uniformly sampled ModelNet40 data. The D2D method
performed very well in the first trial due to the high degree of symmetry in
the data, but performed far worse on the subsequent trials. Our DNN
performance, on the other hand, remained relatively insensitive to using full
or partial point clouds, achieving similar accuracy between trials 1 and 2.
This is likely because the feature vector formed by our PointNet can infer the
contours of an object without requiring total overlap between the two scans.

Figure 8: KITTI Raw 0005 (top), KITTI Raw 0027 (center), and KITTI Odometry 03
(bottom)

The results in Table (2) demonstrate that performance enhancements provided by
our DNN filter are highly dependant on the driving environment. For the urban
case, there is little benefit to removing voxels with our network, while the
forest and residential scenes achieve substantial benefit from our proposed
technique. This is possibly due to the fact that in urban settings, the
influence of self-occlusion prone surfaces on constraining translation in the
forward direction is dwarfed by other large, highly-localizable surfaces, such
as the sides of buildings. The 2$\sigma$ rejection method does a good job of
identifying moving objects when surfaces are well defined (such as in the city
trial) because the global standard deviation of residuals at convergence is
small. In more ill-defined situations (such as those containing lots of dense
foliage), the standard deviation of residual error at a converged solution is
so large that it is unable to reject obviously biased voxels.

Figure 9: Sequential point clouds from the city (top) and forest (bottom)
trajectories registered by our algorithm. Occupied spherical voxels drawn
containing distributions from the keyframe scan. Voxels identified as
containing moving objects (via the 2$\sigma$ rejection filter) are highlighted
in yellow. Voxels flagged as biased by our DNN filter are highlighted in
purple.

The scanning noise normal to surfaces (either due to time of flight error from
the sensor or surface roughness) also plays an important role in the
performance benefits of our proposed bias mitigation strategy. In our testing,
we noted that the ”fuzzier” the surfaces, the more beneficial our bias
mitigation procedure became. Recall, D2D registration techniques use the
mahalanobis distance to weigh residuals between corresponding distribution
ellipsoids [22]. Scenes with higher surface normal noise lead to situations
where flat surfaces (such the side of the building in Fig. (8) are not oblate
enough to sufficiently resist the influence of high perspective shift bias in
other regions of the frame. The surfaces that tend to produce the highest
perspective shift error are either concave or convex, and thus generate
distribution ellipsoids that are far less oblate (and therefore weighted less)
than their well conditioned counterparts, by virtue of a gaussian’s poor
ability to characterize a non-planar surfaces. When overall noise is low,
these features are heavily deweighted. On the other hand, higher overall noise
effectively raises the lower bound on compactness for any distribution of
points, thereby deweighting the contributions of flat, well-conditioned
surfaces and allowing the biased portions of the scan to dominate
registration.

Bias from self-occlusion is also highly dependant on vehicle speed. This is
because sequential point clouds are recorded from further away from one
another as the vehicle speeds up. At large distances between frames, small
objects close to the vehicle are observed from different angles– capturing
different faces of the same object violates the assumption of a static scene
and can cause the voxel to incorrectly attempt to register non-corresponding
parts of the object. This effect is particularly exaggerated in the forest
experiment.

Fig. (9) presents a visualization of two pairs of point clouds and highlights
the voxels flagged as potentially biased. In the forest scene, the flagged
voxels are primarily aligned close to the sides of the vehicle and normal to
the direction of travel. This is where perspective shift bias is most severe–
the fact that these voxels are being suppressed explains the performance
benefits of our network in forested scenes. In the city scene, however, the
voxels flagged for removal are fewer in number and spread much more uniformly
throughout the frame. It is likely a large number of these identified voxels
are false positives rather than systemic error so removing them from the
registration estimate ignores a small amount of useful information, which
explains the slight decrease in solution accuracy. Trading solution accuracy
in well structured environments for solution integrity in unstructured and
novel environments is useful for safety critical applications. Furthermore,
future work in increasing the accuracy of the registration network can
decrease the number of erroneously removed voxels.

The implementation demonstrated in this paper makes use of a simple
registration network trained on a relatively small dataset using a single
NVIDIA RTX 3090. Our network structure only considers the cluster of points
inside each voxel, ignoring any potential queues from their context in the
larger frame. Despite these limitations, our method manages to reduce error in
D2D scan registrations in unstructured forest and residential scenes, and does
not significantly impact accuracy when applied in a well-structured city
environment. Future work in enhancing network structure, especially
aggregating contextual information between voxels has the potential to further
increase overall registration accuracy.

##  7 Conclusion

In this paper we discussed a new method for reducing bias in Distribution-to-
Distribution scan registration. Our proposed technique removes portions of
LIDAR scans that violate the assumption of a static scene by identifying and
discarding voxels that locally disagree with a PointNet-based registration
network. Our method achieves results similar to existing classical outlier
rejection techniques on well structured city scenes and outperforms existing
techniques in unstructured forest and residential environments. Compared to
other machine learning-based bias mitigation techniques, our strategy is
unique in that it does not directly provide a correction factor, and instead
merely removes data from the classical optimization routine that is likely to
introduce error. This allows the registration algorithms to adjust estimates
for solution error covariance as data is removed from the frame, and allows
the scan matching process to remain fully interpretable.

## References

  * [1] Yasuhiro Aoki, Hunter Goforth, Rangaprasad Arun Srivatsan, and Simon Lucey.  Pointnetlk: Robust & efficient point cloud registration using pointnet.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [2] P. J. Besl and N. D. McKay.  A method for registration of 3-d shapes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(2):239–256, 1992.
  * [3] P. Biber and W. Straßer.  The normal distributions transform: A new approach to laser scan matching.  In IEEE International Conference on Intelligent Robots and Systems, volume 3, pages 2743 – 2748 vol.3, 11 2003.
  * [4] Andrea Censi.  An accurate closed-form estimate of icp’s covariance.  In Proceedings 2007 IEEE international conference on robotics and automation, pages 3167–3172. IEEE, 2007.
  * [5] Ali Hassani and Mathieu Joerger.  A new point-cloud-based lidar/imu localization method with uncertainty evaluation.  In Proceedings of the 34th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2021), pages 636–651, 2021.
  * [6] Ashwin Vivek Kanhere and Grace Xingxin Gao.  Lidar slam utilizing normal distribution transform and measurement consensus.  In Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2019), pages 2228–2240, 2019.
  * [7] Hyungjin Kim, Seungwon Song, and Hyun Myung.  Gp-icp: Ground plane icp for mobile robots.  IEEE Access, 7:76599–76610, 2019.
  * [8] Lin Li, Fan Yang, Haihong Zhu, Dalin Li, You Li, and Lei Tang.  An improved ransac for 3d point cloud plane segmentation based on normal distribution transformation cells.  Remote Sensing, 9(5):433, 2017.
  * [9] You Li, Clément Le Bihan, Txomin Pourtau, Thomas Ristorcelli, and Javier Ibanez-Guzman.  Coarse-to-fine segmentation on lidar point clouds in spherical coordinate and beyond.  IEEE Transactions on Vehicular Technology, 69(12):14588–14601, 2020\.
  * [10] Zhichao Li and Naiyan Wang.  Dmlo: Deep matching lidar odometry.  In 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 6010–6017. IEEE, 2020.
  * [11] Weixin Lu, Guowei Wan, Yao Zhou, Xiangyu Fu, Pengfei Yuan, and Shiyu Song.  Deepvcp: An end-to-end deep neural network for point cloud registration.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 12–21, 2019.
  * [12] Bruce D. Lucas and Takeo Kanade.  An iterative image registration technique with an application to stereo vision.  In IJCAI, 1981.
  * [13] Martin Magnusson.  The three-dimensional normal-distributions transform: an efficient representation for registration, surface analysis, and loop detection.  PhD thesis, Örebro universitet, 2009.
  * [14] Matthew McDermott and Jason Rife.  Enhanced laser-scan matching with online error estimation for highway and tunnel driving.  In Proceedings of the 2022 International Technical Meeting of The Institute of Navigation, pages 643–654, 2022.
  * [15] Matthew McDermott and Jason Rife.  Mitigating shadows in lidar scan matching using spherical voxels.  IEEE Robotics and Automation Letters, 7(4):12363–12370, 2022.
  * [16] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  arXiv preprint arXiv:1612.00593, 2016.
  * [17] Vinit Sarode, Animesh Dhagat, Rangaprasad Arun Srivatsan, Nicolas Zevallos, Simon Lucey, and Howie Choset.  Masknet: A fully-convolutional network to estimate inlier points.  In 2020 International Conference on 3D Vision (3DV), pages 1029–1038. IEEE, 2020.
  * [18] Vinit Sarode, Xueqian Li, Hunter Goforth, Yasuhiro Aoki, Rangaprasad Arun Srivatsan, Simon Lucey, and Howie Choset.  Pcrnet: Point cloud registration network using pointnet encoding.  arXiv preprint arXiv:1908.07906, 2019.
  * [19] Tixiao Shan and Brendan Englot.  Lego-loam: Lightweight and ground-optimized lidar odometry and mapping on variable terrain.  In 2018 IEEE/RSJ Int. Conf. Intelligent Robots and Systems (IROS), pages 4758–4765, 10 2018.
  * [20] Shinsuke Shimojo, Gerald H Silverman, and Ken Nakayama.  Occlusion and the solution to the aperture problem for motion.  Vision research, 29(5):619–626, 1989.
  * [21] Sarvesh Kumar Singh, Simit Raval, and Bikram Pratap Banerjee.  Automated structural discontinuity mapping in a rock face occluded by vegetation using mobile laser scanning.  Engineering Geology, 285:106040, 2021.
  * [22] Todor Stoyanov, Martin Magnusson, Henrik Andreasson, and Achim J Lilienthal.  Fast and accurate scan registration through minimization of the distance between compact 3d ndt representations.  The International Journal of Robotics Research, 31(12):1377–1393, 2012.
  * [23] Tim Y. Tang, David J. Yoon, Francois Pomerleau, and Timothy D. Barfoot.  Learning a bias correction for lidar-only motion estimation.  In 2018 15th Conference on Computer and Robot Vision (CRV), pages 166–173, 2018.
  * [24] Mikaela Angelina Uy and Gim Hee Lee.  Pointnetvlad: Deep point cloud based retrieval for large-scale place recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4470–4479, 2018.
  * [25] Yue Wang and Justin M Solomon.  Prnet: Self-supervised learning for partial-to-partial registration.  Advances in neural information processing systems, 32, 2019.
  * [26] Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao.  3d shapenets: A deep representation for volumetric shapes.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2015.
  * [27] Qiangeng Xu, Yiqi Zhong, and Ulrich Neumann.  Behind the curtain: Learning occluded shapes for 3d object detection.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 2893–2901, 2022.
  * [28] Zhiyuan Zhang, Yuchao Dai, and Jiadai Sun.  Deep learning based point cloud registration: an overview.  Virtual Reality & Intelligent Hardware, 2(3):222–246, 2020.  3D Visual Processing and Reconstruction Special Issue.
  * [29] Yin Zhou and Oncel Tuzel.  Voxelnet: End-to-end learning for point cloud based 3d object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4490–4499, 2018.
  * [30] Zhicheng Zhou, Cheng Zhao, Daniel Adolfsson, Songzhi Su, Yang Gao, Tom Duckett, and Li Sun.  Ndt-transformer: Large-scale 3d point cloud localisation using the normal distribution transform representation.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 5654–5660, 2021.


# Improved Training of Physics-Informed
Neural Networks with Model Ensembles

Katsiaryna Haitsiukevich & Alexander Ilin
Department of Computer Science
Aalto University
Espoo, Finland
<EMAIL_ADDRESS>

###### Abstract

Learning the solution of partial differential equations (PDEs) with a neural
network is an attractive alternative to traditional solvers due to its
elegance, greater flexibility and the ease of incorporating observed data.
However, training such physics-informed neural networks (PINNs) is notoriously
difficult in practice since PINNs often converge to wrong solutions. In this
paper, we address this problem by training an ensemble of PINNs. Our approach
is motivated by the observation that individual PINN models find similar
solutions in the vicinity of points with targets (e.g., observed data or
initial conditions) while their solutions may substantially differ farther
away from such points. Therefore, we propose to use the ensemble agreement as
the criterion for gradual expansion of the solution interval, that is
including new points for computing the loss derived from differential
equations. Due to the flexibility of the domain expansion, our algorithm can
easily incorporate measurements in arbitrary locations. In contrast to the
existing PINN algorithms with time-adaptive strategies, the proposed algorithm
does not need a pre-defined schedule of interval expansion and it treats time
and space equally. We experimentally show that the proposed algorithm can
stabilize PINN training and yield performance competitive to the recent
variants of PINNs trained with time adaptation.

_K_ eywords Label propagation $\cdot$ Model ensembles $\cdot$ Partial
differential equations $\cdot$ Physics-informed neural networks

##  1 Introduction

Partial differential equations (PDEs) are a powerful tool for modeling many
real-world phenomena (Evans, 1998, Courant and Hilbert, 1989). When derived
from the first principles, partial differential equations can serve as
predictive models which do not require any data for tuning. When learned from
data, they often outperform other models by incorporating the inductive bias
of the continuity of the modeled domain (time or space) (Chen et al., 2018,
Rubanova et al., 2019, Iakovlev et al., 2021). Inference in this type of
models is done by solving partial differential equations, that is by finding a
trajectory that satisfies the model equations and a set of initial and
boundary conditions. Since analytic solutions exist only for a limited number
of models (most likely derived from the first principles), inference is
typically done by numerical solvers of differential equations.

One way of solving differential equations is to approximate the solution by a
neural network which is trained to satisfy a given set of differential
equations, initial and boundary conditions. This approach is known in the
literature under the name of _physics-informed neural networks_ (PINNs,
Lagaris et al., 1998, Raissi et al., 2019) and it can be seen as a machine
learning alternative to classical numerical solvers. Despite the conceptual
simplicity and elegance of the method, training PINNs is notoriously difficult
in practice (Wang et al., 2021a, Krishnapriyan et al., 2021). It requires
balancing of multiple terms in the loss function (Wang et al., 2021a; c) and
the commonly used neural network architectures and parameter initialization
schemes may not work best for PINNs (Wang et al., 2021b, Sitzmann et al.,
2020).

Early stage

Intermediate results

Final results

Figure 1: The training procedure for the convection equation (Eqs. 9– 10) when
the solution is known for $t=0$ and $t=2$. Gradual expansion of the solution
domain (represented by the blue and black points) from both ends of the
interval is possible due to the flexibility of the proposed approach. The dot
color scheme is from Fig. 4.

It has been noted by many practitioners that training PINNs often results in
convergence to bad solutions (see, e.g., Krishnapriyan et al., 2021, Sitzmann
et al., 2020, Wang et al., 2022). Several recent works (Wight and Zhao, 2020,
Krishnapriyan et al., 2021, Mattey and Ghosh, 2022) address this problem by
splitting the time interval into sub-intervals and sequentially training PINNs
on each sub-interval. This idea is often referred in the literature as time
adaptation or time marching. The time-adaptive strategies make the training
procedure of PINNs similar to classical numerical solvers which compute the
solution gradually moving from the initial conditions towards the other end of
the time interval. Similarly to classical solvers, many existing PINN
algorithms are based on a pre-defined schedule of time adaptation.

Figure 2: The violin plots of the relative error $r$ for the proposed method
(blue violins on the right) and strong PINN baselines in solving three systems
(Eqs. 9–14). The proposed method yields stable training.

In this work, we follow the idea of the gradual expansion of the solution
interval when training PINNs. We propose to train an ensemble of PINNs and use
the ensemble agreement (confidence) as the criterion for expanding the
solution interval to new areas. Due to the flexibility of the domain
expansion, in contrast to existing baselines, our algorithm can easily
incorporate measurements in arbitrary locations: the algorithm will propagate
the solution from all the locations where supervision targets are known (see
Fig. 1 as an example). The proposed algorithm does not need a pre-defined
schedule of interval expansion and it treats time and space equally. We
experimentally show that the proposed algorithm can stabilize PINN training
(see Fig. 2) and yield performance competitive to the recent variants of PINNs
which use time adaptation.

##  2 Background

###  2.1 Physics-informed neural networks

Physics-informed neural networks are neural networks which are trained to
approximate the solution of a partial differential equation

$\frac{\partial{u}(x,t)}{\partial
t}=f\left(\frac{\partial^{2}{u}(x,t)}{\partial{x}^{2}},\frac{\partial{u}(x,t)}{\partial
x},x,t\right)$ (1)

on the interval $x\in[X_{1},X_{2}],t\in[T_{1},T_{2}]$ with initial conditions

${u}(x,T_{1})={u}_{0}(x),\quad x\in[X_{1},X_{2}]$ (2)

and boundary conditions

${\bm{g}}({u}(x,t))=0,\quad x\in\\{X_{1},X_{2}\\}\,.$ (3)

Functions $f$ and ${u}_{0}$ in Eqs. 1–2 are assumed to be known, function
${\bm{g}}$ in Eq. 3 is known as well and it can represent different types of
boundary conditions (e.g., Neumann, Dirichlet, Robin or periodic boundary
conditions). In the PINN approach, one approximates the solution with a neural
network which takes inputs $x$ and $t$ and produces ${u}(x,t)$ as the output.
The network is trained by minimizing a sum of several terms:

$\mathcal{L}=\mathcal{L}_{\text{S}}+\mathcal{L}_{\text{B}}+\mathcal{L}_{\text{PDE}}\text{.
}$ (4)

$\mathcal{L}_{\text{S}}$ is the standard supervised learning loss which makes
the neural network fit the initial conditions:

$\displaystyle\mathcal{L}_{\text{S}}=\sum_{i=1}^{K_{\text{S}}}w_{i}\left({u}(x_{i},t_{i})-{u}_{i}\right)^{2}\,.$
(5)

where $w_{i}$ are point-specific weights, $t_{i}=T_{1}$,
${u}_{i}={u}_{0}(x_{i})$ and $x_{i}$ are sampled from $[X_{1},X_{2}]$.
$\mathcal{L}_{\text{B}}$ is the loss computed to satisfy the boundary
conditions:

$\displaystyle\mathcal{L}_{\text{B}}=\sum_{j=1}^{K_{\text{B}}}w_{j}||{\bm{g}}({u}(x_{j},t_{j}))||^{2}\,,$
(6)

where $w_{j}$ are point-specific weights, $x_{j}\in\\{X_{1},X_{2}\\}$ and
$t_{j}$ are sampled from $[T_{1},T_{2}]$. $\mathcal{L}_{\text{PDE}}$ is the
loss derived from the PDE in 1:

$\displaystyle\mathcal{L}_{\text{PDE}}=\sum_{k=1}^{K_{\text{PDE}}}w_{k}\left(\frac{\partial{u}_{k}}{\partial
t}-f\left(\frac{\partial^{2}{u}_{k}}{\partial{x}^{2}},\frac{\partial{u}_{k}}{\partial
x},x_{k},t_{k}\right)\right)^{2}\,,$ (7)

where $w_{k}$ are point-specific weights and the partial derivatives
$\frac{d{u}_{k}}{dt}$, $\frac{\partial^{2}{u}_{k}}{\partial{x}^{2}}$,
$\frac{\partial{u}_{k}}{\partial x}$ are computed at collocation points
$(x_{k},t_{k})$ sampled from the interval
$x_{k}\in[X_{1},X_{2}],t_{k}\in[T_{1},T_{2}]$. Classical PINNs use shared
weights $w_{\text{S}}=w_{i},\forall i$, $w_{\text{B}}=w_{j},\forall j$,
$w_{\text{PDE}}=w_{k},\forall k$. The values of the weights vary depending on
the implementation and in the simplest case they are
$w_{\text{S}}=1/K_{\text{S}}$, $w_{\text{B}}=1/K_{\text{B}}$,
$w_{\text{PDE}}=1/K_{\text{PDE}}$.

The method can easily be extended to fit a sequence of observations
$\\{((x^{*}_{i},t^{*}_{i}),{u}^{*}_{i})\\}_{i=1}^{N}$ by including the
observed data to the supervision loss in Eq. 5. In this case, the method can
be seen as fitting a neural network to the training data (containing the
observations) while regularizing the solution using the PDE loss in Eq. 7. One
advantage of the PINN method compared to traditional numerical solvers is the
ability to work with ill-posed problems, for example, if the initial
conditions are known only in a subset of points. More details on existing
extensions of PINNs can be found in Appendix A.1.

###  2.2 PINNs with expansion of the time interval

Recently, several papers have proposed to solve the initial-boundary value
problem by gradually expanding the interval from which points $(x_{k},t_{k})$
in Eq. 7 are sampled. These versions of PINNs are closest to our approach.

#### Time-adaptive strategies

Wight and Zhao (2020) and Krishnapriyan et al. (2021) show the effectiveness
of the time marching strategy: the time interval $[T_{1},T_{2}]$ is split into
multiple sub-intervals and the equation is solved on each individual sub-
interval sequentially by a separate PINN. The solution at the border of the
previous sub-interval is used as the initial condition for the next one.

A similar approach is based on progressive expansion of the time interval by
gradually increasing the end point $T_{2}$ during training (Wight and Zhao,
2020, Mattey and Ghosh, 2022). Backward compatible PINNs (bc-PINNs, Mattey and
Ghosh, 2022) implement this idea such that the solution found during the
previous interval extension is used as the PINN targets during training on a
newly expanded interval to prevent catastrophic forgetting.

All these time adaptation strategies need a pre-defined schedule for the time
interval expansion, which makes them similar to classical numerical solvers
which use pre-defined discretization schemes.

#### Causality training

The idea of time adaptation is closely related to the adaptive weighting
scheme that respects causality (Wang et al., 2022). The authors use adaptive
weights for the individual terms in Eq. 7 such that the weights are computed
using the cumulative PDE loss for the preceding points:

$\displaystyle\textstyle
w_{k}=\exp\left(-\sum_{k^{\prime}|t_{k^{\prime}}<t_{k}}\epsilon\mathcal{L}_{\text{PDE}}(t_{k^{\prime}})\right),$
(8)

where $\mathcal{L}_{\text{PDE}}(t_{k^{\prime}})$ denotes an individual term in
Eq. 7 that corresponds to time point $t_{k^{\prime}}$. The idea is to zero out
the effect of the points that are far away from the initial conditions until
the solution is approximated well on all the points before them. The method
assumes that the boundary conditions are enforced as hard constraints and thus
the total loss consists of $\mathcal{L}_{\text{S}}$ and
$\mathcal{L}_{\text{PDE}}$.

##  3 Ensembles of PINNs

###  3.1 Motivation

To motivate our approach, we demonstrate failure cases of PINNs using an
example from (Krishnapriyan et al., 2021) on solving a convection equation.
Fig. 3 shows the ground-truth solution of the equation (Fig. 3a) and two
inaccurate solutions (Figs. 3b, c) found by PINNs trained with the same
network architecture but different random seeds for weight initialization. The
found PINN solutions can be seen as a combination of two solutions: the
correct solution near the initial conditions (for small $t$) and a simpler
solution farther away from the initial conditions (for large $t$). The second
row in Fig. 3 illustrates that the second, simpler solution satisfies well the
solved PDEs. This example illustrates that simple solutions can be attractive
for PINNs, which can cause problems for the optimization procedure. Once a
PINN finds a simple, locally consistent solution in some areas (for example,
far away from the initial conditions), it may be very difficult to change it.
This leads to a final solution which is a combination of the correct solution
and a wrong one.

(a) ground truth solution

(b) PINN solution 1

(c) PINN solution 2

(d) PDE error for PINN solution 1

(e) PDE error for PINN solution 2

Figure 3: Example of solving the convection system (Eqs. 9–10) for $\beta=30$,
$t\in[0,1]$ using PINNs with the LFBGS optimizer. (a): The ground truth
solution. (b)-(c): Two solutions found with different initializations of
PINNs. (d)-(e): The logarithm of the individual terms in Eq. 7 as a function
of $x$ and $t$ for the solutions.

This example provides the following intuition: when points located far away
from the initial conditions contribute to the PDE loss in Eq. 7, it may hurt
the optimization procedure by pulling the solution towards a bad local
optimum. On the other hand, including those points in the PDE loss at the
beginning of training hardly brings any benefits: it makes little sense to
regularize the solution using the loss in Eq. 7 before we know its approximate
shape. The ability to escape from wrong solutions largely depends on the
design choices made for the PINN training such as the optimizer type, the use
of mini-batches, type of the sampling utilized for points in
$\mathcal{L}_{\text{PDE}}$ term in Eq. 7, normalization of the inputs, hard or
soft encoding of the initial and boundary conditions and so on. While some of
these tricks can be beneficial for the PINN accuracy on particular systems, it
is quite difficult to select a common set of settings beneficial across a
wider range of PDEs.

The illustrated problem is avoided by the classical numerical solvers because
they usually “propagate” the solution from the initial and boundary conditions
to cover the entire interval using a schedule determined by the discretization
scheme.

###  3.2 Method

In this paper, we propose to gradually expand the areas from which we sample
collocation points $(x_{k},t_{k})$ to compute loss $\mathcal{L}_{\text{PDE}}$.
Our approach is based on training an ensemble of PINNs: a set of neural
networks initialized with different weights but trained using the same loss
function. Since PINN ensemble members typically converge to the same solution
in the vicinity of observed data but may favor distinct wrong solutions
farther away from the observations (Fig. 3b-c), we can use the ensemble
agreement as the criterion for including new points for computing loss
$\mathcal{L}_{\text{PDE}}$. Optionally, if all ensemble members agree on the
solution in a particular point, we can create a pseudo-label (taken as the
median of the ensemble predictions) for that point and make this point
contribute to the supervision loss $\mathcal{L}_{\text{S}}$ in Eq. 5.

At initialization | After 5000 updates |  After 24000 updates
---|---|---
|  |
|  |
Figure 4: Illustration of the proposed algorithm on solving the convection
system with $\beta=20$ on $t\in[0,1]$. Top row: ensemble median, bottom row:
log-variance of ensemble predictions. Dark blue dots: set $D_{\text{L}}$ of
points which contribute to loss $\mathcal{L}_{\text{S}}$; light blue dots: set
$D_{\text{PL}}$ of points with pseudo-labels (contribute either to
$\mathcal{L}_{\text{S}}$ or $\mathcal{L}_{\text{PDE}}$); black dots: points
that contribute to loss $\mathcal{L}_{\text{PDE}}$; red dots: points that
contribute to loss $\mathcal{L}_{\text{B}}$.  Algorithm 1 Training PINNs with
label propagation

1:$\Delta_{\text{PDE}}$, $\Delta$, $\sigma^{2}$ and $\epsilon$

2:$D_{\text{L}}=\\{((x_{i},t_{i}),{u}_{i})\\}$ $\triangleright$ Points with
targets (blue dots)

3:$D_{\text{PL}}=\\{\\}$ $\triangleright$ Points with pseudo-labels (light
blue dots)

4:$I_{\text{B}}=\\{(x_{j},t_{j})\\}$ $\triangleright$ Sample candidate points
for computing $\mathcal{L}_{\text{B}}$ (empty red circles)

5:$I_{\text{PDE}}=\\{(x_{m},t_{m})\\}$ $\triangleright$ Sample candidate
points for computing $\mathcal{L}_{\text{PDE}}$ (empty black circles)

6:while not converge do

7: $D=D_{\text{L}}\cup D_{\text{PL}}$

8: $I_{\text{B}}^{\prime}=\\{(x_{j},t_{j})\in
I_{\text{B}}\mid\textsc{distance}((x_{j},t_{j}),D)<\Delta_{\text{PDE}}\\}$
$\triangleright$ red dots

9: $I_{\text{PDE}}^{\prime}=\\{(x_{m},t_{m})\in
I_{\text{PDE}}\mid\textsc{distance}((x_{m},t_{m}),D)<\Delta_{\text{PDE}}\\}$
$\triangleright$ black dots

10: Train $L$ networks $f_{l}$ for $N$ iterations using $D_{\text{L}}\cup
D_{\text{PL}}$ (’PL’ version) or $D_{\text{L}}$ (’Ens’ version),
$I_{\text{B}}^{\prime}$, $I_{\text{PDE}}^{\prime}$ to compute
$\mathcal{L}_{\text{S}}$, $\mathcal{L}_{\text{B}}$,
$\mathcal{L}_{\text{PDE}}$, respectively.

11: for $(x_{m},t_{m})\in I_{\text{PDE}}^{\prime}$ do

12: $\hat{{u}}_{l}=f_{l}(x_{m},t_{m}),\quad\forall l\in 1,...,L$
$\triangleright$ Prediction of each ensemble network

13: $v_{m}=\text{variance}(\hat{{u}}_{1},\dots,\hat{{u}}_{L})$
$\triangleright$ Variance of predictions

14: ${\bar{{u}}}_{m}=\text{median}(\hat{{u}}_{1},\dots,\hat{{u}}_{L})$
$\triangleright$ Median of predictions

15: if $v_{m}<\sigma^{2}\And\textsc{distance}((x_{m},t_{m}),D)<\Delta$ then

16: $D_{\text{PL}}\leftarrow D_{\text{PL}}\cup((x_{m},t_{m}),{\bar{{u}}}_{m})$
$\triangleright$ Add a point with a pseudo-label

17: end if

18: end for

19:end while

20:

21:function distance($(x,t),D$)

22: $D^{\prime}=\\{(x_{i},t_{i})\in
D\mid||\frac{1}{L}\sum_{l}f_{l}(x_{i},t_{i})-{u}_{i}||<\epsilon\\}$
$\triangleright$ Points with good fit to targets

23: return $\min_{(x_{i},t_{i})\in D^{\prime}}||(x_{i},t_{i})-(x,t)||$

24:end function

The proposed algorithm is illustrated in Fig. 4. At the beginning of training,
the supervision loss $\mathcal{L}_{\text{S}}$ is computed using points sampled
at the initial conditions (the blue dots in Fig. 4a) and losses
$\mathcal{L}_{\text{PDE}}$ and $\mathcal{L}_{\text{B}}$ are computed using
only points that are close enough to the initial conditions (the black and red
dots in Fig. 4a respectively). The proximity is measured by thresholding the
Euclidean distance to the closest point among the blue dots. After $N$
training iterations, we compute the median and the variance of the ensemble
predictions (see Fig. 4b-c). If the variance in a particular location is small
enough, we use the median of the ensemble predictions at that point as a
pseudo-label and add that point to the data set which is used to compute the
supervision loss $\mathcal{L}_{\text{S}}$ (the light blue dots in Fig. 4b).
Points for pseudo-labeling are selected among the collocation points (the
black dots in Fig. 4a). Locations which are close enough to the data points
with labels or pseudo-labels are added to the set which contributes to
$\mathcal{L}_{\text{PDE}}$ and $\mathcal{L}_{\text{B}}$ (the black and red
dots in Fig. 4b). Then, we train the ensemble of PINNs for a fixed number of
iterations and again increase the sets of inputs which are used to compute the
losses in a similar way (Fig. 4d-e). The iterations continue until all
collocation points (which are pre-sampled at the beginning of the training
procedure) are included in the loss computations.

More formally, the training procedure is presented in Algorithm 1. The dots
and the circles in the algorithm refer to Fig. 4. In the experiments, we test
two versions of the proposed algorithm:

  1. 1.

Pseudo-labels (PL): a version with pseudo-labels in which the points with a
high degree of ensemble agreement are used to compute both
$\mathcal{L}_{\text{PDE}}$ and $\mathcal{L}_{\text{S}}$;

  2. 2.

PINN Ensemble (Ens): a version without pseudo-labels, in which the points with
a high degree of ensemble agreement are used to compute
$\mathcal{L}_{\text{PDE}}$ but not $\mathcal{L}_{\text{S}}$.

Each member of the ensemble is trained to minimize the loss in Eq. 4 with
shared weights $w_{\text{B}}=1/|I_{\text{B}}|$,
$w_{\text{PDE}}=1/|I_{\text{PDE}}|$ and $w_{\text{S}}=1/|D|$ for the PINN
Ensemble version and $w_{\text{S}}=1/|D\cup I_{\text{PDE}}|$ for the version
with pseudo-labels. Other weighting strategies (e.g., Wight and Zhao, 2020,
Wang et al., 2021a; c) could be combined with our approach as well.

###  3.3 Related work

#### PINNs and their extensions

The proposed algorithm is built on the idea of the gradual expansion of the
solution interval, which makes it similar to the time-adaptive techniques
(Wight and Zhao, 2020, Krishnapriyan et al., 2021, Mattey and Ghosh, 2022) and
the adaptive weighting method that respects causality (Wang et al., 2022). The
advantage of the proposed algorithm is its greater flexibility in the way of
expanding the area covered by collocation points: instead of expanding the
time interval with a pre-defined (Wight and Zhao, 2020, Krishnapriyan et al.,
2021, Mattey and Ghosh, 2022) or an automatic (Wang et al., 2022) schedule,
our algorithm considers each collocation point individually during the
expansion and it treats time and space equally. This feature allows the
application of the algorithm to datasets with an arbitrary layout of the
points with known targets. We illustrate this by solving the convection system
(Eqs. 9– 10) on the interval $t\in[0,2]$ when the solution is known for $t=0$
and $t=2$. As illustrated in Fig. 1, the algorithm finds a reasonable schedule
for expanding the area starting from both ends of the interval.

#### Label propagation and ensembles

Training an ensemble of PINNs with pseudo-labeling is related to how label
propagation is done in semi-supervised classification tasks (see, e.g, Lee et
al., 2013, Sohn et al., 2020): when a classifier becomes confident in the
predicted class of an unlabeled example, that example is added to the labeled
set. Model ensembles (Laine and Aila, 2017) or prediction ensembles (Berthelot
et al., 2019) are often used in those tasks to generate better targets. Since
PINNs are trained on real-valued targets, one can view PINNs as regression
models regularized by the loss in Eq. 7. Label propagation in regression tasks
is less studied with only a few existing works on semi-supervised regression
(Jean et al., 2018, Li et al., 2017). In our algorithm, we use the confidence
of the ensemble predictions to decide whether the solution interval can be
extended and which points can be assigned pseudo-labels.

##  4 Experiments

We test the proposed algorithm on finding the solutions of the following
differential equations:

  * •

convection equation used to model transport phenomena

$\displaystyle\frac{\partial u}{\partial t}=$
$\displaystyle-\beta\frac{\partial u}{\partial x}\text{,}\quad
x\in[0,2\pi]\text{,}\quad t\in[0,1]\text{,}\quad\beta=\text{const}$ (9)
$\displaystyle u(x,0)=$ $\displaystyle\sin(x)\text{,}\quad u(0,t)=u(2\pi,t)$
(10)

  * •

reaction system for modelling chemical reactions

$\displaystyle\frac{\partial u}{\partial t}=$ $\displaystyle\rho
u(1-u)\text{,}\quad x\in[0,2\pi]\text{,}\quad
t\in[0,1]\text{,}\quad\rho=\text{const}$ (11) $\displaystyle u(x,0)=$
$\displaystyle\exp\left(-{8(x-\pi)^{2}}/{\pi^{2}}\right)\text{,}\quad
u(0,t)=u(2\pi,t)$ (12)

  * •

reaction-diffusion equation that models reactions together with diffusion of
substances

$\displaystyle\frac{\partial u}{\partial t}=$
$\displaystyle\nu\frac{\partial^{2}u}{\partial{x}^{2}}+\rho
u(1-u)\text{,}\quad x\in[0,2\pi]\text{,}\quad t\in[0,1]\text{,}\quad\nu\text{,
}\rho=\text{const}\text{,}$ (13) $\displaystyle u(x,0)=$
$\displaystyle\exp\left(-{8(x-\pi)^{2}}/{\pi^{2}}\right)\text{,}\quad
u(0,t)=u(2\pi,t)\text{,}\quad u_{x}\left(0,t\right)=u_{x}\left(2\pi,t\right)$
(14)

  * •

diffusion equation with periodic boundary conditions

$\displaystyle\frac{\partial u}{\partial t}$
$\displaystyle=\frac{1}{d^{2}}\frac{\partial^{2}u}{\partial{x}^{2}}\text{,}\quad
x\in[0,2\pi]\text{, }\quad t\in[0,1]\text{,}\quad d=\text{const}$ (15)
$\displaystyle u(x,0)$ $\displaystyle=\sin(dx)\quad
u(0,t)=u(2\pi,t)\text{,}\quad u_{x}\left(0,t\right)=u_{x}\left(2\pi,t\right)$
(16)

  * •

diffusion equation in 15 with boundary conditions of the Dirichlet type:

$\displaystyle u(0,t)=u(2\pi,t)=0$ (17)

Start of training

Intermediate results

Final results

(a) reaction equation (Eqs. 11–12), $\rho=8$

(b) reaction-diffusion equation (Eqs. 13–14), $\rho=5$, $\nu=5$

(c) diffusion equation with periodic boundary conditions (Eqs. 15–16), $d=10$

Figure 5: The training progress of PINN ensembles with pseudo-labels. The dot
color scheme is from Fig. 4.

In all the experiments, we use a multi-layer perceptron with four hidden
layers with 50 neurons and the tanh activation in each hidden layer as a
backbone PINN. Our ensemble of PINNs contains five such networks. The two
inputs $x$ and $t$ of the network are normalized to $[-1,1]$, which has a
positive effect on the accuracy in our experiments. The models are trained
either with the Adam optimizer (Kingma and Ba, 2014) with learning rate 0.001
or with Adam followed by fine-tuning with LBFGS (Liu and Nocedal, 1989). The
LBFGS fine-tuning is done before each update of the collocation points which
contribute to loss $\mathcal{L}_{\text{PDE}}$ and at the very end of training.

Illustrations of the training procedure for the considered systems can be
found in Figs. 4 and 5. The plots show that the proposed algorithm finds
accurate solutions for the considered systems.

Table 1:  The means (and standard deviations) across 10 random seeds of the
relative $l_{2}$ errors in Eq. 18 for the convection (Eqs. 9–10), reaction
(Eqs. 11–12) and reaction-diffusion (Eqs. 13–14) systems. All numbers should
be multiplied by $10^{-3}$. The rows marked with ∗ show the results obtained
with $K_{\text{PDE}}=25600$ collocation points, otherwise
$K_{\text{PDE}}=1000$. In

bold

: the best three results.  Method | Convection, $\beta$ | Reaction, $\rho$ | Reaction-diffusion, $\nu$
---|---|---|---
$\beta=30$ | $\beta=40$ | $\rho=5$ | $\rho=7$ | $\nu=3$ | $\nu=4$
PINN, LBFGS | 827 (39) | 887 (21) | 985 (7) | 997 (2) | 726 (346) | 799 (257)
PINN, Adam | 13.4 (3.2) | 16.2 (4.1) | 11.9 (6.2) | 887 (269) | 252 (375) | 437 (433)
PINN, Adam+LBFGS |  9.02 (2.03) |  13.1 (2.7) | 10.8 (3.5) | 967 (26) | 170 (328) | 426 (427)
SA-PINN | 925 (35) | 980 (58) |  4.48 (0.94) | 11.5 (2.3) |  5.11 (0.53) | 5.94 (1.63)
Causality | 1425 (4) | 1393 (3) | 41.4 (4.2) | 220 (23) | 700 (4) | 713 (3)
Causality∗ | 436 (496) | 1104 (228) | 6.01 (1.87) | 43.8 (8.0) | 6.51 (0.24) | 6.56 (0.17)
bc-PINN, Adam | 212 (251) | 405 (331) | 5.37 (2.69) |  8.45 (2.14) | 6.80 (2.70) | 8.94 (4.44)
bc-PINN, Adam+LBFGS | 21.2 (40.8) | 580 (301) |  1.86 (0.81) |  3.23 (2.99) | 6.65 (0.11) |  6.51 (0.13)
PL, Adam (our) | 21.3 (1.5) | 51.0 (4.8) | 21.3 (2.5) | 55.7 (4.1) | 8.35 (0.42) | 7.95 (0.47)
PL, Adam+LBFGS (our) |  5.42 (1.23) |  6.72 (1.37) |  2.41 (0.44) |  6.02 (0.85) |  6.59 (0.01) |  6.45 (0.03)
Ens, Adam (our) |  7.00 (2.74) |  12.6 (3.0) | 7.21 (1.57) | 13.4 (3.9) | 6.68 (0.37) | 6.63 (0.30)
Ens, Adam+LBFGS (our) | 8.69 (10.5) | 14.0 (11.8) | 6.83 (1.34) | 7.07 (2.14) |  6.62 (0.01) |  6.48 (0.02)
Table 2:  The means (and standard deviations) across 10 random seeds of the
relative $l_{2}$ errors in Eq. 18 for the diffusion system (Eq. 15) with
periodic and Dirichlet (Eqs. 16– 17) boundary conditions. All numbers should
be multiplied by $10^{-3}$. The rows marked with ∗ show the results obtained
with $K_{\text{PDE}}=20000$, otherwise $K_{\text{PDE}}=1000$. In

bold

: the best three results.  Method | periodic boundary conditions | Dirichlet boundary conditions
---|---|---
$d=5$ | $d=7$ | $d=10$ | $d=5$ | $d=7$ | $d=10$
PINN, LBFGS | 997 (.1) | 999 (.1) | 999 (.1) | 998 (.1) | 999 (.1) | 999 (.1)
PINN, Adam | 7.16 (5.61) | 20.1 (18.2) | 36.6 (18.1) | 10.6 (12.0) | 15.4 (14.8) | 17.9 (10.8)
PINN, Adam+LBFGS |  0.30 (0.04) |  0.60 (0.42) |  1.00 (0.27) |  0.32 (0.07) |  0.39 (0.08) |  0.54 (0.08)
SA-PINN, Adam | 5.28 (1.68) | 5.94 (1.00) | 9.52 (2.61) | 4.85 (1.38) | 6.76 (1.46) | 8.08 (2.69)
SA-PINN, Adam+LBFGS | 2.59 (1.07) | 3.10 (1.62) | 6.46 (3.29) | 2.71 (0.96) | 3.53 (1.28) | 5.68 (2.85)
bc-PINN, Adam | 16.4 (6.1) | 82.2 (176) | 163 (212) | 13.6 (6.1) | 22.0 (7.6) | 50.4 (50.7)
bc-PINN∗, Adam | 10.7 (5.1) | 18.7 (6.9) | 48.1 (38.1) | 16.9 (12.2) | 20.5 (8.8) | 29.7 (16.1)
bc-PINN, Adam+LBFGS | 4.49 (4.67) | 70.0 (190) | 108 (257) | 3.77 (2.87) | 6.14 (2.44) | 28.1 (44.2)
PL, Adam (our) | 8.26 (3.52) | 11.5 (3.4) | 27.0 (8.0) | 9.15 (3.59) | 10.3 (4.0) | 25.9 (13.5)
Ens, Adam (our) | 5.67 (3.02) | 10.1 (4.0) | 18.2 (6.7) | 8.07 (4.95) | 10.1 (5.2) | 14.5 (7.3)
PL, Adam+LBFGS (our) |  0.30 (0.05) |  0.53 (0.04) |  1.64 (0.59) |  0.29 (0.06) |  0.55 (0.07) |  1.73 (0.58)
Ens, Adam+LBFGS (our) |  0.12 (0.02) |  0.41 (0.13) |  1.59 (0.43) |  0.11 (0.01) |  0.33 (0.14) |  2.03 (1.25)

To evaluate the accuracy of the proposed algorithm, we compute metrics used in
the previous works (Krishnapriyan et al., 2021, Wang et al., 2022): we report
the relative $l_{2}$ error

$\displaystyle\textstyle r=l_{2}(\hat{{\bm{u}}}-{\bm{u}})/l_{2}({\bm{u}})\,,$
(18)

where ${\bm{u}}$ is a vector of the ground-truth values in the test set,
$\hat{{\bm{u}}}$ is the corresponding PINN predictions and $l_{2}()$ denotes
$l_{2}$ norm. Our test set consists of points on a regular $256\times 100$
grid.

In Table 1, we compare the accuracy of the proposed method with several strong
baselines proposed recently in the literature. We present the means and
standard deviations of the solution errors obtained in 10 runs with different
initializations. Since most of the considered methods benefit from longer
training, we limit the number of network updates for all the comparison
methods. We use the following comparison methods: (1) Vanilla PINN trained
with the LBFGS optimizer, (2) Vanilla PINN trained with the Adam optimizer,
(3) Causality (Wang et al., 2022), (4) SA-PINN (McClenny and Braga-Neto,
2020), (5) bc-PINN (Mattey and Ghosh, 2022). More details on the methods and
the selection of hyperparameters can be found in Appendix A.2.

The results in Table 1 show that vanilla PINNs trained with LBFGS struggle to
find accurate solutions for all the considered systems. Vanilla PINNs trained
with Adam yield more accurate results but the results are unsatisfactory for
the reaction and reaction-diffusion systems. The causality weighting scheme
seems to require many more collocation points to find satisfactory solutions.
SA-PINN and bc-PINN work very well on the reaction and reaction-diffusion
systems but do not find good solutions for the convection system. The proposed
ensemble method provides stable training (see the summary of the results in
Fig. 2) and shows competitive performance in all the considered systems. We
can also observe that the LBFGS fine-tuning has positive effect on the
accuracy when a good approximation is found during pre-training with Adam.

In Table 2, we compare the accuracy of the proposed algorithm on solving the
diffusion system with the two types of boundary conditions: periodic (Eq. 16)
and Dirichlet-type (Eq. 17) boundary conditions. In this comparison, we omit
the causality-motivated baseline (Wang et al., 2022) because the periodic
boundary conditions in the existing implementations are enforced as hard
constraints. The results in Table 2 show that training ensembles of PINNs
yields competitive performance in these settings as well.

In Appendix A.3, we study the sensitivity of the ensemble training of PINNs to
its hyperparameters. The results show that the PINN Ensemble algorithm is
generally stable at solving the considered systems with little sensitivity to
the hyperparameter values. We note the importance of hyperparameter
$\Delta_{\text{PDE}}$ which determines how far the points considered for
inclusion can be from the points already included in the loss calculations.
The performance of the algorithm drops when $\Delta_{\text{PDE}}$ is too
large. This happens because ensemble members may be attracted by the same
trivial solution in areas too far away from the initial conditions, the effect
caused by commonly used network architectures and initialization schemes (see,
e.g., Wong et al., 2022, Rohrhofer et al., 2022), which may negatively affect
the ensemble diversity.

##  5 Discussion and future work

In this paper, we propose to stabilize training of PINNs by gradual expansion
of the solution interval based on the agreement of an ensemble of PINNs. The
obtained results suggest that the proposed approach can reduce the number of
failure cases during PINN training. Another potential advantage of the
proposed algorithm is that the PINN ensemble produces confidence intervals
which can be viewed as uncertainty estimates of the found solution (see Fig.
4c, e). Although the proposed algorithm is more computationally expensive
compared to vanilla PINNs, ensemble training can be effectively parallelized
in which case the wall clock time of training does not grow significantly. The
method shows good results for simple systems such as convection and reaction-
diffusion.

This work can be extended in a number of ways. One potential direction is to
use different ways of creating model ensembles, for example, by dropout
(Hinton et al., 2012) or pseudo-ensembles (Bachman et al., 2014). It is
interesting to investigate how the proposed algorithm can be combined with
other tricks from the PINN literature, for example, adaptive balancing of the
loss terms (Wang et al., 2021c). Another direction is to find alternative ways
of the solution interval expansion (e.g., update sets
$\mathcal{I}_{\text{PDE}}^{\prime}$, $\mathcal{I}_{\text{B}}^{\prime}$ more
frequently), which may increase the convergence speed. It should also be
possible to improve the PINN architecture such that the knowledge of a found
(local) solution in one area could be used in finding the solution in another
area. Using neural networks with the right inductive bias (e.g., similar to
the ones proposed by Sanchez-Gonzalez et al., 2020, Iakovlev et al., 2021,
Brandstetter et al., 2022), might provide a solution for that.

## Acknowledgment

We thank CSC (IT Center for Science, Finland) for computational resources and
the Academy of Finland for the support within the Flagship programme: Finnish
Center for Artificial Intelligence (FCAI).

## References

  * Bachman et al. (2014) P. Bachman, O. Alsharif, and D. Precup.  Learning with pseudo-ensembles.  _Advances in Neural Information Processing Systems_ , 27, 2014.
  * Berthelot et al. (2019) D. Berthelot, N. Carlini, I. Goodfellow, N. Papernot, A. Oliver, and C. A. Raffel.  Mixmatch: A holistic approach to semi-supervised learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Brandstetter et al. (2022) J. Brandstetter, D. E. Worrall, and M. Welling.  Message passing neural PDE solvers.  _International Conference on Learning Representations_ , 2022.
  * Chen et al. (2018) R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Courant and Hilbert (1989) R. Courant and D. Hilbert.  _Methods of mathematical physics: partial differential equations_.  John Wiley & Sons, Ltd, 1989.  ISBN 9783527617234.
  * Daw et al. (2022) A. Daw, J. Bu, S. Wang, P. Perdikaris, and A. Karpatne.  Rethinking the importance of sampling in physics-informed neural networks.  _arXiv preprint arXiv:2207.02338_ , 2022.
  * Dong and Ni (2021) S. Dong and N. Ni.  A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks.  _Journal of Computational Physics_ , 435:110242, 2021.
  * Evans (1998) L. C. Evans.  _Partial Differential Equations_ , volume 19 of _Graduate studies in mathematics_.  American Mathematical Society, 1998.  ISBN 9780821807729.
  * Hennigh et al. (2021) O. Hennigh, S. Narasimhan, M. A. Nabian, A. Subramaniam, K. Tangsali, Z. Fang, M. Rietmann, W. Byeon, and S. Choudhry.  NVIDIA SimNet: An AI-accelerated multi-physics simulation framework.  In _International Conference on Computational Science_ , pages 447–461. Springer, 2021.
  * Hinton et al. (2012) G. E. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors.  _arXiv preprint arXiv:1207.0580_ , 2012.
  * Iakovlev et al. (2021) V. Iakovlev, M. Heinonen, and H. Lähdesmäki.  Learning continuous-time PDEs from sparse data with graph neural networks.  _International Conference on Learning Representations_ , 2021.
  * Jean et al. (2018) N. Jean, S. M. Xie, and S. Ermon.  Semi-supervised deep kernel learning: Regression with unlabeled data by minimizing predictive variance.  _Neural Information Processing Systems_ , 2018.
  * Kingma and Ba (2014) D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Krishnapriyan et al. (2021) A. S. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, and M. W. Mahoney.  Characterizing possible failure modes in physics-informed neural networks.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Lagaris et al. (1998) I. E. Lagaris, A. Likas, and D. I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  _IEEE transactions on neural networks_ , 9(5):987–1000, 1998.
  * Laine and Aila (2017) S. Laine and T. Aila.  Temporal ensembling for semi-supervised learning.  _International Conference on Learning Representations_ , 2017.
  * Lee et al. (2013) D.-H. Lee et al.  Pseudo-label: The simple and efficient semi-supervised learning method for deep neural networks.  In _Workshop on challenges in representation learning, ICML_ , volume 3, page 896, 2013.
  * Li et al. (2017) Y.-F. Li, H.-W. Zha, and Z.-H. Zhou.  Learning safe prediction for semi-supervised regression.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 31(1), Feb. 2017.
  * Liu and Wang (2021) D. Liu and Y. Wang.  A dual-dimer method for training physics-constrained neural networks with minimax architecture.  _Neural Networks_ , 136:112–125, 2021.
  * Liu and Nocedal (1989) D. C. Liu and J. Nocedal.  On the limited memory BFGS method for large scale optimization.  _Mathematical programming_ , 45(1):503–528, 1989\.
  * Lu et al. (2021a) L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  DeepXDE: A deep learning library for solving differential equations.  _SIAM Review_ , 63(1):208–228, 2021a.
  * Lu et al. (2021b) L. Lu, R. Pestourie, W. Yao, Z. Wang, F. Verdugo, and S. G. Johnson.  Physics-informed neural networks with hard constraints for inverse design.  _SIAM Journal on Scientific Computing_ , 43(6):B1105–B1132, 2021b.
  * Markidis (2021) S. Markidis.  The old and the new: Can physics-informed deep-learning replace traditional linear solvers?  _Frontiers in big Data_ , page 92, 2021.
  * Mattey and Ghosh (2022) R. Mattey and S. Ghosh.  A novel sequential method to train physics informed neural networks for Allen-Cahn and Cahn-Hilliard equations.  _Computer Methods in Applied Mechanics and Engineering_ , 390:114474, 2022.  ISSN 0045-7825.
  * McClenny and Braga-Neto (2020) L. McClenny and U. Braga-Neto.  Self-adaptive physics-informed neural networks using a soft attention mechanism.  _arXiv preprint arXiv:2009.04544_ , 2020.
  * Raissi et al. (2019) 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–707, 2019\.
  * Rohrhofer et al. (2022) F. M. Rohrhofer, S. Posch, C. Gößnitzer, and B. C. Geiger.  Understanding the difficulty of training physics-informed neural networks on dynamical systems.  _arXiv preprint arXiv:2203.13648_ , 2022.
  * Rubanova et al. (2019) Y. Rubanova, R. T. Chen, and D. K. Duvenaud.  Latent ordinary differential equations for irregularly-sampled time series.  _Advances in neural information processing systems_ , 32, 2019.
  * Sanchez-Gonzalez et al. (2020) A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia.  Learning to simulate complex physics with graph networks.  In _International Conference on Machine Learning_ , pages 8459–8468. PMLR, 2020.
  * Sitzmann et al. (2020) V. Sitzmann, J. Martel, A. Bergman, D. Lindell, and G. Wetzstein.  Implicit neural representations with periodic activation functions.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 7462–7473. Curran Associates, Inc., 2020.
  * Sohn et al. (2020) K. Sohn, D. Berthelot, N. Carlini, Z. Zhang, H. Zhang, C. A. Raffel, E. D. Cubuk, A. Kurakin, and C.-L. Li.  Fixmatch: Simplifying semi-supervised learning with consistency and confidence.  _Advances in Neural Information Processing Systems_ , 33:596–608, 2020.
  * Sukumar and Srivastava (2022) N. Sukumar and A. Srivastava.  Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks.  _Computer Methods in Applied Mechanics and Engineering_ , 389:114333, 2022.
  * Tancik et al. (2020) M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. Barron, and R. Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 7537–7547. Curran Associates, Inc., 2020.
  * Wang et al. (2021a) S. Wang, Y. Teng, and P. Perdikaris.  Understanding and mitigating gradient flow pathologies in physics-informed neural networks.  _SIAM Journal on Scientific Computing_ , 43(5):A3055–A3081, 2021a.
  * Wang et al. (2021b) S. Wang, H. Wang, and P. Perdikaris.  On the eigenvector bias of fourier feature networks: From regression to solving multi-scale PDEs with physics-informed neural networks.  _Computer Methods in Applied Mechanics and Engineering_ , 384:113938, 2021b.
  * Wang et al. (2021c) S. Wang, X. Yu, and P. Perdikaris.  When and why PINNs fail to train: A neural tangent kernel perspective.  _Journal of Computational Physics_ , page 110768, 2021c.
  * Wang et al. (2022) S. Wang, S. Sankaran, and P. Perdikaris.  Respecting causality is all you need for training physics-informed neural networks.  _arXiv preprint arXiv:2203.07404_ , 2022.
  * Wight and Zhao (2020) C. L. Wight and J. Zhao.  Solving Allen-Cahn and Cahn-Hilliard equations using the adaptive physics informed neural networks.  _arXiv preprint arXiv:2007.04542_ , 2020.
  * Wong et al. (2022) J. C. Wong, C. Ooi, A. Gupta, and Y.-S. Ong.  Learning in sinusoidal spaces with physics-informed neural networks.  _IEEE Transactions on Artificial Intelligence_ , 2022.
  * Wu et al. (2022) C. Wu, M. Zhu, Q. Tan, Y. Kartha, and L. Lu.  A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks.  _arXiv preprint arXiv:2207.10289_ , 2022.
  * Zubov et al. (2021) K. Zubov, Z. McCarthy, Y. Ma, F. Calisto, V. Pagliarino, S. Azeglio, L. Bottero, E. Luján, V. Sulzer, A. Bharambe, et al.  NeuralPDE: Automating physics-informed neural networks (PINNs) with error approximations.  _arXiv preprint arXiv:2107.09443_ , 2021.

##  Appendix A Appendix

###  A.1 PINN extensions

Despite the elegance of the PINN approach, the method is known to be prone to
failures, especially when the solution has a complex shape on the considered
interval (see, e.g., Wang et al., 2021a, Krishnapriyan et al., 2021, Wang et
al., 2021c; b; 2022, Wight and Zhao, 2020). The optimization problem solved
while training PINNs is very hard and the accuracy of the found solution is
very sensitive to the hyperparameters of PINNs (see, e.g., Markidis, 2021):
the weights $w_{\text{S}}$, $w_{\text{B}}$, $w_{\text{PDE}}$ of the loss
terms, the parameterization used in the neural network and the strategy for
sampling the points in which the loss terms are computed.

Balancing shared weights $w_{\text{S}}$, $w_{\text{B}}$, $w_{\text{PDE}}$ of
the individual loss terms in Eqs. 5–7 is a popular way to improve the accuracy
of the PINN solution. Wight and Zhao (2020) propose using larger weights
$w_{\text{S}}$ relative to $w_{\text{B}}$ and $w_{\text{PDE}}$ because the
initial conditions largely determine the shape of the solution. Several works
propose different schemes for dynamically adapting $w_{\text{S}}$,
$w_{\text{B}}$, $w_{\text{PDE}}$ during training such that the weights that
correspond to problematic loss terms get higher values. Wang et al. (2021a; c)
propose to adjust the weights either based on the magnitudes of the gradients
of the corresponding loss terms or based on the eigenvalues of the limiting
Neural Tangent Kernel. Liu and Wang (2021) propose to treat weights
$w_{\text{S}}$, $w_{\text{B}}$, $w_{\text{PDE}}$ as trainable parameters and
adjust them jointly with the PINN parameters solving a minimax optimization
problem. Self-Adaptive PINNs (SA-PINNs, McClenny and Braga-Neto, 2020)
increase _point-wise_ weights $w_{i}$, $w_{j}$ and $w_{k}$ of the loss terms
in Eqs. 5–7 during training such that the changes of the weights are
proportional to the corresponding loss terms.

The accuracy of PINNs can also be improved by using a different
parameterization for the trained neural network instead of the most standard
multilayer perceptron architecture. Several works (see, e.g., Lagaris et al.,
1998, Dong and Ni, 2021, Lu et al., 2021b, Sukumar and Srivastava, 2022)
propose neural network parameterizations which guarantee that the initial or
boundary conditions are satisfied exactly, thus eliminating terms
$\mathcal{L}_{\text{S}}$ and $\mathcal{L}_{\text{B}}$ from Eq. 4. Wang et al.
(2021b) propose to use Fourier features as the inputs of the neural network,
which is motivated by the success of Fourier feature networks (Tancik et al.,
2020) in preserving high-frequencies in the modeled solution. The method
called SIREN (Sitzmann et al., 2020) proposes to use a multilayer perceptron
with periodic activation functions and adjusts the weight initialization
schemes to work better with such activation functions. This parameterization
can also improve learning of high frequencies in the modeled solution. Wang et
al. (2021a) propose to add multiplicative connections to the standard
multilayer perceptron architecture to account for possible multiplicative
interactions between different input dimensions.

Another popular way of improving the accuracy of PINNs is to use adaptive
strategies for sampling collocation points $(x_{k},t_{k})$ to compute
$\mathcal{L}_{\text{PDE}}$ in Eq. 7. Wight and Zhao (2020) propose to sample
more collocation points in the areas with large values of
$\mathcal{L}_{\text{PDE}}$, which helps to learn solutions with fast
transitions. Daw et al. (2022) propose an evolutionary strategy for sampling
the collocation points: the points with large contribution to
$\mathcal{L}_{\text{PDE}}$ are kept for the next iteration while the rest of
the points are re-sampled uniformly from the domain. These approaches resemble
the strategy of the classical solvers to reduce the discretization interval
when the solution cannot be estimated accurately. More details on the sampling
strategies for PINN and empirical comparison can be found in (Wu et al.,
2022).

Many of the improvements proposed to PINNs can be combined, which is supported
by existing software libraries (Lu et al., 2021a, Zubov et al., 2021, Hennigh
et al., 2021).

###  A.2 Training details

For all comparison methods (similarly to Krishnapriyan et al. (2021)), we use
$K_{\text{PDE}}=1000$ collocation points randomly sampled from a regular
$256\times 100$ grid on the solution interval. We use $K_{\text{S}}=256$
points to fit the initial conditions, these points are selected from a regular
grid on $x\in[0,2\pi]$ with $t=0$. We use $K_{\text{B}}=100$ for the boundary
condition loss $\mathcal{L}_{\text{BC}}$. These points are selected form a
regular grid on $t\in[0,1]$.

For the proposed ensemble method, we use the following hyperparameters:
$\sigma^{2}=0.0004$, $\epsilon=0.001$, $\Delta=0.05$ and
$\Delta_{\text{PDE}}=0.1$. We perform $N=1000$ gradient updates before
extending sets $I_{\text{PL}}^{\prime}$, $I_{\text{PDE}}^{\prime}$,
$I_{\text{B}}$, except for the first set extension which is done after
$N_{1}=5000$ training steps. We also set weighting coefficient
$w_{s}=64/|D\cup I_{\text{PDE}}|$ for runs with pseudo-labels trained with
Adam and LBFGS similar to the closest baseline.

For the baseline methods we adapted the authors’ implementations with the
following adjustments.

  * •

Causality (Wang et al., 2022): PINN with an adaptive weighting scheme that
respects causality. We report results for the values of hyperparameter
$\epsilon$ in Eq. 8 that worked best for the considered systems:
$\epsilon=0.01$ for the convection system, $\epsilon=100$ for the reaction
system and $\epsilon=1$ for the reaction-diffusion systems. We train the model
until all the adapted weights become greater than 0.99 (the stopping criterion
used by Wang et al. (2022)). As the existing implementation requires a regular
grid, we use a grid of $K_{\text{PDE}}=32\times 32=1024$ collocation points to
compute loss $\mathcal{L}_{\text{PDE}}$. We also report results with a denser
grid of $K_{\text{PDE}}=256\times 100=25600$ points for this algorithm. We do
not add normalization for network inputs $x$ and $t$ to $[-1,1]$ and follow
input encoding of the original implementation.

  * •

SA-PINN (McClenny and Braga-Neto, 2020) with trainable point-specific weights.
We switch off finetuning with the LBFGS optimizer after Adam due to observed
training instabilities. We also add normalization of network inputs to
$[-1,1]$ similar to our method.

  * •

bc-PINN (Mattey and Ghosh, 2022): a technique of expanding the time interval
with a pre-defined schedule. We expand the time interval four times and sample
250 collocation points on each interval to compute loss
$\mathcal{L}_{\text{PDE}}$ and 25 points per interval to compute
$\mathcal{L}_{\text{B}}$. We report the results obtained with the Adam
optimizer and with the Adam optimizer followed by LBFGS on each time interval.

For all of the methods, we train convection equation (Eqs. 9–10) with
$\beta=30$ for 104000 gradient updates and with $\beta=40$ for 154000 updates,
reaction (Eqs. 11–12) and reaction-diffusion (Eqs. 13–14) equations are
trained for 64000 updates, diffusion (Eqs. 15–17) equation with $d=5$ and
$d=7$ for 84000 updates and with $d=10$ for 104000 (unless an automatic
stopping criterion is used).

###  A.3 Hyperparameter sensitivity

In Table 3, we test the sensitivity of the ensemble training of PINNs to its
hyperparameters. In these experiments, we use the number of updates as
described in Appendix A.2. We consider that a run has not converged if less
than 95% of the sampled points have been added to the set which is used to
compute $\mathcal{L}_{\text{PDE}}$ by the last update. Such runs are excluded
from the statistics reported in Table 3. This step is done only in the
hyperparameter sensitivity experiments.

The most important hyperparameters are the variance threshold $\sigma^{2}$ and
the distance parameters $\Delta$ and $\Delta_{\text{PDE}}$ which determine how
quickly the algorithm expands the solution interval. When $\sigma^{2}$ is
small then the interval is expanded more slowly and the algorithm may require
more iterations to converge. Too large values of $\sigma^{2}$ can result in
adding new regions in the training procedure too early. The distance
hyperparameters $\Delta$ and $\Delta_{\text{PDE}}$ have a similar effect.
However, for the considered systems, the model performance is stable with
little sensitivity to the values of these hyperparameters as well as to the
number $K_{\text{PDE}}$ of collocation points and the number of networks in
the ensemble.

Table 3:  The results obtained with the proposed ’PINN Ensemble’ algorithm
using different hyperparameter values. We report the same metric for the same
systems as in Table 1. The superscript (n) indicates how many runs among the
10 runs did not converge after a fixed number of updates. We assume that the
algorithm has converged if more than 95% of the sampled points have been added
to the set used to compute loss $\mathcal{L}_{\text{PDE}}$.

| Convection, $\beta$ | Reaction, $\rho$ | Reaction-diffusion, $\nu$
---|---|---|---
$\beta=30$ | $\beta=40$ | $\rho=5$ | $\rho=6$ | $\rho=7$ | $\nu=2$ | $\nu=3$ | $\nu=4$
default | 7.00 (2.7) | 12.6 (3) | 7.21 (1.6) | 11.0 (2) | 13.4 (4) | 6.93 (.3) | 6.68 (.4) | 6.63 (.3)
Varying the threshold for the ensemble disagreement, default
$\sigma^{2}=0.0004$
$\sigma^{2}=10^{-4}$ | 8.15 (2.7)(1) | 14.0 (2) | 7.16 (1.6) | 11.3 (2) | 12.3 (3)(3) | 6.67 (.4) | 6.73 (.6) | 6.60 (.5)
$\sigma^{2}=10^{-3}$ | 8.11 (2.0) | 14.5 (4) | 7.16 (1.5) | 11.3 (2) | 13.0 (4) | 6.68 (.3) | 6.79 (.4) | 6.63 (.4)
Varying the distances which determinate candidate collocation points, default
$\Delta_{\text{PDE}}=0.1$, $\Delta=0.05$
$\Delta_{\text{PDE}}=.07$ | 10.1 (4.4) | 17.6 (10)(2) | 6.80 (1.1) | 11.5 (3) | 10.7 (3)(2) | 7.03 (.3) | 6.76 (.6) | 6.68 (.6)
$\Delta_{\text{PDE}}=.2$ | 5.33 (0.9) | 11.5 (2) | 7.61 (2.) | 12.7 (3) | 11.6 (4) | 7.13 (.4) | 6.89 (.4) | 6.49 (.4)
$\Delta=.07$ | 7.82 (1.6) | 13.2 (2) | 7.67 (1.5) | 12.0 (2) | 15.3 (4) | 6.70 (.6) | 6.64 (.5) | 6.63 (.4)
| $\Delta_{\text{PDE}}=.25$
---
& $\Delta=.125$
6.25 (1.6) | 13.0 (2) | 8.42 (2.) | 14.2 (3) | 14.5 (8) | 6.97 (.4) | 6.88 (.6) | 6.75 (.4)
$\Delta=10^{5}$ | 8.98 (1.9) | 13.7 (2) | 8.33 (1.6) | 13.6 (2) | 18.4 (3) | 6.69 (.2) | 6.91 (.6) | 6.62 (.4)
$\Delta_{\text{PDE}}=10^{5}$ | 11.1 (1.2) | 16.9 (3) | 8.22 (1.5) | 56 (71) | 36.3 (53) | 455 (143) | 242 (164) | 74 (82)
Varying the prediction error which determinate candidate collocation points,
default $\epsilon=0.001$
$\epsilon=5\cdot 10^{-4}$ | 8.94 (2.3) | 13.9 (6) | 7.13 (1.5) | 11.2 (2) | 12.9 (4)(1) | 6.66 (.3) | 6.93 (.7) | 6.61 (.4)
$\epsilon=10^{-2}$ | 7.65 (1.7) | 12.3 (2) | 7.15 (1.5) | 11.4 (2) | 13.0 (4)(1) | 6.72 (.4) | 6.93 (.7) | 6.63 (.3)
$\epsilon=10^{5}$ | 7.19 (1.4) | 12.2 (2) | 7.11 (1.5) | 11.4 (2) | 12.9 (4)(1) | 6.72 (.4) | 6.93 (.7) | 6.61 (.4)
Varying the number of the collocation points, default $K_{\text{PDE}}=1000$
$K_{\text{PDE}}=5000$ | 8.02 (1.9) | 14.1 (2) | 6.75 (1.4) | 10.7 (2) | 10.4 (3) | 6.88 (.4) | 6.52 (.5) | 6.55 (.5)
$K_{\text{PDE}}=10^{4}$ | 7.71 (1.7) | 16.8 (4) | 7.00 (1.3) | 9.89 (2) | 10.3 (4) | 6.94 (.4) | 6.71 (.5) | 6.55 (.4)
Number of training epochs for the first iteration, default $N_{1}=5000$
$N_{1}=10^{3}$ | 7.12 (1.0) | 15.6 (3) | 7.86 (1.3) | 12.6 (2) | 14.8 (2)(2) | 6.89 (.3) | 6.45 (.4) | 6.58 (.2)
$N_{1}=10^{4}$ | 10.8 (7) | 15.0 (8)(1) | 6.26 (1.4) | 9.18 (2) | 9.73 (3.2) | 6.79 (.3) | 6.67 (.4) | 6.66 (.5)
Number of training epochs between adding new points, default $N=1000$
$N=700$ | 7.64 (1.4) | 15.1 (7) | 7.53 (1.7) | 12.1 (2) | 14.3 (3)(1) | 6.91 (.4) | 6.66 (.4) | 6.65 (.2)
$N=2000$ | 9.95 (1.7) | 14.4 (5) | 6.62 (1.3) | 9.99 (2) | 11.8 (3) | 6.96 (.4) | 6.87 (.4) | 6.54 (.5)
Number of networks, default $L=5$
$L=3$ | 8.96 (3.7) | 12.2 (4) | 6.45 (1.4) | 11.6 (4) | 14.4 (6)(3) | 6.83 (.3) | 6.60 (.3) | 6.75 (.6)
$L=7$ | 6.68 (1.6) | 12.3 (2) | 6.07 (.81) | 11.4 (2) | 11.9 (3)(5) | 6.88 (.2) | 6.78 (.4) | 6.57 (.2)
$L=10$ | 7.94 (1.4) | 13.6 (2) | 6.48 (.58) | 10.6 (.9) | 12.5 (3) | 7.05 (.3) | 6.72 (.1) | 6.59 (.2)

###  A.4 More experimental results

In this section we report additional experimental results for reaction (Eqs.
11–12) with $\rho=6$ and reaction-diffusion (Eqs. 13–14) with $\nu=2$.

Table 4:  The means (and standard deviations) across 10 random seeds of the
relative $l_{2}$ errors in Eq. 18 for reaction (Eqs. 11–12) and reaction-
diffusion (Eqs. 13–14) systems. All numbers should be multiplied by $10^{-3}$.
The rows marked with ∗ show the results obtained with $K_{\text{PDE}}=25600$
collocation points otherwise $K_{\text{PDE}}=1000$. In

bold

: the best three results.

Method | Reaction, $\rho$ | Reaction-diffusion, $\nu$
---|---|---
$\rho=6$ | $\nu=2$
PINN, LBFGS | 993 (2) | 608 (469)
PINN, Adam | 592 (459) | 7.01 (0.74)
PINN, Adam+LBFGS | 678 (433) |  6.83 (0.03)
SA-PINN | 8.12 (1.40) | 6.97 (0.54)
Causality | 102 (35) | 678 (5)
Causality∗ | 18.9 (4.1) | 7.06 (0.27)
bc-PINN, Adam |  5.66 (2.02) | 7.11 (1.35)
bc-PINN, Adam+LBFGS |  2.96 (1.25) | 6.89 (0.08)
PL, Adam (our) | 38.5 (3.5) | 8.83 (0.42)
PL, Adam+LBFGS (our) |  4.17 (0.58) |  6.83 (0.02)
Ens, Adam (our) | 11.0 (2.3) | 6.93 (0.32)
Ens, Adam+LBFGS (our) | 9.73 (2.02) |  6.85 (0.02)