Documents Live, a web authoring and publishing system
If you see this, something is wrong
Table of contents
First published on Sunday, Jun 7, 2026 and last modified on Wednesday, Jun 10, 2026 by François Chaplais.
Graduate School of Information Science and Technology Email
Graduate School of Information Science and Technology and RIKEN Center for Advanced Intelligence Project (AIP) Email
Abstract
Physics-Informed Neural Networks (PINNs) are a common class of machine learning-based partial differential equation (PDE) solvers which train a network to represent a solution by minimizing a residual loss that encodes the PDE. Despite their successes, they are known to fail on certain simple equations, converging to an incorrect solution despite low loss. These failure modes have garnered significant attention in the literature over the past several years, motivating both architectural and optimization based solutions. By directly visualizing the residual, we show that failure modes are the result of overfitting: the loss is minimized on the collocation points, but not elsewhere. Applying regularization causes the failure modes to vanish. Finally, we extend double backpropagation over the full set of residuals, and use it to achieve state-of-the-art performance on four standard failure mode equations with up to \( 23\times\) fewer collocation points and a vanilla architecture.
1 Introduction
Using neural networks to represent the solution of partial differential equations has garnered much attention over the past several years. Many different approaches have emerged, from methods using networks as one component of the solution [1, 2], to those based on classical methods like the Ritz and Galerkin methods [3, 4]. Other approaches have found success in approximating the solution operators instead of the solutions themselves [5, 6, 7, 8, 9, 10]. However, the formalism of Physics-Informed Neural Networks (PINNs) [11] remains one of the most studied. PINNs use a neural network to directly approximate the solution of the partial differential equation, and train it by minimizing the residual of the PDE operator applied to the network in tandem with losses representing the boundary and initial conditions. The architecture is simple to implement, and has shown great success.
Yet despite this, PINNs can struggle on relatively simple problems. Most notably is the convection / linear advection equation, which acts to translate the initial condition. While PINNs learn the solution easily for weak convection speeds, they collapse for strong convection, learning instead the trivial zero solution over most of the domain. This convergence to the incorrect solution is commonly referred to as a failure mode [12].
Failure modes have proven an immense challenge for PINNs, and various diagnoses and fixes have been proposed. Reasons for failure have been broadly attributed to spectral bias and loss imbalances between the domain and boundary terms [13, 14], propagation issues from the initial condition into the domain [15, 16, 17], as well as the complexity of the loss landscape [18, 19, 20]. Proposed solutions include loss reweighting, adaptive sampling techniques, novel optimization methods, new architectures, and notably, a simple upgrade from floating point 32 (FP32) to floating point 64 (FP64) [13, 14, 16, 17, 18, 19, 20].
While the focus on the optimization difficulties has delivered many useful architectures, and improved our understanding of PINNs and machine learning optimization more generally, even the state-of-the-art method (FP64 optimization) [20], still struggles with a relatively high error, leading to a suboptimal performance when compared to traditional methods. Interestingly, as we will show, the high-precision method can still suffer from failure modes.
If a robust optimization procedure cannot fully eliminate the collapse of PINNs into failure modes, then what is the underlying cause inhibiting convergence to the true solution? We argue that there are two fundamental issues, that have been conflated under the broad umbrella of optimization difficulties. The first is the genuine optimization problems that have been diagnosed in detail [18, 19, 20]. These inhibit convergence in higher order problems, such as the wave or Allen-Cahn equations. Previous literature on failure modes has mostly targeted this perspective.
The second, which has remained undiagnosed, is overfitting. By directly visualizing the PDE residual of a failure mode, we show that in regions of high residual the loss has been minimized only at the collocation points, leading to a system well converged to the wrong solution. This diagnosis is consistent across equations, and is further confirmed by a divergence between a training and test loss when the system enters a failure mode.
We then demonstrate that the success of increased precision relies on using a dense collocation grid. Reducing the number of points causes failure even at higher precision, implying that FP64 training in part relies on over-determining the problem. This too is consistent with overfitting: more collocation points makes it harder for the network to memorize the data. Yet increasing the points directly increases the cost of training.
By incorporating standard \( L^2\) regularization used for reducing overfitting in broader machine learning, we show the failure modes vanish even at FP32, while also allowing us to use fewer collocation points over the domain (\( <1000\) , compared to \( >10000\) in Refs [21, 22, 20]). Then, we generalize double backpropagation [23] to PINNs, applying it to each of the residual terms, which both improves training efficiency and final accuracy.
Combining this double-backpropagation gradient regularization with FP64 optimization, then allows us to outperform state-of-the-art by over an order of magnitude in error on the full benchmark suite of failure mode equations, while using far fewer points.
1.1 Contributions.
Our main contributions are as follows:
- We reveal that while PINNs failure modes are considered to be caused by the optimization landscape, they are in fact caused by overfitting, where the optimization minimizes the loss only at the collocation points in regions of high residual. This changes the understanding of why and when PINNs fail.
- Using this understanding, we show that adding a standard \( L^2\) regularizer takes the optimization process from failure to success, even at FP32.
- Finally, we extend double backpropagation (gradient regularization) to PINNs, which fits in naturally into the PINNs framework, mitigates failure modes, improves optimization efficiency and reduces collocation requirements. This allows us to achieve state-of-the-art performance on the standard set of failure mode equations with a vanilla architecture.
2 Related Works
2.1 Physics Informed Neural Networks and Their Failure Modes
A large body of literature has focused on applying machine learning methods to the solution of PDEs [24, 25, 4, 3, 26, 2]. PINNs are one of the most popular formulations, showing how the methods can be applied to both forward and inverse problems, as well as the data-driven discovery of PDEs themselves [11]. Our focus here is on the forward problem, meaning using the neural network itself as an approximation to the solution. This is in contrast to the distinct but related field of neural operators, which learn a solution functional [5, 6, 8, 7, 10].
Despite their promise and success in high-dimensional problems, PINNs have been known to exhibit optimization difficulties, collapsing and failing to learn the solution for relatively simple low dimensional problems. These have been termed failure modes, and result in a seemingly successful optimization, with a loss near zero, yet having learned a solution far from the desired result [12, 13]. These failure modes have been extensively studied, and several techniques and architectures have been proposed with the goal of ameliorating them [12, 13, 27, 28, 17, 22, 16, 19, 29, 30, 15, 31, 32]. Notably, it has been recently proposed that the optimization difficulties stem from insufficient numerical precision inhibiting optimization, rather than from successful optimization to an incorrect solution [20]. We will show that while optimization can be improved with precision, the origin of the failure modes lies in overfitting, not optimization difficulty.
2.2 Regularization
Regularization is a standard tool in deep learning [33, 34]. Its goal is to improve the generalizability of the machine learning model by reducing overfitting to the training dataset. Common versions include applying a soft constraint on the weights in the form of an addition to the loss, such as with \( L^2\) and \( L^1\) style regularizations [35, 36], and derivative methods like double backpropagation and Sobolev training [23, 37]. Regularization has rarely been applied to PINNs. Occasionally, the PINN losses have been referred to as regularizers themselves, particularly when data is also being fit [37, 11]. However, this is a different meaning from how we use the terminology: without the existing data, they are simply the PINNs objectives, rather than a tool used to reduce overfitting, which is the sense with which we employ the term "regularizer" here. Gradient regularization has been applied to PINNs in a style similar to double backpropagation, though only on the domain objective [37]. Sobolev training was considered in Ref. [37]. While gradient regularized PINNs are simple and improve performance, they kept the theoretical analysis of why as an open question [37]. Sobolev PINNs are powerful, but require separate structure on a PDE by PDE basis [37]. Neither were applied to PINNs failure modes.
3 Experimental Setup
This section describes the problem formulation of PINNs, explaining how the network is trained to predict the solution to the PDE. Our network architecture itself follows closely to the original PINNs formulation, and uses 4 hidden layers, each with 128 neurons and \( \tanh\) activation functions. The network is initialized using the Glorot normal initialization [37], and training was done using the Limited Memory BFGS (L-BFGS) algorithm [37], which is also standard [11, 12]. Matching current failure modes literature, we use a grid collocation method over the domain, and linearly spaced points on the boundaries and initial conditions [20]. Experiments using FP32 were run on a Nvidia Titan RTX card, while experiments with FP64 precision were run on a Nvidia A100 GPU. In addition to the training collocation set, a much larger set of random collocation points over the domain is used as a test set, in addition to denser boundaries and initial conditions. (see Appendix A for details). We wrote our implementation using the JAX library [37].
3.1 Physics Informed Neural Networks
PINNs aim to train a deep neural network \( u_\theta\) , parametrized by \( \theta\) , to approximate the solution to a PDE. The PDE is represented by a differential operator \( \mathcal{F}\) acting on the true solution \( u\) , which for the time dependent problems considered here gives
(1)
where \( x\) is the \( d\) -dimensional spatial variable representing the input over the domain \( \Omega\) , and \( t\) the 1d temporal variable. The network learns to approximate \( u(x,t)\) by \( u_{\theta}(x,t)\) by minimizing the loss \( \mathcal{L}(u_\theta)\) , which consists of a domain loss \( \mathcal{L}_{\mathcal{F}}(u_\theta)\) , set of boundary losses \( \{\mathcal{L}_{\mathcal{B}i}(u_\theta)\}\) , and set of initial condition losses \( \{\mathcal{L}_{\mathcal{I}}(u_\theta)\}\) :
(2)
where \( N_\mathcal{B}\) , \( N_\mathcal{I}\) , are the number of boundary and initial conditions, which depend on the PDE of interest, and
(3)
are the domain, boundary, and initial losses, with \( n_d\) , \( n_b\) , and \( n_i\) referring to the number of collocation points sampled at domain, boundaries, and initial conditions. Where multiple boundary or initial conditions are required, we use the same number of collocation points for each of them. We follow the original PINNs formulation and set \( \lambda_\mathcal{F} = \lambda_\mathcal{B} = \lambda_\mathcal{I} = 1\) [11].
If a subset of data representing the real solution \( u\) is known, the above loss values are sometimes known as physics-based regularizers, as they enforce the PDE prior knowledge during optimization in addition to an extra loss term minimizing the error between the network and known solution points \( \| u(x_i,t_i) - u_\theta(x_i,t_i)\|^2\) [11, 37]. This is true while known data is present. However, as is common in modern PINNs, we work without prior knowledge or data of the solution, and the loss terms serve to enforce the correct domain, boundary, and initial conditions that the solution requires in order to satisfy the PDE. When discussing regularization later on, we use it in the sense of improving the performance of the model, rather than referring to the role of the core loss terms themselves.
Throughout this paper we use a regularizer we call double PINN, our extension of double backpropagation [23] to the full set of PINN losses:
(4)
Double-backpropagation regularizes the gradient of the loss with respect to the network inputs, and double PINN extends this to each PINN loss. By flattening the residuals, double PINN penalizes the network for reducing the losses at the collocation points at the expense of the surrounding region. We define it here for reference; the motivation is developed in Sec. 4.1 and a comparison to alternatives in Sec. 5.
4 Failure Modes are Caused by Overfitting
There are two leading explanations for the root cause of the failure modes. The first, is that there is an insurmountable loss-barrier, which traps the optimizer in a local minimum with high loss [12]. Recently a second proposal emerged, suggesting that instead of a local minimum the issue is unsuccessful optimization due to limited precision inhibiting navigation of the loss landscape during optimization [20].
In this section, we challenge both interpretations. By examining the residual as well as the loss during training, we show that the network has successfully been optimized to a global minimum with low loss, not a high loss basin, but that this has occurred by overfitting to the collocation points used during training. We then show that an upgrade to FP64 does not rescue the optimization.
4.1 Optimization Difficulty or Failure Mode
We consider two classic failure mode equations: the convection equation, and the wave equation. Both are relatively simple equations, standard in failure mode literature, and have known closed form solutions.
4.1.1 Convection Equation (FP32 Failure)
First we start with the 1D convection equation with periodic boundary conditions, which acts to transport the initial condition \( u(x,t=0)\) , with advection speed \( \beta\) . The problem is notable, because it is first order, has a simple solution, and can be readily solved by PINNs for small \( \beta\) , while failing if \( \beta\) is made large. The exact problem formulation can be found in App. B.1. For the advection speed \( \beta\) , a common choice is \( \beta =50\) , which we use here.
To examine the loss basin hypothesis, we first consider a typical training run that results in a failure mode. Figure 1(a) shows the expected exact solution, while Fig. 1(b) shows the result of the PINNs training. This is the prototypical failure mode, with the PINN outputting zero over most of the domain. The loss-basin hypothesis asserts that the PINNs is stuck in a local minimum, and that a true global minimum exists elsewhere [12]. Yet looking at the loss in Fig. 1(d) shows smooth optimization to \( \sim 10^{-1}\) , indicating successful optimization to a near-global minimum. It is the test loss, rather than the training loss, that remains large. While we have not seen a direct comparison between a test and train loss in the PINNs literature before, we consider this divergence the hallmark of a PINNs failure mode, and it is something that was consistently observed across all of our tests.
Fig. 1(c) further shows that even though the residual is large in the boundary region, it has been minimized at the collocation points. Rather than a local minimum, this indicates overfitting, where the PINN has learned the dataset (the training collocation points) at the expense of generality (satisfying the solution everywhere). The prescription is then regularization.
The same system with added regularizer is shown in Fig. 1 (e)-(g), where we no longer see overfitting at the collocation points, and the test and train loss decrease in tandem.
4.1.2 Wave Equation (64 bit failure)
What about the same-basin hypothesis, where the failure mode is suspected to be caused by insufficient numerical precision [20]? We consider this by examining the wave equation at both FP32 and FP64. The wave equation describes the time evolution of waves, and is detailed in App. B.2
Reference [20] reported successful optimization on the wave equation using a \( 101\times101=10201\) grid over the domain. A large amount of points makes it difficult for the network to memorize the data. If our hypothesis of overfitting is correct, we should expect reduced performance, and failure, as the number of collocation points is reduced, even at FP64.
Figure 3 shows FP64 PINNs trained on the wave equation with a reduced \( 12\times 12=144\) points. While successful on the \( 10201\) grid, FP64 PINNs fail on the \( 144\) grid. However, regularizing the same training run allows the optimization to succeed (Fig. 3(d)), even on such a small grid. The PDE residual loss plotted against a time slice of the collocation points is shown in (c) and (e), with the grey lines indicating their position. The unregularized PINN has notable minima exactly at the collocation points. The regularized version shows no such correlation between collocation and residual minima.
4.2 Regularization Combined with Precision is Optimal
Section 4.1.2 showed that precision is not sufficient to avoid failure, and that problems with a low amount of data remain susceptible to failure modes caused by overfitting. Yet increasing the precision can still be helpful, especially for stiff equations.
Consider the Allen-Cahn equation (App. B.4), which is sometimes considered a failure mode equation. We push back on this characterization, and instead suggest that optimization dynamics cause large error without hitting a failure mode. Figure 4 shows numerical solution of the Allen-Cahn equation in (a), along with a PINNs training run at FP32 (b). Here, we use a denser \( 4096\) point grid, in order to be able to capture the sharp interfaces in the solution. The FP32 PINN does not optimize well. Yet examining train and test losses show they move in lockstep, completely different behaviour from the observed failure modes of the convection and wave equations. This suggests an optimization difficulty, rather than a classic failure mode. If this is true, an upgrade to FP64 should indeed rescue optimization, without the need for a regularizer. This is exactly the behaviour seen in (c), with the unregularized FP64 PINN successfully modelling the equation. However, a late-stage increase in the test loss suggests that after learning the solution the model starts to memorize the data as optimization continues. In this case, a regularizer should be beneficial, even in the absence of a failure mode. Figure 4(d) shows this, with the best performance obtained through a combination of FP64 PINN and double PINN.
4.3 Optimizing the Performance of Double PINN
While we have shown that switching to higher precision is neither necessary nor sufficient for eliminating failure modes, it remains true that using higher precision will improve the accuracy of the PINN. We test double PINN at FP64 on the four benchmark equations: Convection, reaction, wave, and Allen-Cahn. We use a total number of \( 400\) , \( 256\) , \( 144\) , and \( 4096\) grid domain points for each equation, respectively. The convection equation uses \( 200\) points for each boundary and initial condition, while the Reaction, Wave, and Allen-Cahn equations use \( 100\) . Furthermore, we compare against the top three performing PINNs in the literature: unregularized vanilla FP64 PINN [20], PINNsFormer [21], and PINNMamba [22], all of which use a denser \( 10201\) grid.
| Model | 3Loss | rMAE | rRMSE | 2\( N\) |
| Convection | ||||
| 1PINNsFormer [21] | 9.0e-4 \( \pm\) 1.0e-4 | 3.3e-2 \( \pm\) 6.8e-3 | 4.4e-2 \( \pm\) 7.3e-3 | 10403 |
| 1PINNMamba [22] | 1.0e-4 \( \pm\) 2.0e-5 | 1.8e-2 \( \pm\) 3.7e-3 | 2.0e-2 \( \pm\) 3.8e-3 | 10403 |
| PINN_FP64 [20] | 5.0e-6 \( \pm\) 1.0e-6 | 5.9e-3 \( \pm\) 1.3e-3 | 7.2e-3 \( \pm\) 1.7e-3 | 10403 |
| double PINN (Ours) | 4.0e-11 \( \pm\) 4.1e-11 | 4.9e-6 \( \pm\) 4.4e-6 | 5.4e-6 \( \pm\) 4.6e-6 | 800 |
| Reaction | ||||
| 1PINNsFormer [21] | 3.0e-6 \( \pm\) 1.0e-6 | 1.5e-2 \( \pm\) 1.3e-3 | 3.0e-2 \( \pm\) 2.7e-3 | 10403 |
| 1PINNMamba [22] | 1.0e-6 \( \pm\) 1.0e-6 | 9.2e-3 \( \pm\) 1.7e-3 | 2.1e-2 \( \pm\) 3.6e-3 | 10403 |
| PINN_FP64 [20] | 1.0e-5 \( \pm\) 5.0e-6 | 2.7e-2 \( \pm\) 6.3e-3 | 5.0e-2 \( \pm\) 1.1e-2 | 10403 |
| double PINN (Ours) | 3.5e-12 \( \pm\) 3.5e-13 | 1.1e-5 \( \pm\) 2.3e-6 | 2.9e-5 \( \pm\) 8.0e-6 | 456 |
| Wave | ||||
| 1PINNsFormer [21] | 2.3e-2 \( \pm\) 1.7e-3 | 3.5e-1 \( \pm\) 8.7e-2 | 3.6e-1 \( \pm\) 8.7e-2 | 10504 |
| 1PINNMamba [22] | 2.0e-4 \( \pm\) 3.0e-5 | 1.9e-2 \( \pm\) 3.3e-3 | 2.0e-2 \( \pm\) 3.3e-3 | 10504 |
| PINN_FP64 [20] | 4.2e-5 \( \pm\) 1.6e-5 | 8.0e-3 \( \pm\) 3.2e-3 | 8.1e-3 \( \pm\) 3.1e-3 | 10504 |
| double PINN (Ours) | 3.5e-7 \( \pm\) 8.9e-8 | 3.8e-4 \( \pm\) 1.0e-4 | 3.8e-4 \( \pm\) 1.0e-4 | 444 |
| Allen–Cahn | ||||
| 1PINNsFormer [21] | 4.6e-1 \( \pm\) 2.9e-1 | 9.9e-1 \( \pm\) 4.0e-2 | 9.9e-1 \( \pm\) 4.2e-2 | 10504 |
| 1PINNMamba [22] | 2.7e-3 \( \pm\) 2.0e-4 | 1.4e-1 \( \pm\) 1.2e-2 | 2.7e-1 \( \pm\) 2.0e-2 | 10504 |
| PINN_FP64 [20] | 1.3e-5 \( \pm\) 4.0e-6 | 1.6e-2 \( \pm\) 3.6e-3 | 5.5e-2 \( \pm\) 1.1e-2 | 10504 |
| double PINN (Ours) | 4.0e-5 \( \pm\) 6.1e-10 | 9.6e-5 \( \pm\) 2.0e-5 | 3.6e-4 \( \pm\) 9.8e-5 | 4396 |
| 1 PINNsFormer and PINNMamba data, and high-resolution numerical solution for Allen-Cahn taken from Ref. [20] | ||||
| 2 \( N\) is the total number of collocation points (domain + initial + boundary). | ||||
| 3 Loss is not directly comparable across the methods due to the regularization penalty, but is included here for reference. | ||||
Figure 5 shows the result of this training for \( 10^6\) L-BFGS iterations on the four benchmark equations, all of which achieve excellent error. The final loss values, relative Mean Absolute Error (rMAE), and relative Root Mean Square Error (rRMSE) are shown in Table 1, along with the total number of collocation points \( N\) used in each. Double PINN at FP64 performs better than the compared methods in both error metrics by at least an order of magnitude, achieving a new state-of-the-art. Note the Allen-Cahn requires significantly more points, as it developed sharp interfaces that are difficult to resolve.
Finally, a note should be given on the computational cost. Double backpropagation requires taking the gradients of the loss terms with respect to the inputs, which theoretically doubles the computational cost. This is significant. Yet, as shown in Figure 6, double PINN is able to converge much faster as well. This, combined with the computational savings from using fewer collocation points, means that training of the regularized PINN becomes much cheaper in practice. This is especially true compared to architectures that use transformers or state space models, such as PINNsFormer and PINNMamba [21, 22].
5 Regularization of PINNs
Analysis in Sec. 4 points to overfitting being the underlying cause of PINNs failure modes, and that regularization with double PINN (Eq. (4)) can avoid such collapse while allowing the use of a vanilla PINN architecture with far fewer collocation points than normal. There are, however, several common choices of regularizers in the literature that might also be expected to work. This section discusses their performance when applied to PINNs in comparison to double PINN.
5.1 Regularization Strategies
We consider three regularization methods here. One common category are regularizers that shrink the size of the weights. Two of the most common of this type are \( L^2\) and \( L^1\) regularization, which penalize the norm of the combined sum of the weights. Specifically, if the weight parameters are denoted by \( \mathbf{w}\) , and the unregularized total loss as \( \mathcal{L}_0\) , then \( L^2\) regularizer adds the \( L^2\) norm of the weights as an auxiliary term, giving \( \mathcal{L}_\theta = \mathcal{L}_0 + (\lambda_r/2)\| \mathbf{w} \|^2_2\) where \( \lambda_r\) controls the regularization strength [33]. The \( L^1\) regularizer has the same form, swapping out the \( L^2\) norm for the \( L^1\) norm: \( \mathcal{L}_\theta = \mathcal{L}_0 + \lambda_r\| \mathbf{w} \|_1\) . Using an \( L^1\) norm promotes sparsity in the network. Our final and preferred method is double PINN, defined in Eq. (4).
The performance of each of these regularizers is shown on the convection equation with \( \beta=50\) in Figure 6. The unregularized case shows failure caused by overfitting. All three of the regularization methods remove the overfitting, which can be seen in both the test loss tracking the train loss, and in the residual plots. double PINN performs the best, achieving the lowest loss and fastest convergence, while the \( L^1\) regularization performs the worst. We suspect that the sparsity enforced by \( L^1\) reduces the expressiveness of the network too much. The \( L^2\) regularizer enforces no such sparsity, and overcomes the failure mode. We attribute the better performance of double PINN to the fact that it does not limit the network weights directly, but rather acts only on the residual and so does not limit the expressivity of the network.
6 Conclusions
Failure modes in PINNs have received considerable attention over the past several years, with research attempting to determine both their underlying causes and effective solutions. Here we have reanalyzed the problem, and shown that the underlying cause is overfitting. This changes the current understanding of failure modes, and suggests regularization as the straightforward prescription. We show that without it, even systems that aim to solve the failure modes (such as high precision PINNs), implicitly rely on overdetermining the problem, adding significant computational cost in the form of a large number of collocation points. By extending double-backpropagation to PINNs, we show that we can train efficiently, and with top tier performance, using a fraction of the points. We expect the advances here will have broad benefits in PINNs, particularly on complex and high dimensional PDEs.
Acknowledgments and Disclosure of Funding This work was partially supported by JST CREST (JPMJCR24Q5), ASPIRE (JPMJAP2329), and PRESTO (JPMJPR24TB), and JSPS KAKENHI (24K15105), Japan.
A Experimental Setup
Here the experimental setup is discussed, along with the training procedure. For the most part, we follow the standard vanilla description of PINNs [11].
A.1 Network Setup
Our implementation is written in JAX, using the Flax NNX library for the base network, the Optax library for the L-BFGS implementation, and the Orbax library for saving the model [37, 37, 37]. The network itself contains 4 hidden layers, each with 128 neurons. We use a hyperbolic tangent activation, and initialize the weights using the Glorot normal initializer [37]. Inputs to the network are normalized from \( [-1,1]\) , matching the original vanilla PINNs implementation.
The defaults for Optax’s L-BFGS optimizer vary slightly from those of the popular Pytorch version, with PyTorch’s optimizer containing a notably higher memory setting of 100, instead of Optax’s 10 [37, 37]. We increase the memory to 100 to match the Pytorch implementation, but leave the rest of Optax’s defaults alone.
The L-BFGS optimizer also contains no inherent stopping criteria. In our training algorithm, we therefore take the following approach: training proceeds with the relative difference in the current train loss noted every 100 iterations. We then compare it against the square root of the machine epsilon for whatever precision we were using on the run (either 32 bit floating point or 64 bit floating point). Taking the root of machine epsilon is standard in numerical optimization, and we find our method gives a good indication to when the optimizer is no longer making progress [37]. We use the difference over 100 runs to avoid pre-mature exiting while the optimizer is traversing the loss landscape. Together, this means we optimize as long as there is a relative difference in the loss of \( 3.453 \times 10^{-4}\) at FP32, and \( 1.490 \times 10^{-8}\) at FP64 over 100 iterations. Finally, we also specify a maximum number of epochs. This is relevant mostly for floating point 64, as the increased precision allows the optimizer to continuously make improvements for many iterations without stalling.
We use two sets of collocation points during training: a training and a test set. The training set consists of a small grid of collocation points. For example, Fig. 1 uses \( 200\) points for each of the boundaries, and \( 400\) points as a grid over the domain. The test set consists of a much larger set of random points, meant to sample the solution over a larger area to indicate its broader performance without relying on knowing the analytical solution. This test set uses \( 10000\) points randomly selected over the domain, and \( 10000\) linearly spaced points for each of the boundary and initial conditions. We expect the train and test set to move together for a well-trained model, while a lower train loss compared to the test loss indicates overfitting. A divergence between the two indicates a catastrophic failure, where the loss is minimized only at the collocation points, at the expense of the solution elsewhere. We consider this divergence the hallmark of a PINNs failure mode.
B Partial Differential Equations
The following section contains the problem formulations for each of the partial differential equations described in the main work.
B.1 1D Convection or Linear Advection
The 1D convection equation (sometimes called linear advection) is
(5)
with initial and boundary conditions
(6)
The analytic, closed form solution to the PDE is
(7)
For the advection speed \( \beta\) , a common choice in the failure mode literature is \( \beta =50\) , which we use here.
B.2 Wave Equation
The 1D wave equation with Dirichlet boundary conditions is:
(8)
with initial and boundary conditions
(9)
where \( \beta\) controls the frequency of the harmonic in the initial condition. We use \( \beta=3\) here as is standard. The PDE has a closed form analytic solution
(10)
B.3 Reaction Equation
1D reaction equation with periodic boundary conditions:
(11)
with initial and boundary conditions
(12)
The closed-form analytic solution is
(13)
a common choice for the \( \rho\) , the growth-rate coefficient, is \( \rho=5\) , which we adopt here as well. As with the case of linear advection, the problem becomes more difficult for PINNs the larger the value of \( \rho\) .
B.4 Allen-Cahn Equation
1D Allen-Cahn equation with periodic boundary conditions:
(14)
with initial and boundary conditions
(15)
The Allen-Cahn equation above is the only PDE we consider that does not have an exact solution. Instead, for comparison, we use a high-resolution solution from a spectral solver, obtained from the source code of Ref. [20].
C Limitations
Our method directly targets the mechanism of failure modes in PINNs. The \( L^2\) and \( L^1\) regularizers, while simple and cheap to implement, are not the strongest in terms of absolute performance and convergence speed. Our double backpropagation method, double PINN, drastically reduces the number of iterations needed for convergence, however requires a theoretical doubling of the computational cost due to the need to compute the extra gradients. This is made up for by the ability to use fewer points, and the more efficient training, but is still a noticeable limitation. While the method here removes failure modes, even when using single precision, it does not target the optimization procedure itself, meaning no benefit is gained on stiff equations that do not fail (though our method still improves late stage accuracy on these equations if they can be optimized effectively). Stiff equation and their complicated loss landscapes are different from failure modes, and have been discussed elsewhere in the literature.
D Broader Impacts
Our method improves the training of PINNs by enabling failure-free training with few collocation points, and without needing to go beyond single precision neural networks. We expect this to be important for high-dimensional equations, which cannot be tackled with solvers that require dense collocation (the curse of dimensionality). Many modern graphics processing units are unoptimized for double-precision throughput, so being able to use single-precision is important for utilizing current hardware to its fullest. As PDE solving is important in many disciplines, e.g. engineering, the methods introduced here could have indirect broader impacts if used in those contexts. However, we do not feel the need to highlight anything specific beyond its direct application to PDE solving.
References
[1] A hybrid neural network-first principles approach to process modeling AIChE J. 1992 38 10 1499–1511
[2] Solving high-dimensional partial differential equations using deep learning Proc. Natl. Acad. Sci. U.S.A. 2018 115 34 8505–8510
[3] Artificial neural networks for solving ordinary and partial differential equations IEEE Trans. Neural Netw. Learn. Syst. 1998 9 5 987–1000
[4] The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems Commun. Math. Stat. 2018 6 1 1–12
[5] Neural operator: Graph kernel network for partial differential equations arXiv preprint arXiv:2003.03485 2020
[6] Fourier neural operator for parametric partial differential equations International Conference on Learning Representations 2020
[7] Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems IEEE Trans. Neural Netw. Learn. Syst. 1995 6 4 911–917
[8] Physics-informed neural operator for learning partial differential equations ACM/IMS J. Data Sci. 2024 1 3 1–27
[9] Functional multi-layer perceptron: a non-linear tool for functional data analysis Neural Netw. 2005 18 1 45–60
[10] Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators Nat. Mach. Intell. 2021 3 3 218–229
[11] Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations J. Comput. Phys. 2019 378 686–707
[12] Characterizing possible failure modes in physics-informed neural networks Advances in Neural Information Processing Systems 2021 34 26548–26560
[13] When and why PINNs fail to train: A neural tangent kernel perspective J. Comput. Phys. 2022 449 110768
[14] A unified hard-constraint framework for solving geometrically complex pdes Advances in Neural Information Processing Systems 2022 35 20287–20299
[15] ProPINN: Demystifying propagation failures in physics-informed neural networks arXiv preprint arXiv:2502.00803 2025
[16] Mitigating Propagation Failures in Physics-informed Neural Networks using Retain-Resample-Release (R3) Sampling International Conference on Machine Learning 2023 7264–7302 PMLR
[17] Respecting causality for training physics-informed neural networks Comput. Methods Appl. Mech. Eng. 2024 421 116813
[18] Challenges in Training PINNs: A Loss Landscape Perspective International Conference on Machine Learning 2024 42159–42191 PMLR
[19] Preconditioning for physics-informed neural networks arXiv preprint arXiv:2402.00531 2024
[20] Fp64 is all you need: Rethinking failure modes in physics-informed neural networks Advances in Neural Information Processing Systems 2025 38 142949–142970
[21] PINN The Twelfth International Conference on Learning Representations 2024
[22] Sub-Sequential Physics-Informed Learning with State Space Model International Conference on Machine Learning 2025 69507–69525 PMLR
[23] Improving generalization performance using double backpropagation IEEE Trans. Neural Netw. 1992 3 6 991–997
[24] DGM: A deep learning algorithm for solving partial differential equations J. Comput. Phys. 2018 375 1339–1364
[25] Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations Commun. Math. Stat. 2017 5 4 349–380
[26] Neural-network-based approximations for solving partial differential equations Commun. Numer. Methods Eng. 1994 10 3 195–201
[27] On the eigenvector bias of Fourier feature networks: From regression to solving multi-scale PDEs with physics-informed neural networks Comput. Methods Appl. Mech. Eng. 2021 384 113938
[28] Self-adaptive physics-informed neural networks J. Comput. Phys. 2023 474 111722
[29] PINNACLE The Twelfth International Conference on Learning Representations 2024
[30] Solving differential equations with constrained learning International Conference on Learning Representations 2025 2025 44906–44948
[31] BWLer: Barycentric Weight Layer Elucidates a Precision-Conditioning Tradeoff for PINNs arXiv preprint arXiv:2506.23024 2025
[32] Dual cone gradient descent for training physics-informed neural networks Advances in Neural Information Processing Systems 2024 37 98563–98595
[33] Deep learning MIT press Cambridge 2016 1 2
[34] Probabilistic machine learning: an introduction MIT press 2022
[35] A simple weight decay can improve generalization Advances in Neural Information Processing Systems 1991 4
[36] Regression shrinkage and selection via the lasso J. R. Stat. Soc. Ser. B Stat. Methodol. 1996 58 1 267–288
[37]