{'fraction_non_alphanumeric': 0.031340544717017575, 'fraction_numerical': 0.007919644581804134, 'mean_word_length': 4.779793469159922, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Transport Based Refinement of Physics-Informed Neural Networks

In this paper, we propose a refinement strategy to the well-known Physics-Informed Neural Networks (PINNs) for solving partial differential equations (PDEs) based on the concept of Optimal Transport (OT). Conventional black-box PINNs solvers have been found to suffer from a host of issues: spectral bias in fully-connected architectures, unstable gradient pathologies, as well as difficulties with convergence and accuracy. Current network training strategies are agnostic to dimension sizes and rely on the availability of powerful computing resources to optimize through a large number of collocation points. This is particularly challenging when studying stochastic dynamical systems with the Fokker-Planck-Kolmogorov Equation (FPKE), a second-order PDE which is typically solved in high-dimensional state space. While we focus exclusively on the stationary form of the FPKE, positivity and normalization constraints on its solution make it all the more unfavorable to solve directly using standard PINNs approaches. To mitigate the above challenges, we present a novel training strategy for solving the FPKE using OT-based sampling to supplement the existing PINNs framework. It is an iterative approach that induces a network trained on a small dataset to add samples to its training dataset from regions where it nominally makes the most error. The new samples are found by solving a linear programming problem at every iteration. The paper is complemented by an experimental evaluation of the proposed method showing its applicability on a variety of stochastic systems with nonlinear dynamics.


Introduction
The problem of uncertainty propagation through nonlinear dynamical systems with stochastic excitations or uncertain parameters arises in a variety of fields, ranging from statistical physics to economy and finance [1,2]. The exact description of the time evolution of the probability density function (pdf) ρ(t, x) corresponding to the state x, can be analysed by studying the Fokker-Planck-Kolmogorov equation (FPKE) [3][4][5]. It is a partial differential equation (PDE) whose analytical solutions for different dynamical systems are difficult to obtain except under very special circumstances, for eg., systems with linear dynamics and a small class of nonlinear low-dimensional systems [6,7]. In this paper, we are concerned with solving the steady-state form of the FPKE for any nonlinear system.
There are several numerical techniques to solve the FPKE like the finite element method [8][9][10], finite difference method [11,12], path integral methods [13,14], that have emerged from the literature on solving PDEs. While these methods perform reliably for systems with low dimensions , they do not scale well and are known to become computationally intractable for high-dimensional systems. Refinement techniques [15,16] to combat the curse of dimensionality have been introduced in the past. The orthodox way to study uncertainty propagation is by using the Monte Carlo method [17,18]: to generate a large number of sample paths from a family of test points, propagate the nonlinear system forward exactly, and study the statistics at the time-steps of choice. While this certainly provides a reliable solution to the governing PDE, it too depends heavily on the amount of computing resources available to generate large amounts of data, plagued again by the curse of dimensionality. Alternative methods to analyse uncertainty propagation involve linearization of the nonlinear system in consideration [19], and the choice to ignore the higher-order terms can lead to significant errors. This is where the universal approximation power of neural networks is coming to the fore. With novel developments in architecture design and training strategies, they have been tremendously successful in the data-rich domains of computer vision, speech processing and the like, revolutionizing the way we build predictive models of complex systems today. In the function approximation world of physics-informed scientific computing, there has been remarkable progress in the way these black-box models are being harnessed to solve solve problems by endowing them with prior domain knowledge and appropriate inductive biases [20][21][22].
Of late, the dominant school of thought to train neural networks to solve differential equations is by introducing penalty terms into the loss functions that constrain the obtained solutions to satisfy some sets of governing conditions [23]. Specifically, physics-informed neural networks (PINNs) solve PDEs by enforcing a composition of constraints coming from the PDE and associated initial and boundary conditions. Unfortunately however, PINNs have been found to suffer from a host of issues: spectral bias from fully-connected architectures, instabilities during training, convergence and accuracy issues with "stiff" PDEs [24][25][26]. Numerous modifications to the baseline PINNs algorithm have been proposed and most of them revolve around adaptation mechanisms to strategically weight the different loss functions in the composite loss to account for differing gradients used during training.
Using the baseline PINNs for solving the FPKE is also not as straightforward. The FPKE is a formidable challenge to solve using a generic neural network because the resulting pdf has to satisfy two unique properties: positivity over the entire domain and a normalization condition.While positivity can be ensured through an appropriate choice of architecture, enforcing normalization through a penalty function [27] introduces another gradient pathology that can be hard to resolve.
To overcome these issues, we transform the stationary FPKE to an equivalent PDE that allows us to use the conventional PINNs framework. A state-independent diffusion vector allows us to generalize the solution to the true FPKE as the exponential of a potential function [4], i.e., ρ(x) = N 0 exp(−η(x)), where η(x) is the potential and N 0 is a normalizing coefficient. This transformation implicitly preserves positivity, and normalization can be performed using numerical integration over the entire domain after the PDE in η(x) has been solved.
Further, to overcome any instability issues that may arise during training, we propose a refinement strategy based on the theory of optimal transport (OT) [28] to actively choose collocation points that improve the training data regime over which the solution is learnt. The OT framework allows the synthesis of optimal mappings between functions on some measure spaces with respect to a variety of cost functions, depending on the nature of the problem. The theory links very well with stochastic dynamical systems, and has been used with great success in the past for solving Bayesian filtering problems [29].
Typically, collocation points are chosen through sampling from a pre-determined distribution or placed uniformly across the grid in consideration. Uniform grid distributions do not scale well with increasing dimensionality, nor do they account for the efficacy of the obtained solution over the entire domain. Alternatively, we seek to generate these points by examining the underlying distribution of the error regime from a nominally trained network. We solve a linear programming problem based on the theory of optimal transport to obtain these new points at the regions where the network fails the most, and iteratively train the network to a reasonable degree of accuracy. This can work in tandem with the adaptive weight tuning strategies proposed in the last few years, remarkably improving the utility of a small training set and reducing the computational overhead of training a highly complex network over each and every point in a very fine mesh.
The primary contributions of this paper are: • A novel iterative training scheme based on the theory of Optimal Transport to generate collocation points for PINNs in regions of high error, thereby improving the sample efficiency of a conventional PINNs solver • A PINNs framework to solve the stationary Fokker-Planck equation for any dynamical system through a convenient transformation, thus ensuring the requisite conditions of the obtained probability density function The paper is organized as follows: in section 2, we briefly describe the original PINNs formulation and the Fokker-Planck-Kolmogorov equation for any stochastic dynamical system. Section 3.1 details the reduced formulation for the FPKE that allows the conventional PINNs framework to yield satisfactory results. In section 3.2, we elaborate upon the proposed sampling strategy using Optimal Transport. This is succeeded by numerical simulations on a few systems, demonstrating the prowess of the formulations proposed. A helpful algorithm has been provided for reference. The paper ends with a summary of the findings and some questions posed for further research.

Physics Informed Neural Networks
Consider any nonlinear PDE of the form: where N t,x is a spatio-temporal differential operator. x ∈ Ω is a n-dimensional spatial vectorvariable operating in a domain Ω ∈ R n and t is time. ∂Ω refers to the boundary of the domain Ω. Following the work of [21], the objective is to approximate u(t, x) by a fully-connected neural network u θ (t, x) with inputs t and x, parameterized by θ. Define an additional residual network r θ (t, x) with the same weights θ such that: The optimal set of weights is obtained by training the neural network on a composite loss function of the form: where Here, L r (θ) corresponds to the loss function on the residual from the original PDE, and L b (θ) and L i (θ) relate to losses on the boundary and initial conditions respectively.

Fokker-Planck Equation
The probability density function of stochastic differential equations is governed by the Fokker-Planck equation. The stationary response of nonlinear dynamical systems can be obtained by solving the steady-state Fokker-Planck equation. It is a linear partial differential equation (PDE) whose analytical solutions have been derived for a restricted class of systems.
Assuming an n-dimensional stochastic differential equation (SDE) of the forṁ with x ∈ R n , Λ is an n × m matrix and Γ is a vector of uncorrelated standard Gaussian white-noise processes. We are interested in the evolution of the state PDF (probability density function) ρ(t, x) given by the Fokker-Planck equation, where D = 1 2 ΛΛ T . The probability density function ρ(t, x) gives the probability of being in a differential element (x, x + dx) of the phase plane at time t. It is subject to the normalization condition: The boundary conditions are given by ρ → 0 as |x i | → ∞ ∀i = 1, 2, . . . , n. This is inferred from the normalization condition, but alternative conditions like zero flux at infinity can be imposed as well [30].
Further, assuming a constant Λ reduces (4) to: In the steady state, the equation is: We seek to find the solution ρ θ (x) to (7) using PINNs. Alternatively, (7) can be expressed as It is a linear second-order parabolic PDE for which, a typical PINNs solver will return a trivial zero solution for ρ(x) during training. We circumvent this problem by transforming the PDE to an equivalent PDE after a change of variables. The resulting PDE can be solved in the conventional PINNs framework. Since we are solving the steady-state FPKE, the composite loss function is simply, 3 Proposed Formulation

Steady-state FPKE in η
From [4], the stationary solution to the standard FPKE can be expressed as: where η(x) is the potential function associated with the stationary solution to the Fokker-Planck equation. N 0 is a normalizing constant that ensures: lim t→∞ x ρ(x)dx = 1.
This implies, Consequently, (7) transforms to The solution to this PDE can now be obtained by using conventional PINNs methods.

Obtaining Collocation Points using Optimal Transport
The baseline PINNs solver has several modes of failure: instability during training [25,26], stiffness in gradient flow dynamics arising from imbalance in loss gradient magnitudes [24]. Most recent work in the realm of PINNs has been to mitigate issues of this sort by introducing modifiable weights to these loss functions, thus altering the loss function during training to arrive at a reasonably good approximation.
Typically, collocation points used to train the network are generated either on a uniform grid or sampled from a chosen distribution. The distributions are pre-determined and do not account for the efficacy of the obtained solution across the entire domain. We propose an iterative formulation in which we use the theory of Optimal Transport (OT) to add points to the training dataset in regions where the nominal PINNs solution fails the most.
All OT problems involve two probability measures µ and ν defined on some space Ω and a nonnegative cost function c(.) on Ω × Ω. The choices of Ω and c(.) subjectively depend on the problem being considered. Let D µ denote the support of the measure µ. The objective is to find a map φ : D µ → D ν such that the total cost transporting a single element x to φ(x) is minimized. This cost is given by: inf φ(.) c(x, φ(x))µ(dx). In the context of our proposed formulation, µ is a uniform distribution over the domain, whereas ν denotes the equation error distribution computed from the nominal solution at each step.
This transport map can be approximated by a linear map Φ, represented as a matrix Φ := [φ ij ] such that: where φ ij = p(I ij ), i.e., the probability that x − i,k is mapped to x + j,k and where the indicator I ij = 1 if x − i,k is mapped to x + j,k and 0 otherwise. N denotes the number of samples being transformed. x − i,k indicates the i th sample at iteration k and x + i,k denotes the transformed sample. Compactly, is the equally weighted ensemble from uniform distribution µ at iteration k. Similarly, the equally weighted transformed ensemble from the error distribution at iteration k is: To ensure that the matrix T := [t ij ] is measure-preserving, we need to enforce the following constraints: where w i is the equation error squared at each point x − i,k . The optimization problem can be rewritten as a linear programming problem of size N 2 , where D(x − i,k ,x + j,k ) is the Euclidean distance between x − i,k and x + j,k . The transformed samples x + j,k have weights equal to w i . The step, X + k = X − k Φ = X − k N T can be interpreted as a resampling.
Our algorithm detailed in 1 in section 4.2 makes use of this optimization to transform a uniformly weighted distribution to one that corresponds to the distribution of the square of equation error. We then add a chosen number of samples from the transformed distribution to the training dataset and re-train our PINNs network to find the optimal solution.

Numerical Simulations
We first demonstrate the PINNs formulation to solve the steady-state FPKE with the PDE in η(x) as the solution variable. In the succeeding subsection, we apply the OT-based sampling strategy to iteratively train our neural network with far fewer points than the chosen baseline and arrive at a similar degree of accuracy.
We use the same architecture for all the examples that follow: [d, 48, 1] where d is the size of the state-space of the system. A hyperbolic tangent activation function (tanh) is used between the input layer and the hidden layer, where as commonly observed in regression problems, a linear layer is used between the hidden layer and the final output layer. We used Julia [31][32][33][34][35] to run our simulations on the High Performance Research Computing facilities at Texas A&M University. To obtain the normalization constant N 0 , we rely on numerical integration schemes. We used the BFGS optimizer [36,37] with maximum number of training iterations set to 10000 for all cases.
There are two ways to measure the accuracy of the returned solution. For examples in which the analytical solution is available, the mean squared solution error evaluated over a finely discretized evaluation grid may be a fair criterion for judging correctness. We shall refer to this as ρ . While this remains suitable for most purposes, the true test of accuracy comes from measuring the residual equation error pde over the interior of the testing grid. Assuming N U is the number of points in the interior of the testing grid U, For systems without analytical solutions, we use pde to gauge the accuracy of the solution. As will be seen in the results, the stationary FPKE boundary conditions are easily satisfied.

Steady-state FPKE using Conventional PINNs
Consider a 2-dimensional Van der Pol -Rayleigh oscillator with dynamics given by: This can be rewritten as a system of differential equations in x = [x 1 , x 2 ] as: where Γ is standard Gaussian white noise. The FPKE in η(x) can be obtained from equation 11. The analytical solution to the stationary FPKE is available [6]: We consider the domain: 2]. We discretize the region with dx = [0.05, 0.05]. This corresponds to a training dataset of 6241 points inside the boundary of the training dataset. After running the neural network training to get the optimal sets of weights and biases minimising the composite loss function, we obtain a solution (figure 1) with pde = 1.92e −4 and ρ = 3.5e −5 .

OT-Based Sampling for Refinement
In this section of the paper, we demonstrate the effectiveness of the proposed OT-based refinement strategy on two examples. The first is the two-dimensional Van der Pol -Rayleigh oscillator, described above in section 4.1 whose solution η(x) is a polynomial in state space. The second is the classic Van der Pol oscillator, for which no analytical solution exists.
In the following algorithm, nOT refers to the number of OT iterations we wish to carry out with M points discovered in each iteration. f (x) is the dynamics of the nonlinear system in consideration without the diffusion term. S and U denote the interiors of the training and testing set respectively.  Transform U err to S OT using equation (14) ; S ← S ∪ S OT ; i = i + 1 ; θ i ← Train NN with S and θ i−1 as the initial guess ; end

Van der Pol -Rayleigh oscillator
We shall assume the solution obtained above using a uniform grid of dx = [0.05, 0.05] as the baseline. Our objective is to train a neural network with fewer points and obtain similar levels of accuracy without compromising on computational time.
For our nominal solution, we set up our network with the aforementioned architecture with dx = 0.25, which corresponds to 225 training points, a measly 3.6% of the baseline training dataset.
As expected, the solution obtained after training performs poorly when its accuracy is measured w.r.t the equation error on the finely discretized testing grid.
We then solve the Optimal Transport linear programming problem to obtain 200 new data points from regions in which the solution performs the worst. We add these points to our training dataset and subsequently train the neural network again.

Van der Pol Oscillator
Consider a 2-dimensional Van der Pol oscillator with dynamics given by: This can be rewritten as a system of differential equations in x = [x 1 , x 2 ] as: We shall solve for ρ(x) with σ 2 = 0.1. This is a more challenging example than the system described in   Table 1 summarizes the key benefits of the proposed method in terms of sample efficiency. Clearly, the OT-PINNs algorithm is significantly more data-efficient than the baseline PINNs framework and produces comparable results.

Conclusion
In this paper, we solved the stationary Fokker-Planck-Kolmogorov equation using Physics-Informed Neural Networks. We proposed and demonstrated an iterative strategy based on the concept of Optimal Transport to generate training datasets. This dramatically reduced the size of dataset required   to produce a solution with high accuracy. This will pay off when we seek to study uncertainty propagation for high-dimensional systems. Moreover, we bypassed the typical difficulties of solving the FPKE by solving an equivalent PDE in terms of the associated potential function of the probability density current for any given system.
However, there are still a few questions that remain unanswered in the proposed framework. Is there a requisite degree of accuracy that the nominal network must achieve for the OT-based strategy to succeed? How many points should the user choose to transform and how many OT iterations would it take to return a meaningfully accurate solution? Would increasing the number of hidden layers in the nominal network worsen the refinement strategy? We shall investigate these questions in our future works.