Monte Carlo Simulation of Stochastic Differential Equation to Study Information Geometry

Abstract

Information Geometry is a useful tool to study and compare the solutions of a Stochastic Differential Equations (SDEs) for non-equilibrium systems. As an alternative method to solving the Fokker–Planck equation, we propose a new method to calculate time-dependent probability density functions (PDFs) and to study Information Geometry using Monte Carlo (MC) simulation of SDEs. Specifically, we develop a new MC SDE method to overcome the challenges in calculating a time-dependent PDF and information geometric diagnostics and to speed up simulations by utilizing GPU computing. Using MC SDE simulations, we reproduce Information Geometric scaling relations found from the Fokker–Planck method for the case of a stochastic process with linear and cubic damping terms. We showcase the advantage of MC SDE simulation over FPE solvers by calculating unequal time joint PDFs. For the linear process with a linear damping force, joint PDF is found to be a Gaussian. In contrast, for the cubic process with a cubic damping force, joint PDF exhibits a bimodal structure, even in a stationary state. This suggests a finite memory time induced by a nonlinear force. Furthermore, several power-law scalings in the characteristics of bimodal PDFs are identified and investigated.

1. Introduction

Stochastic Differential Equations (SDEs) (Equation (5)) are used to model various phenomena in nature, including Brownian motion, asset pricing, population dynamics, COVID-19 spread and interaction [1,2,3,4,5,6], and various other non-equilibrium processes. Due to their stochasticity, SDEs do not have an unique solution, but a distribution of solutions. A Fokker–Planck Equation (FPE) [7] is a Partial Differential Equation (PDE) that describes how the probability density of solutions of a SDE evolves with time.

Comparing solutions of different SDEs can be achieved by looking at different statistics of the solutions like mean and variance. However, when we are interested in large fluctuations and extreme events in the solutions, simple statistics might not suffice. In such cases, quantifying and comparing the time evolution of probability density functions (PDFs) of solutions will provide us with more information [8]. The time evolution of PDFs can be studied and compared through the framework of information geometry [9], wherein PDFs are considered as points on a Riemannian manifold and their time evolution can be considered as a motion on this manifold. In general, in order to have a manifold structure on the probability space in information geometry, we need to define a metric. Several different metrics can be defined on a probability space [10,11,12,13].

Different metrics have different physical and mathematical significance. For example, the Wasserstein metric (also known as the Earth mover’s distance) naturally comes up in optimal transport problems [14]; the Ruppeiner metric is based on the geometry of equilibrium thermodynamics [13]. In this work, we use a metric related to the Fisher Information [15], known as the Fisher Information metric [16,17].

$$g_{jk}(\theta ):={\int }_X\frac{\partial logp(x;\{\theta \})}{\partial {\theta }_j}\frac{\partial logp(x;\{\theta \})}{\partial {\theta }_k}p(x;\{\theta \})dx$$

(1)

Here, $p(x;\{\theta \})$ denotes a continuous family of PDFs parametrized by parameters $\{\theta \}$. This metric was used to physically represent the number of statistically distinguishable states [11,18,19]. Note that two Gaussians with same standard deviation but different means are statistically indistinguishable, if the difference in their means are much smaller than their standard deviation.

If a time-dependent PDF $p(x,t)$ is considered as a continuous family of PDFs parameterized by a single parameter, time t, the scalar metric can be given by:

$$g(t)=\int d\mathbf{x}\frac{1}{p(\mathbf{x},t)}{\left[\frac{\partial p(\mathbf{x},t)}{\partial t}\right]}^2$$

(2)

However, time in classical mechanics is a passive quantity that cannot be changed by an external control. The infinitesimal distance $d\mathcal{L}$ on the manifold is then given by $d{\mathcal{L}}^2=g(t)d{t}^2$. Here, $\mathcal{L}$ is the Information Length defined by:

$$\mathcal{L}(t):={\int }_0^{t}d{t}_1\sqrt{\int d\mathbf{x}\frac{1}{p(\mathbf{x},t_1)}{\left[\frac{\partial p(\mathbf{x},t_1)}{\partial t_1}\right]}^2}$$

(3)

The Information Length $\mathcal{L}$ represents the dimensionless distance, which measures the total distance traveled on the manifold, or the total number of distinguishable states a system passes through during the course of its evolution. It has previously been used in optimization problems [20]. It was also used to study dynamical systems, thermodynamics, phase transitions, memory effects, and self-organization [21,22,23,24,25,26,27,28].

The gradient of $\mathcal{L}$, ${lim}_{dt\to 0}d\mathcal{L}/dt\equiv \mathsf{\Gamma }$, then represents a velocity on this manifold

$$\mathsf{\Gamma }(t):=\sqrt{\int d\mathbf{x}\frac{1}{p(\mathbf{x},t)}{\left[\frac{\partial p(\mathbf{x},t)}{\partial t}\right]}^2}$$

(4)

Note that we use the notation $\mathsf{\Gamma }$ instead of $\sqrt{g}$ to make it clear that it is a quantity defined for a time-dependent PDF. $\mathsf{\Gamma }$ represents the rate of change of statistically distinguishable states in a time-evolving PDF, and is sometimes referred to in the literature [29,30,31] as the Information Rate. Note that in information theory [32], the term “Information Rate” is used for the rate at which information is passed over the channel [33].

As for the physical significance of $\mathsf{\Gamma }$, in a non-equilibrium thermodynamic system, $\mathsf{\Gamma }$ is related to the entropy production [30]. $\mathsf{\Gamma }$ has also been used to study causality [29] and abrupt changes in the dynamics of a system [8]. ${\mathsf{\Gamma }}^2$ is equivalent to the (symmetric) KL divergence of infinitesimally close PDFs, as shown in Appendix E. It should also be noted that $\mathsf{\Gamma }(t)$ defined by Equation (4) has the dimensions of $t^{-1}$, and the time-integral of $\mathsf{\Gamma }(t)$ gives a dimensionless distance $\mathcal{L}$ in Equation (3).

Due to the lack of general mathematical techniques to solve SDEs or its associated FPE, analytical study of SDEs using Information Geometry ($\mathsf{\Gamma }$ and $\mathcal{L}$) has been limited to a few special cases [24,25,34,35]. To date, numerical studies have relied on solving the associated FPE [24,34], which has the advantage of generating smooth time-dependent PDFs and information diagnostics, but has the limitations outlined in

Table 1. Comparison between grid-based FPE solver and MC SDE simulation.

	Grid-Based FPE Solver	MC SDE Simulation
Accuracy	Depends on grid size. No well-defined prescription on choosing grid-size.	Depends on number of samples(n) [36]. Typically less accurate for practical sample sizes.
Boundary condition	Requires carefully chosen non-trivial boundary conditions. Cannot handle discontinuous initial conditions such as Dirac delta function.	Requires only an initial distribution as boundary condition.
Memory Usage & Runtime	Scales exponentially with dimension d. $\mathcal{O}\left(n_1{n}_2\dots n_d\right)$. Here, $n_1,n_2,\dots ,n_d$ are the number of grid points along each dimension.	Scales linearly with dimension d. $\mathcal{O}\left(nd\right)$. Here, n is the number of samples.
Correlation Study	Cannot study correlations and associated memory effects using FPE.	Can study correlations. See Section 4 for unequal time joint PDF estimates.

. To overcome the limitations of a FPE solver, in this work, we develop a Monte Carlo (MC) method to study time-dependent PDFs and the Information Geometry of SDEs.

The main aims of this paper are twofold. The first aim is to develop a new MC SDE simulation method and validate it by recovering the previous results obtained using the FPE method. The second is to calculate unequal time joint PDFs and investigate the effect of nonlinear forces on PDF form and various power-scaling relations. The remainder of the paper is organized as follows: Section 2.1 gives a brief introduction of the theory of MC SDE simulation. Section 2.2 develops the methods to measure Information Geometry from the simulation. Using this method, we compare a linear and a nonlinear SDE in Section 3. In Section 4, we showcase the measurement of joint PDFs of the same variable but at unequal times, which is not possible using FPE solvers. We then study a type of phase transition in the joint PDF of the nonlinear SDE and numerically verify its theoretically calculated scaling relations. Discussions are found in Section 5.

2. Methods

2.1. SDE Simulation

A typical SDE for the variable x in d dimensions has the following form:

$$d{\overrightarrow{\mathbf{x}}}_t=\overrightarrow{\mathbf{\mu }}({\overrightarrow{\mathbf{x}}}_t,t)dt+\mathbf{\sigma }({\overrightarrow{\mathbf{x}}}_t,t).d{\overrightarrow{\mathbf{W}}}_t$$

(5)

Here, $\overrightarrow{\mathbf{\mu }}$ is known as the drift vector (drift coefficient in 1D), $\mathbf{D}:=\mathbf{\sigma }.{\mathbf{\sigma }}^T/2$ the diffusion tensor (diffusion coefficient in 1D) and ${\overrightarrow{\mathbf{W}}}_t$ is the Weiner process [37]. $d{\overrightarrow{\mathbf{W}}}_t$ represents an infinitesimal random noise term, making the equation stochastic. Generalizations to SDEs can be achieved with more general noise terms and higher-order derivative terms, but are not pursued here. The associated FPE describes the time evolution of the PDF $p(\overrightarrow{\mathbf{x}},t)$ of solutions of the SDE.

$$\frac{\partial p(\overrightarrow{\mathbf{x}},t)}{\partial t}=-\sum _{i=1}^d\frac{\partial }{\partial x_i}\left[{\mathbf{\mu }}_i(\overrightarrow{\mathbf{x}},t)p(\overrightarrow{\mathbf{x}},t)\right]+\sum _{i=1}^d\sum _{j=1}^d\frac{{\partial }^2}{\partial x_i\partial x_j}\left[D_{ij}(\overrightarrow{\mathbf{x}},t)p(\overrightarrow{\mathbf{x}},t)\right]$$

(6)

Instead of numerically solving the FPE, in this work, we use a Monte Carlo (MC) method to estimate $p(\overrightarrow{\mathbf{x}},t)$ by simulating a large number of instances of an SDE. Explicit numerical solution of an SDE involves choosing an initial position ${\overrightarrow{\mathbf{x}}}_0$ and iteratively updating it to get the value at time t, ${\overrightarrow{\mathbf{x}}}_t$. For the MC simulation, we start with a set of initial positions $\{{\overrightarrow{\mathbf{x}}}_0\}$ sampled from the desired initial distribution and numerically solve each of them independently. We can then use the set of samples at time t, $\{{\overrightarrow{\mathbf{x}}}_t\}$, to compute the desired statistics. We can formally write this iteration step as follows:

$${\overrightarrow{\mathbf{x}}}_{t+\mathsf{\Delta }t}=\mathcal{M}\left({\overrightarrow{\mathbf{x}}}_t,\overrightarrow{\mathbf{\mu }}({\overrightarrow{\mathbf{x}}}_t,t),\mathbf{\sigma }({\overrightarrow{\mathbf{x}}}_t,t),\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_t,\mathsf{\Delta }t\right)$$

(7)

Here, $\mathcal{M}$ denotes an arbitrary update scheme. There are several methods [38] by which we can create this update. Throughout this work, we will be working with autonomous SDEs ($\overrightarrow{\mathit{\mu }}$ and $\mathit{\sigma }$ do not explicitly depend on time) with diagonal noise ($\mathit{\sigma }$ is a diagonal matrix). Therefore, we will use a simple update scheme called the Milstein scheme [39]:

$${\overrightarrow{\mathbf{x}}}_{t+\mathsf{\Delta }t}={\overrightarrow{\mathbf{x}}}_t+\overrightarrow{\mathbf{\mu }}({\overrightarrow{\mathbf{x}}}_t)\mathsf{\Delta }t+\mathbf{\sigma }({\overrightarrow{\mathbf{x}}}_t).\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_t+\frac{1}{2}\mathbf{\sigma }({\overrightarrow{\mathbf{x}}}_t).{\mathbf{\sigma }}^{\prime }({\overrightarrow{\mathbf{x}}}_t).\left({\left(\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_t\right)}^2-\mathsf{\Delta }t\right)+\mathcal{O}\left(\mathsf{\Delta }{t}^{1.5}\right)$$

(8)

Here, $\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_t$ is a random sample from ${\mathcal{N}}^d\left(0,\mathsf{\Delta }t\right)$. Note that the Milstein scheme is only accurate up to $\mathcal{O}(\mathsf{\Delta }t)$, and this error can accumulate through the course of the simulation. Therefore, to control the error in the numerical update scheme, we will adaptively choose the time step $\mathsf{\Delta }t$ for updates by setting a local error tolerance. Local error at time t is defined as the deviation between a single step update made with time step $\mathsf{\Delta }t$ and a two-step update made with the time step $\mathsf{\Delta }t/2$ each.

$$err({\overrightarrow{\mathbf{x}}}_t,\mathsf{\Delta }t):=\frac{∥\left(\mathcal{M}\left({\overrightarrow{\mathbf{x}}}_t,\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_1,\frac{\mathsf{\Delta }t}{2}\right)+\mathcal{M}\left({\overrightarrow{\mathbf{x}}}_{t+\frac{\mathsf{\Delta }t}{2}},\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_2,\frac{\mathsf{\Delta }t}{2}\right)\right)-\mathcal{M}\left({\overrightarrow{\mathbf{x}}}_t,\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_t,\mathsf{\Delta }t\right)∥}{\sqrt{d}}$$

(9)

Here, we omitted the dependence of $\mathcal{M}$ on $\overrightarrow{\mathit{\mu }}$ and $\mathit{\sigma }$ for brevity. $\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_1$ and $\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_2$ are chosen in such a way that $\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_1+\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_2=\mathsf{\Delta }{\overrightarrow{\mathbf{W}}}_t$. To satisfy the local error tolerance $Tol$, we need to choose $\mathsf{\Delta }t$, such that $err({\overrightarrow{\mathbf{x}}}_t,\mathsf{\Delta }t)<Tol$. The exact prescription on how this choice is made can be found here [40].

Computing $\mathsf{\Gamma }$ requires computing the derivative of $P({\overrightarrow{\mathbf{x}}}_t,t)$ of the numerical solutions. Since derivatives are sensitive to numerical noise, we need accurate estimates of the probability distribution. This is achieved by numerically solving a large number of SDEs (we use at least $2\times {10}^7$ samples in this work). This is an impractically large number for most computers. Numerically finding $2\times {10}^7$ solutions, with each requiring around 10,000 time steps (depending on the required accuracy) will require at least 1600 GB of memory, assuming 64-bit floating point values. Additionally, assuming a fast 2 ms per solution, the entire simulation will take around 11 h if carried out serially. To solve these problems, we use GPU computing. With GPU computing [41], we can perform updates on a set of values $\{{\overrightarrow{\mathbf{x}}}_t\}$, simultaneously using Equation (7) to get $\{{\overrightarrow{\mathbf{x}}}_{t+\mathsf{\Delta }t}\}$. $\{{\overrightarrow{\mathbf{x}}}_t\}$ can then be removed from memory after computing the desired statistics, making it memory efficient. As for the runtime, a GPU-based parallel implementation [42] in Python takes around 4 min to simulate $2\times {10}^7$ samples for 10,000 time steps on a consumer laptop equipped with Nvidia RTX 2080 GPU. See Appendix A for detailed scaling relations of simulation runtime.

2.2. Estimating $\mathsf{\Gamma }$

The form of $\mathsf{\Gamma }$ in Equation (4) makes it unsuitable for numerical computation due to the presence of $p(\mathbf{x},t)$ term in the denominator, which can become zero. We therefore rewrite the equation using the redefinition $q(\mathbf{x},t):=\sqrt{p(\mathbf{x},t)}$.

$${\mathsf{\Gamma }}^2(t)=4\int d\mathbf{x}{\left[\frac{\partial q(\mathbf{x},t)}{\partial t}\right]}^2$$

(10)

To calculate $\mathsf{\Gamma }$, we first estimate the PDF using histograms. (It would be more accurate to use kernel density estimators [43], but that is more computationally expensive). The derivative in Equation (10) is approximated as a finite difference and the integral as a Riemann sum. The specific methods and the error estimates are provided in Appendix B. After computing $\mathsf{\Gamma }$, information length can be computed from the numerical integration of Equation (3).

A major source of error in the estimation of $\mathsf{\Gamma }$ is due to the compact support and overlap of numerically estimated probability densities. (This will also be a source of error with numerical FPE solvers). During the simulation, time steps are chosen adaptively by setting a tolerance for the local numerical error. However, in some regimes, this results in the system evolving too fast, so that there is little or no overlap between the probability densities at adjacent time steps (Figure 1). Theoretically, these densities (Figure 1) have the support of the entire real line, and the only error will be caused by approximating the derivative and the integral. However, in practice, we are running the MC simulation with a finite sample size, and computers have finite numerical precision; as such, the estimated densities have compact support. Figure 2 shows how the amount of overlap between the densities affects the error in $\mathsf{\Gamma }$ estimate. There are two main sources of sub-optimal overlap: changes in mean and changes in standard deviation.

Consider $\mathsf{\Delta }t$ chosen adaptively to satisfy the local error tolerance. Now assume for some $\widetilde{\mathsf{\Delta }t}$, we get optimal overlap. If $\widetilde{\mathsf{\Delta }t}>\mathsf{\Delta }t$, we can wait for a few steps before estimating the $\mathsf{\Gamma }$. However, if $\widetilde{\mathsf{\Delta }t}<\mathsf{\Delta }t$, we can choose a $\widetilde{\mathsf{\Delta }t}$ as a temporary time step and generate a collection of points $\{{\overrightarrow{\mathbf{x}}}_{t+\widetilde{\mathsf{\Delta }t}}\}$ from $\{{\overrightarrow{\mathbf{x}}}_t\}$ using Equation (7), which results in optimal overlap between densities. The local error will still be smaller than the tolerance since $\mathsf{\Delta }t<\widetilde{\mathsf{\Delta }t}$. In order to derive the value of $\widetilde{\mathsf{\Delta }t}$, we use the Milstein update scheme and restrict ourselves to the one-dimensional case for simplicity. First, we look at the change in mean value.

$$x_{t+\mathsf{\Delta }t}=x_t+\mu (x_t)\mathsf{\Delta }t+\sigma (x_t)(x_t)\mathsf{\Delta }{W}_t+\frac{1}{2}\sigma (x_t){\sigma }^{\prime }(x_t)\left(\left(\mathsf{\Delta }{W}_t^2\right)-\mathsf{\Delta }t\right)+\mathcal{O}\left(\mathsf{\Delta }{t}^{1.5}\right)$$

(11)

Taking the expectation value of both sides, we get:

$$\mathbb{E}\left[x_{t+\mathsf{\Delta }t}\right]=\mathbb{E}\left[x_t\right]+\mathbb{E}\left[\mu (x_t)\right]\mathsf{\Delta }t+\mathcal{O}\left(\mathsf{\Delta }{t}^{1.5}\right)$$

(12)

Here, we use the fact that $\mathbb{E}\left[\mathsf{\Delta }{W}_t\right]=0$ and $\mathbb{E}\left[{\left(\mathsf{\Delta }{W}_t\right)}^2\right]=\mathsf{\Delta }t$. Note that $x_t$ and $\mathsf{\Delta }{W}_t$ are independent random variables. For X and Y, the independent random variables are $\mathbb{E}\left[XY\right]=\mathbb{E}\left[X\right]\mathbb{E}\left[Y\right]$.

In order to achieve optimal overlap (Figure 2 (left)), we need $\mathbb{E}\left[x_{t+\mathsf{\Delta }t}\right]-\mathbb{E}\left[x_t\right]=0.2Std\left[x_t\right]$.

$${\widetilde{\mathsf{\Delta }t}}_{mean}=\frac{0.2Std\left[x_t\right]}{\mathbb{E}\left[\mu (x_t)\right]}$$

(13)

Now to derive the effect on change in standard deviation on $\widetilde{\mathsf{\Delta }t}$, a similar calculation for $Var\left[x_{t+\mathsf{\Delta }t}\right]$ can be performed, which yields:

$$Var\left[x_{t+\mathsf{\Delta }t}\right]=Var\left[x_t\right]+\left(2Cov\left[\mu (x_t),x_t\right]+\mathbb{E}\left[{\sigma }^2(x_t)\right]\right)\mathsf{\Delta }t+\mathcal{O}\left(\mathsf{\Delta }{t}^{1.5}\right)$$

(14)

Now, in order to achieve optimal overlap, we need $Var\left[x_{t+\mathsf{\Delta }t}\right]/Var\left[x_t\right]=0.9^{\pm 2}$ (Figure 2 (right)). $0.9^{+2}$ when $Var\left[x_{t+\mathsf{\Delta }t}\right]<Var\left[x_t\right]$ and $0.9^{-2}$ when $Var\left[x_{t+\mathsf{\Delta }t}\right]>Var\left[x_t\right]$.

$${\widetilde{\mathsf{\Delta }t}}_{std}=Var\left[x_t\right]\left|\frac{0.9^2-1}{2Cov\left[\mu (x_t),x_t\right]+\mathbb{E}\left[{\sigma }^2(x_t)\right]}\right|$$

(15)

Note that we have only considered the first and second moments here. It is potentially possible to improve the accuracy of the $\mathsf{\Gamma }$ estimate by considering higher moments. However, this improvement will be marginal, since overlap between two distributions is most affected by its first and second moments.

After calculating both the time steps ${\widetilde{\mathsf{\Delta }t}}_{mean}$ and ${\widetilde{\mathsf{\Delta }t}}_{std}$, we take the minimum of the two to perform the update on $\{x_t\}$ to get $\{x_{t+\widetilde{\mathsf{\Delta }t}}\}$, where $\widetilde{\mathsf{\Delta }t}=min\left({\widetilde{\mathsf{\Delta }t}}_{mean},{\widetilde{\mathsf{\Delta }t}}_{std}\right)$. After the estimation, the $\mathsf{\Gamma }$$\{x_{t+\widetilde{\mathsf{\Delta }t}}\}$ is discarded and we continue the simulation with $\{x_{t+\mathsf{\Delta }t}\}$. This prevents any significant slowing in the simulation if $\widetilde{\mathsf{\Delta }t}\ll \mathsf{\Delta }t$. Note that calculating ${\widetilde{\mathsf{\Delta }t}}_{mean}$ and ${\widetilde{\mathsf{\Delta }t}}_{std}$ using the entire set of points $\{x_t\}$ will be computationally expensive. Only a small subset is used to perform this calculation.

After $\mathsf{\Gamma }$ is estimated, we can calculate the information length by approximating the integral as a Riemann sum.

Looking at Figure 3, it can be seen than the percentage error in $\mathsf{\Gamma }$ blows up towards the end of the simulation. This is when the system approaches a stationary state and the probability density stops evolving. The exact $\mathsf{\Gamma }$ reaches zero, whereas the numerical estimate will have a small nonzero error (Figure 4). Even though the percentage error in $\mathsf{\Gamma }$ becomes large, the absolute error remains small (Figure 2 (left)) and will have minimal contribution to the error in Information Length calculation. However, when the initial distribution is ’closer’ to the stationary distribution (small $x_0$ values in Figure 3), the absolute value of Information Length will be small, and the error in $\mathsf{\Gamma }$ estimate will have a more significant contribution to the Information Length calculation.

It is to be noted that measuring $\mathsf{\Gamma }$ of multi-dimensional problems is a significant computational challenge which requires further investigation. MC SDE simulation is better at handling such problems compared with FPE solvers due to its linear scaling with number of dimensions. However, this comes at the price of accuracy when estimating PDFs of higher dimensional problems. It is still possible to accurately study $\mathsf{\Gamma }$ from marginal PDFs of multi-dimensional problems using MC SDE simulation.

A python implementation of SDE Simulation, along with the $\mathsf{\Gamma }$ measurement used in this work, can be found here [42].

3. Linear vs. Cubic Statistics

In this section, we will use methods developed in the previous section and study nonlinear damping of the Information Geometry of a stochastic process. We will compare the Ornstein–Uhlenbeck process [44], a model for prototypical linear driven dissipative process, defined by linear SDE:

$$d{x}_t=-\theta x_{t}dt+\sigma d{W}_t$$

(16)

with a cubic SDE defined by the equation:

$$d{x}_t=-\theta x_t^{3}dt+\sigma d{W}_t$$

(17)

The cubic damping term has been previously used to model ${\mathrm{CO}}_2$ emissions [45], phase transitions [24], and self-organized shear flows [46].

Throughout this section, unless otherwise stated, we restrict ourselves to the values $\theta =1.0$ and $\sigma =0.1$. The initial distribution is always $\mathcal{N}\left(x_0,{10}^{-10}\right)$, and we run the SDE simulation from $t=0$ to $t=100$. The data in this section were generated using $2\times {10}^7$ samples, and PDF was approximated by a histogram with 703 bins, which was chosen using an empirical formula $2.59\sqrt[3]{n}$. The time steps were adaptively chosen by setting local error tolerance at $5\times {10}^{-4}$.

Two of the simplest statistics that can be measured are the mean and the standard deviation of the distribution (Figure 5 and Figure 6). The trends in mean can be readily seen by taking expectation value on both sides of Equations (16) and (17). For the linear SDE, we have $d\mathbb{E}\left[x_t\right]=-\theta \mathbb{E}\left[x_t\right]dt$, which can be solved to get $\mathbb{E}\left[x_t\right]=x_0{e}^{-t}$. For the Cubic SDE, we can follow similar steps:

$$\begin{array}{cc}\hfill d\mathbb{E}\left[x_t\right]& =-\theta \mathbb{E}\left[x_t^3\right]dt\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill d\mathbb{E}\left[x_t\right]& \approx -\theta \mathbb{E}{\left[x_t\right]}^{3}dt\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \Rightarrow \mathbb{E}\left[x_t\right]& \approx \frac{1}{\sqrt{\frac{1}{x_0^2}+2\theta t}}\hfill \end{array}$$

(20)

The approximate solution of the mean value of Cubic SDE is only valid when the standard deviation of the distribution is much smaller than the mean (Appendix C). Note that in Equation (20), when $1/x_0^2\ll \theta t$, $\mathbb{E}\left[x_t\right]\approx {(2\theta t)}^{-0.5}$, making the trajectory independent of $x_0$. This can be seen in Figure 5 (right), when all the lines merge into one. However, around $t=5$, the approximation fails, since the value of the standard deviation (Figure 6 (right)) and the mean (Figure 5 (right)) becomes comparable. Therefore, the curves fail to follow the same trajectory $\mathbb{E}\left[x_t\right]\approx {(2\theta t)}^{-0.5}$ for $t>5$.

The trends in standard deviation of the Linear SDE are readily explained by Equation (A14), which does not depend on initial mean position $x_0$, but only on the initial standard deviation, the drift coefficient, and the diffusion coefficient. As we write this paper, no exact analytic solution exists for the Cubic SDE. However, an approximate analytic treatment can be found here [47].

We define the asymptotic Information Length (${\mathcal{L}}_{\infty }$) to be the Information Length it took for the system to reach the stationary state of its PDF.

$${\mathcal{L}}_{\infty }:=\mathcal{L}(t\to \infty )$$

(21)

Analytically, a SDE reaches its stationary state as $t\to \infty$. However, numerically, we see that the probability density stops evolving after a finite time. We can see this from Figure 7, as $\mathcal{L}$ becomes a constant. $\mathcal{L}$ will still continue to increase slightly due to numerical error, but this contribution will be negligible, as shown in Figure 3. ${\mathcal{L}}_{\infty }$ measures how many statistically distinguishable states the system passes through to reach its stationary state. From Figure 6, it is evident that, compared to the linear SDE, the PDF of the cubic SDE undergoes a lot more change before reaching its stationary state, for larger values of $x_0$. For small values of $x_0$, the trend in standard deviation is similar between linear and cubic SDEs, since the initial evolution of the cubic process is Gaussian. This is confirmed in Figure 7, which shows that for large values of $x_0$, ${\mathcal{L}}_{\infty }$ is significantly larger for the cubic SDE for same values of $x_0$ and, for small $x_0$ values, the ${\mathcal{L}}_{\infty }$ values are comparable.

Information Length Scaling

In Figure 7, we have already seen that ${\mathcal{L}}_{\infty }$ depends on $x_0$ and has different behavior for linear and cubic SDEs. In [47], by numerically solving the FPE (also analytically for the linear SDE), it was shown that for large values of $x_0$, ${\mathcal{L}}_{\infty }$ shows different scaling behavior for linear and cubic SDEs. For the linear SDE, ${\mathcal{L}}_{\infty }\sim x_0$. For the cubic SDE, ${\mathcal{L}}_{\infty }\sim x_0^m$, where $1.5<m<1.9$. We reproduce this result in Figure 8. Note that, since ${\mathcal{L}}_{\infty }$ is a dimensionless quantity, it is not possible to derive these values theoretically using dimensional analysis. In the absence of general analytic tools to study this property, numerical methods are indispensable.

4. Unequal Time Joint PDF: Bimodality for Cubic Force

A clear advantage of MC SDE Simulation over solving FPE is the ability to compute the unequal time joint PDFs $P\left(x\left(t_1\right),x\left(t_2\right)\right)$. Unequal time joint PDFs have been previously used for causality and to establish causal relations. In this section, we showcase the ability of MC SDE simulation to estimate $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ and study its properties.

4.1. Unequal Time PDFs in the Stationary State

For the linear SDE $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ is always a Gaussian, but with a covariance matrix that depends on $t_1$ and $t_2$. However, for the cubic SDE, we see the emergence of a bimodal distribution (Figure 9 (right)) depending on $t_1$ and $t_2$ values. This behavior prevails even after the system has reached its stationary state, but now only depends on the difference $\mathsf{\Delta }t:=t_2-t_1$. The bimodality indicates a finite memory induced by the non-linearity. In order to further study this behavior, we first ensure the system has reached a stationary state by simulating from $t=0$ to $t=1500$ with a sample size of $8\times {10}^7$ for the cubic process. After setting $t_1=1500$, $t_2$ values are chosen from the data generated by further simulating the system for 300 time units. The PDF was approximated by a histogram with $60\times 60$ bins. The time steps were adaptively chosen by setting the local error tolerance as $0.01\sigma$.

Note that, because of the symmetry of diffusion term and the drift term, the bimodality of $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ for cubic SDE is always symmetric with respect to the line $x\left(t_1\right)=x\left(t_2\right)$. Therefore to make the numerical study of the bimodality easier, we restrict our attention to the diagonal of the joint probability density $P\left(x\left(t_1\right)=x\left(t_2\right)\right)$, which is then normalized to integrate to one.

To quantify the bimodality of $P\left(x\left(t_1\right)=x\left(t_2\right)\right)$, we use the following fitting function to perform a nonlinear fit and estimate the parameters.

$$P\left(x\left(t_1\right)=x\left(t_2\right)\right)=a\left[exp\left(-{\left(\frac{x}{c}\right)}^4\right)+bexp\left(-{\left(\frac{x-e}{d}\right)}^2\right)+bexp\left(-{\left(\frac{x+e}{d}\right)}^2\right)\right]$$

(22)

To motivate this fitting function, note that for $\mathsf{\Delta }t=0$, we get a purely quartic exponential, since it is nothing but the stationary state $P_s(x)$ of the cubic SDE. For $\mathsf{\Delta }t\to \infty$, the correlation between points reaches zero, $P\left(x\left(t_1\right),x\left(t_2\right)\right)\sim P_s(x_1(t_1))P_s(x_2(t_2))$, the product of two independent stationary distributions. The quadratic exponential terms are motivated by Gaussian distribution of the noise, and are found to describe the data accurately. There are two quadratic exponential terms, since the SDE is symmetric about the point $x=0$ and the bimodal peaks occur symmetrically on opposite sides.

In the fitting function, Parameter a is a measure of the overall height of the density curve. Parameter b represents the ratio of contribution from the quadratic exponential to the quartic exponential, and denotes the degree of bimodality. Parameter e is the location of the peaks. Parameter c is a measure of the overall spread of the density function, while parameter d is a measure of the spread of bimodal peaks. Note that a nonzero value for parameter e and b denotes bimodality in the distribution.

From Figure 10, we can see that the joint PDF becomes bimodal for a range of $\mathsf{\Delta }t$ values, since parameters b and e have nonzero values. The standard deviation of the quadratic term denoted by parameter d is almost a constant for a fixed $\sigma$, whereas the location of the bimodal peak denoted by parameter e changes slightly, but never becomes zero. That means that while transitioning from a bimodal to a unimodal distribution, the bimodal peaks do not continuously merge into one another, but slowly become less prominent and eventually disappear, as inferred from the value of parameter b. The artifacts towards the end of the curves of parameters b and e are due to the fact that it is not possible to consistently fit parameters b and e when there is negligible contribution from the quadratic exponential term as $b\to 0$.

The values of parameter a and c in Figure 10 can be understood by first considering the $\mathsf{\Delta }t\to 0$ and $\mathsf{\Delta }t\to \infty$ limits. For $\mathsf{\Delta }t\to 0$, we have $P\left(x\left(t_1\right),x\left(t_2\right)\right)=\delta \left(x\left(t_1\right)-x\left(t_2\right)\right)P_s(x)$, where $P_s(x)$ is the stationary state. The diagonal of joint PDF then becomes $P\left(x\left(t_1\right)=x\left(t_2\right)\right)=P_s(x):=a_{s}exp\left(-{(x/c_s)}^4\right)$. Here, after integrating x, we can find, $a_s\approx 0.55/c_s$ for all values of $\sigma$, which can be numerically verified. When $\mathsf{\Delta }t\to \infty$, since there is no correlation, $P\left(x\left(t_1\right),x\left(t_2\right)\right)=P_s\left(x\left(t_1\right)\right)P_s\left(x\left(t_2\right)\right)$. The normalized diagonal then becomes $P\left(x\left(t_1\right)=x\left(t_2\right)\right)\sim exp\left(-2{(x/c_s)}^4\right)$. The normalization factor can be derived by integrating out x. Therefore, when $\mathsf{\Delta }t\to \infty$, we have $a=2^{1/4}{a}_s$ and $c=c_s/2^{1/4}$, which agrees with numerical results. For intermediate $\mathsf{\Delta }t$ values, there will be correlation, and the behavior cannot be easily explained.

The trends in Figure 10 for different $\sigma$ values can be explained by looking at scaling relations. For the cubic equation, we have $d{x}_t=-\theta x_t^{3}dt+\sigma d{W}_t$. Looking at the individual terms, we can infer the dimensions, $\theta \sim x^{-2}{t}^{-1}$ and $\sigma \sim x^1{t}^{-1/2}$. Therefore, we expect the following scaling relations:

$$\begin{array}{cc}\hfill t& \sim \frac{1}{{\theta }^{1/2}\sigma }\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill x& \sim {\theta }^{-1/2}{\sigma }^{1/2}\hfill \end{array}$$

(24)

These relations are numerically verified in Figure 11 for a fixed value of $\theta$, for parameters a, c, d, e, and the relationship between noise and $\mathsf{\Delta }t$, corresponding to peak value of parameter b. Note that $a\sim {\sigma }^{0.5}$ in Figure 11, because parameter a is a normalization constant. $\int ap(x)dx=1\Rightarrow a\sim x^{-1}\sim {\theta }^{1/2}{\sigma }^{-1/2}$. The peak value of parameter b seems to be a constant with the value of $1/3$. It cannot be derived from simple scaling arguments alone; further investigation is required to understand its origin.

4.2. Evolution of Bimodality in the Non-Stationary State

In this section, we will look at how the bimodality in the joint PDF evolves qualitatively before $p(x,t)$ has reached a stationary state. Since bimodality occurs only for the cubic process, the following results are for the SDE $d{x}_t=-1.0x_t^{3}dt+0.1d{W}_t$. To this end, we follow the same prescription from Section 4, but with a slightly modified version of the fitting function (Equation (22)).

$$P\left(x\left(t_1\right)=x\left(t_2\right)\right)=aexp\left(-{\left(\frac{x}{c}\right)}^4\right)+bexp\left(-{\left(\frac{x-e}{d}\right)}^2\right)+bexp\left(-{\left(\frac{x+e}{d}\right)}^2\right)$$

(25)

This modification is made since for some $t_1$ and $t_2$ values, there is no contribution from the quartic exponential term, unlike in the case of stationary state, where there is always a quartic exponential contribution to the distribution.

The nonstationary phase exhibits (Figure 12) rich behavior, which asymptotically transitions to the stationary state behavior as $t_1$ becomes large. For small $t_1$ values, Parameter a ≪ Parameter b, since $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ has very little contribution from the quartic exponential term. This is because the initial distribution is a Gaussian distribution, and it is slowly evolving towards the quartic exponential distribution in the stationary state. For larger $t_1$ values, Parameter a starts dominating, since $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ has predominant contribution from quartic exponential term, as expected of the system reaching its quartic exponential stationary state. For some intermediate $t_1$ values, depending on $\mathsf{\Delta }t$ values, we see an interesting behavior where $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ becomes purely a quartic exponential (Parameter b = 0) twice before becoming a mixture of quadratic and quartic exponential terms. Further investigation into this behavior is not undertaken at this time, and only serves to demonstrate the potential capabilities of GPU-accelerated MC SDE simulation for future work.

5. Discussion

In this work, we developed a method for fast and accurate study of the Information Geometry of SDEs using Monte Carlo simulation. We identified the computational challenges and overcame them by using GPU computing. Specifically, the major limitation with MC SDE simulation in the estimation of $\mathsf{\Gamma }$ was the sub-optimal overlap of PDFs at subsequent time points. We solved this problem by developing an interpolation method to compute PDFs with optimal overlap.

As an application of the new method, we compared the Information Geometry of SDEs with a linear and a cubic damping force. We were able to reproduce the analytic results for the linear SDE and the previous numerical results [34] for cubic SDE obtained using FPE solvers. This was particularly true of large values of $x_0$, ${\mathcal{L}}_{\infty }\sim x_0$ for the linear case and ${\mathcal{L}}_{\infty }\sim x_0^m$ for the cubic case, where $m=1.88$ when $\theta =1.0$ and $\sigma =0.1$. We further showcased the advantage of MC SDE simulation over FPE solver by computing the joint PDF $P\left(x\left(t_1\right),x\left(t_2\right)\right)$. Unlike the linear SDE, the cubic SDE led to an interesting bimodal PDF $P\left(x\left(t_1\right),x\left(t_2\right)\right)$, which was observed even after reaching a stationary state. After reaching the stationary state $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ only depends on $\mathsf{\Delta }t=t_2-t_1$. In the stationary state, we further studied the bimodality by quantifying it and looking at the power-law scaling relations with respect to noise levels $\sigma$ and provided theoretical scaling arguments. The maximum value of the ratio of quadratic to quartic contribution (Parameter b) is found to be a constant $1/3$ irrespective of the noise levels, which requires further study. Finally, we qualitatively looked at $P\left(x\left(t_1\right),x\left(t_2\right)\right)$ in the nonstationary state. The MC SDE simulation can be an important tool for further studying this behavior.

It is important to note that the methods that we developed here for one stochastic variable are general, and can be extended for more than one variable, as well as for investigating different metrics or thermodynamic quantities. These will be addressed in future work. Furthermore, it will be of interest to investigate fast implementations of Kernel Density Estimators [48,49] which will improve the accuracy of joint PDF estimates. We note that in numerical experiments, it was seen that compared with histograms, using Kernel Density Estimators to estimate PDFs provided 2–5 times reduction in error with identical PDFs, albeit with a performance trade-off.